Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[GeoMechanicsApplication] Added a 3D dynamic test case #12656

Merged
merged 6 commits into from
Sep 6, 2024

Conversation

mnabideltares
Copy link
Contributor

@mnabideltares mnabideltares commented Sep 2, 2024

📝 Description
A 3D dynamic test case is added to the integration tests.

🆕 Changelog

  • Added a test case
  • Added documentation

Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Dear Mohammed,

Thank you for the elaborate work on testing here. Even if we are not miraculously close to the semi-analytical formulation I think these tests are valuable.

Regards, Wijtze Pieter

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The direction of the force seems really vertical, while the block looks to be at an angle with the vertical. Force arrow should be aligned with the edge of the block.

"body_domain_sub_model_part_list" : ["soil_block"]
},
"output_processes": {
"gid_output" : [{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if this output should be asked for when only the next output process gives the things that are actively checked.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

- $\nu$ = Poisson's ratio $\mathrm{[-]}$
- $\rho$ = Mass density $\mathrm{[kg/m^3]}$

In this test case, $E = 91800 \mathrm{Pa}$ and $\nu = 0.25$. This leads to shear and compression wave velocities of $6 \mathrm{m/s}$ and $10.39 \mathrm{m/s}$, respectively. Hence, in order to avoid the effects of reflecting waves from the boundaries, considering a domain of $10 \mathrm{m} × 10 \mathrm{m} × 10 \mathrm{m}$ with a simulation time of 1 second is sufficient.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With 10.39 m/s and a distance of 10m, the wave starts to reflect at 10/10.39 s which is less than the computed 1 s.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, there may happen a slight reflection at the end, but as it is just a fraction of time, we don't observe it. The reflection may affect the very end points of the figure.
Tried 20m domain, but the results did not change much around 1s

## Case Specification
In this test case, a three-dimentional soil block with dimensions of 10 m x 10 m x 10 m is considered. This is a quart of the full domain, as this case is axisymmetric. The pressure is fixed at 0 Pa throughout the entire simulation to eliminate the effect of water pressure. This is done to allow comparison of our results with published semi-analytical solutions. The bottom side of the domain is fixed, while the vertical sides are allowed to move only in the tangential directions. A sudden force of -250 N is then applied in the vertical direction at the top surface of the block. The simulation spans 1 second. This test is conducted on a mesh with triangular elements of type 3D4N. The deformation at the surface is then compared with semi-analytical results.

The geometry and boundary conditions are shown below:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the Dirichlet boundary conditions are not in the picture ( and would be a lot of work to show clearly ). Therefore I propose to only mention the geometry and loading conditions. Not the boundary conditions.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

- $\nu$ = Poisson's ratio $\mathrm{[-]}$
- $\rho$ = Mass density $\mathrm{[kg/m^3]}$

In this test case, $E = 91800 \mathrm{Pa}$ and $\nu = 0.25$. This leads to shear and compression wave velocities of $6 \mathrm{m/s}$ and $10.39 \mathrm{m/s}$, respectively. Hence, in order to avoid the effects of reflecting waves from the boundaries, considering a domain of $10 \mathrm{m} × 10 \mathrm{m} × 10 \mathrm{m}$ with a simulation time of 1 second is sufficient.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

91800 seems a strange number in [Pa] I would expect something in the order of 80 MPa.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do agree that a such number does not make much sense practically. We used it just to compare.

**Source files:** [Dynamic solution in 3D](https://github.com/KratosMultiphysics/Kratos/tree/master/applications/GeoMechanicsApplication/tests/test_dynamic/test_constant_point_load_3d)

## Case Specification
In this test case, a three-dimentional soil block with dimensions of 10 m x 10 m x 10 m is considered. This is a quart of the full domain, as this case is axisymmetric. The pressure is fixed at 0 Pa throughout the entire simulation to eliminate the effect of water pressure. This is done to allow comparison of our results with published semi-analytical solutions. The bottom side of the domain is fixed, while the vertical sides are allowed to move only in the tangential directions. A sudden force of -250 N is then applied in the vertical direction at the top surface of the block. The simulation spans 1 second. This test is conducted on a mesh with triangular elements of type 3D4N. The deformation at the surface is then compared with semi-analytical results.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

triangular elements are 2D. Here tetrahedral elements are used.

$$ M = \frac{E \left(1 - \nu \right)}{\left( 1 + \nu \right) \left(1 - 2 \nu \right)} $$

- $G$ = Shear modulus $\mathrm{[Pa]}$
- $M$ = P-wave modulus $\mathrm{[Pa]}$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm guessing that the usual name for this is bulk modulus ( mostly denoted as K ) or is it one of the 2 constants of Lame that usually have Greek symbols lambda and mu? If we are following the notation used by Verruit.we shouldn't change it.


In the semi-analytical solution, there is a singulair point around $\tau = 1.1$. However, this behavior, which is caused by the arrival of the Rayleigh wave, is captured by a peak in the numerical solution.

Note: To avoid numerical oscillations, the force is gradually applied over a time span of 0.1 seconds.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

On a total time of 1 s the 0.1 for the increase seems rather long. That may have an influence on what you describe as shift in time for the peak behavior. Is the same ramping time used in the 2D computation?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the case of < 0.1, the results even shift more to the left side. for the case of > 0.1, the peak becomes smoother but not much shifts to the right

Copy link
Contributor

@rfaasse rfaasse left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for adding this test case with a clear readme! I left some suggestions/questions, let me know what you think!

"output_file_name" : "json_output.json",
"output_variables" : ["DISPLACEMENT_Y"],
"gauss_points_output_variables": [],
"time_frequency" : 0.0099999999
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems very specific, what is the reason it can't be 0.01?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the case of 0.01, it somehow skips some of the outputs (don't know the reason). In order to be certain for a correct output, we make it slightly less than 0.01.

Copy link
Contributor Author

@mnabideltares mnabideltares left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@WPK4FEM and @rfaasse
Thank you for review. These are my changes and comments

"output_file_name" : "json_output.json",
"output_variables" : ["DISPLACEMENT_Y"],
"gauss_points_output_variables": [],
"time_frequency" : 0.0099999999
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the case of 0.01, it somehow skips some of the outputs (don't know the reason). In order to be certain for a correct output, we make it slightly less than 0.01.

- $\nu$ = Poisson's ratio $\mathrm{[-]}$
- $\rho$ = Mass density $\mathrm{[kg/m^3]}$

In this test case, $E = 91800 \mathrm{Pa}$ and $\nu = 0.25$. This leads to shear and compression wave velocities of $6 \mathrm{m/s}$ and $10.39 \mathrm{m/s}$, respectively. Hence, in order to avoid the effects of reflecting waves from the boundaries, considering a domain of $10 \mathrm{m} × 10 \mathrm{m} × 10 \mathrm{m}$ with a simulation time of 1 second is sufficient.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, there may happen a slight reflection at the end, but as it is just a fraction of time, we don't observe it. The reflection may affect the very end points of the figure.
Tried 20m domain, but the results did not change much around 1s

- $\nu$ = Poisson's ratio $\mathrm{[-]}$
- $\rho$ = Mass density $\mathrm{[kg/m^3]}$

In this test case, $E = 91800 \mathrm{Pa}$ and $\nu = 0.25$. This leads to shear and compression wave velocities of $6 \mathrm{m/s}$ and $10.39 \mathrm{m/s}$, respectively. Hence, in order to avoid the effects of reflecting waves from the boundaries, considering a domain of $10 \mathrm{m} × 10 \mathrm{m} × 10 \mathrm{m}$ with a simulation time of 1 second is sufficient.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do agree that a such number does not make much sense practically. We used it just to compare.


In the semi-analytical solution, there is a singulair point around $\tau = 1.1$. However, this behavior, which is caused by the arrival of the Rayleigh wave, is captured by a peak in the numerical solution.

Note: To avoid numerical oscillations, the force is gradually applied over a time span of 0.1 seconds.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the case of < 0.1, the results even shift more to the left side. for the case of > 0.1, the peak becomes smoother but not much shifts to the right

## Case Specification
In this test case, a three-dimentional soil block with dimensions of 10 m x 10 m x 10 m is considered. This is a quart of the full domain, as this case is axisymmetric. The pressure is fixed at 0 Pa throughout the entire simulation to eliminate the effect of water pressure. This is done to allow comparison of our results with published semi-analytical solutions. The bottom side of the domain is fixed, while the vertical sides are allowed to move only in the tangential directions. A sudden force of -250 N is then applied in the vertical direction at the top surface of the block. The simulation spans 1 second. This test is conducted on a mesh with triangular elements of type 3D4N. The deformation at the surface is then compared with semi-analytical results.

The geometry and boundary conditions are shown below:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

rfaasse
rfaasse previously approved these changes Sep 6, 2024
Copy link
Contributor

@rfaasse rfaasse left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for incorporating the changes and diving into this problem in the first place! From my perspective this is good to go!

Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you Mohamed,
Almost there. I'm afraid that I caused some confusion about bulk modulus and p wave modulus.

Did we try the computation with both Rayleigh parameters 0? Verruijt's equations have no damping and I'm not sure which frequencies we are damping here. Adding damping usually smears out the peak behavior.

Wijtze Pieter

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now it is consistent, thank you.

$$ M = \frac{E \left(1 - \nu \right)}{\left( 1 + \nu \right) \left(1 - 2 \nu \right)} $$

- $G$ = Shear modulus $\mathrm{[Pa]}$
- $M$ = Bulk modulus $\mathrm{[Pa]}$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

M is the P wave modulus. See e.g. https://en.wikipedia.org/wiki/Lam%C3%A9_parameters that is different from the bulk modulus usually denoted with K


$$ c_p = \sqrt{\frac{M}{\rho}} $$

where $G$ and $M$ are shear modulus and bulk modulus, respectively. They are defined as:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

bulk modulus -> P wave modulus.

Copy link
Contributor

@WPK4FEM WPK4FEM left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for adding the test and documentation.

@rfaasse rfaasse merged commit 4ef6e05 into master Sep 6, 2024
9 of 11 checks passed
@rfaasse rfaasse deleted the geo/12653-3d-dynamic-test branch September 6, 2024 15:22
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

Successfully merging this pull request may close these issues.

[GeoMechanicsApplication] The 3D dynamic test case results do not match the analytical solution
3 participants