Hello Brad,
I am trying to simulate single step elastic deformation with gravity on. I have applied Fault normal tractions along the fault nodeset and all other faces (except the top one) are constrained using roller dirichlet BC.
For some reason I am getting a unusual shear stress distribution at the bottom node in the fault nodeset. Please see the attached fig. My version of Pylith is 2.2.2
Thank you
Some things to check:
If you still need help after checking these things, please supply a diagram/cartoon/sketch showing the domain with the fault and BC labeled.
Thanks for the swift response, Brad.
As you told, I have checked back my inputs and still the results are same. As you can see from the figure below, for some reason my fault normal traction is not balancing out the gravity effect. So, as I run the simulation with gravity on, the tip of the geometry bends downwards.
Could you please explain what am I doing wrong here?
Thank you
I don’t quite understand what you are trying to do. Is the wireframe in the figure you provided the entire domain or just a piece of it? Why are you trying to balance fault normal tractions with gravity?
If the wireframe is your domain, then I think the problem is your boundary conditions along the “fault” and bottom boundary. I assume your “fault” is a Neumann boundary condition, because an actual fault on the boundary of the domain doesn’t make sense numerically. If you don’t have initial stresses in the bulk material, then I don’t think your boundary conditions are consistent and will balance the gravitational body forces; I would expect to see the stress concentration at the intersection of the “fault” and bottom Dirichlet boundary condition. In this case, I think the problem is trying to approximate a fault with a Neumann boundary condition.
Normally, when we use gravity and initial stresses, we use a domain that is a box and impose Dirichlet boundary conditions that constrain the degree(s) of freedom perpendicular to the boundary on the lateral sides and bottom of the domain. Inside the domain we impose initial stresses in the bulk materials equal to the lithostatic stress. If we have a fault, we do not impose initial fault normal tractions because the gravitational body forces will give rise to fault tractions equal to the lithostatic load.
Hi Brad,
I am trying to compute the effect of gravity on the forearc region as depicted by the wireframe figure. I am not using a footwall block, instead am trying to balance gravity out by applying normal traction on the fault interface. It is same as having a surface along the fault that is constrained in the Z-direction with gravity on(not sure how this can be done in PyLith).
At present, with fault normal traction equal to z rho g (with gravity on), the model bends downwards (shown above). As I increase the normal traction at some point it starts bending upwards. I am not sure that why after applying gravity the normal traction is not getting properly balanced, as if a rigid plate interface is existing along with the fault that is frictionless. My shear traction is zero along the ‘fault’ at present, meaning that fault friction is zero in this case.
But in next case, I do want to model a fault with some friction coefficient, which I am planning to do by imposing shear traction along the fault interface.
But before that I need your help to get this first case modeled properly, which clearly is not correct at the moment.
Thank you for the help.
If you don’t impose initial stresses, then turning on gravity with Dirichlet boundary conditions will not give rise to \sigma_xx = \sigma_yy = \sigma_zz = \rho g z. The horizontal stresses (\sigma_xx and \sigma_yy) will differ from \sigma_zz due to Poisson’s ratio due to vertical displacement. You can easily verify this by solving the boundary value problem analytically. If you want \sigma_xx = \sigma_yy = \sigma_zz = \rho g z then you need to impose initial stress equal to \rho g z. In this case, you can impose a Neumann boundary condition along your “fault” boundary equal to \rho g z and you should get zero deformation.
So, I reran my simulations with initial stresses sigma_xx, sigma_yy, sigma_zz as rho g z, with gravity on and with normal traction applied along the fault boundary again equal to rho g z. (g =9.81 m/s2 at all places)
But again my domain is deflecting downwards indicating that the forces are not getting balanced. I have checked my input values multiple times and there doesn’t seem to be a problem there.
Thanks again.
I suggest checking what axial stresses are in the solution. Are they equilibrating to the initial stresses (density * g * z) or something else? Do you get the desired behavior if you use Dirichlet BCs on the bottom and lateral boundaries while maintaining a free surface on the top boundary? If so, then the Neumann BC is likely the issue. Is the normal traction on the Neumann BC compressive?
With such deformation (shown by domain) we can’t expect the axial stresses to give reasonable values (shown below is sigma_zz)?
The normal traction is negative (compressive) normal along fault boundary (AB from figure).
Initial stresses applied along A->C, negative (sigma xx, yy, zz) and (shear stresses = zero).
Sorry if I missed catching this earlier. What is the geometry of your domain? Your sketches only show the xz plane. If it is 2D, then it needs to be in the xy plane and you would impose gravity in the -y direction instead of the -z direction. You would still impose initial stress \sigma_xx = \sigma_yy = \sigma_zz = density * g * y to get the correct out of plane (sigma_zz) stress.
No, my geometry is 3-d, with the y-direction domains (positive and negative y) having Dirichlet BC constraining them in y.
What happens if you replace the Neumann BC on boundary AB with a Dirichlet BC constraining x or z? If you have the correct initial stresses, then they should balance the gravitational body forces and there should be zero deformation.
There sigma_zz stress looks like it is following densitygz on the right side of the domain (+x side). However, sigma_zz on the left side (-x side) looks quite different. There seems to be a horizontal discontinuity near the middle of the domain.
Some possible issues to investigate:
We might be able to provide more guidance if you provide a more detailed sketch/diagram of your 3D domain that includes information about materials, all boundary conditions, etc.
The density is uniform as I am trying to solve a simple problem to begin with. It is possible that the mesh quality can cause weird deformation but still it is difficult to comprehend that it is causing 4-5 orders of magnitude difference to the actual solution.
Since, as this point it is becoming difficult to explain and discuss with figures, I am sharing my pylith input files here, so that you can have a look at them.
Thank again.
SimulationFiles.zip (41.9 KB)
Thanks for providing the input files. There are a few issues that I found when I examined the input files.
Attached is a set of simulation inputs for two 2D problems. One uses a box and one uses a geometry similar to yours but with a flat top surface. In both cases the initial stresses and Neumann BC exactly match the gravitational body forces, so there is no deformation and the solution converges immediately with the expected stress field.
GravityNeumann.zip (12.8 KB)
Hi Brad, thanks for looking into my input files. I am planning to analyze the current problem in the y-direction in future work for the same model geometry, therefore I am having a thin section in that direction. For this problem, however I assume that I should get the same results as if it was in 2-D.
So, it is important for me to have the same geometry at the top surface. As you have explained, it seems that I need to first compute initial stresses with the bottom and slope constrained in vertical direction (step 1) and then use those computed stress values (as initial stress input) in the following step (step 2).
Will these stress values from step 1 be a reasonable approximation to be used in step 2 for this problem?
Also, I use cubit to extract coordinates of node points (which come out as float values) for specifying the input stress values. What if I miss some significant decimal points and end up in between 1st and 2nd nodes instead of the 1st one. For example 1st node (0,1.00) and 2nd node (0,1.02) and I end up specifying stress at position (0,1.01), node at which position does not exist.
How will Pylith computations take care of such situation internally?
Thanks again for the help.
Initial stresses computed from a static simulation can be a reasonable approximation, if the problem is setup with the appropriate boundary conditions. You will likely need to experiment to get things right.
The spatial database implementations try to take into account floating point roundoff errors, so that coordinates do not need to match exactly. The safest thing to do is to move points slightly outside the boundaries or add an extra layer of points outside the boundaries to ensure all points in the domain lie within the set of points in the spatial database.
Hi Brad, your suggestions helped me get zero deformation, with Dirichlet BC (DBC) in vertical direction along A->B->C (shown below step 1).
Hi Brad, I am having an issue with my simulation, which is taking unusually longer time to run.
Problem description: I have one fault (along AB) and am trying to see the elastic response of my system with application of gravity, similar to my other post but this one has a plate underneath.
I have tried both Implicit as well as ImplicitLgDeform problem formulations for my model. Both are taking similar runtime (~30 min) with Intel 11th gen i7 8 core CPU overclocking most of the time, the results being produced are also incorrect.
Additionally, in many instances the simulation runs faster but ends up giving zero deformation (fig below), which is not making sense to me.
Here is the link to my json file.
Please let me know if you need any other information.
Thank you
Viven
