Hello,
I am using PyLith 3.0.1 to model a sill intrusion in the upper mantle. I have a 2D domain crust, mantle 80 km by 60 km depth. The intrusion is located at the center of the mesh at a depth of 40 km and the size is 10 km by 500 m. I am editing the reverse-2d/reference_state example to model this problem.
I want to look at the stress difference the intrusion generates due to its lower density.
I first made the intrusion have the same density as the mantle and assigned the same reference stresses as the mantle. However, I am generating stresses around the intrusion even though it has the same density as the mantle. I would expect stress differences around the intrusion IF it had a density that is different to the mantle.
Why am I generating stresses around the intrusion? Could it be that I am assigning the wrong reference stresses to the mantle/intrusion? Attached you will find the files and below you will see domain displacement, stress_xx, and stress_yy.
sill_intrusion.zip (1.3 MB)
Thank you,
Eddie
Please provide a diagram of the boundary value problem you are trying to solve with boundary conditions, materials, etc. labeled. How are you modeling the intrusion?
Hello Brad,
This is a static model of an intrusion. I want to look at the stresses the intrusion generates due to its density difference.
I found the error. When you create the intrusion curve loop in the Gmsh Python script generate_gmsh.py, you are traversing the curves clockwise, which is inconsistent with the counter-clockwise direction for the other curve loops. This also explains why the cells in the intrusion show up very dark in Gmsh (you are looking at the “back” of the cells). If you fix the script and traverse the curves in a counter-clockwise direction for all of the curve loops, the reference stress for uniform material properties exactly balances the gravitational body forces. I didn’t test the more complicated case of nonuniform material properties.
I don’t know why Gmsh didn’t reject your curve loop that had the incorrect orientation.
Note: PyLith uses 8.60665 m/s for the acceleration for gravity. You can change this in the properties of the GravityField object if you want.
One more thing, you should update from PyLith v3.0.1 to v3.0.3 as there are lots of small bugs that have been fixed with the subsequent releases.
Wow! Thank you! It is weird that Gmsh did not reject my loops. Okay, I will update to v3.0.3 thank you Brad!
We have added an additional topology check in the current PyLith main branch that will catch inverted cells that occur in cases like this. These topology checks occur immediately after the mesh is imported. It will be available (and turned on by default) in the next PyLith release.
