[CIG-SHORT] Bizarre tractions on dynamic fault with gravity

[Brad Aagaard]

Romain,

I believe I understand the source of the problem you are getting with 
gravity, initial stresses, and a fault with friction.

I plotted the fault normal traction at time step 0 and compared it 
against rho*g*(y-500); the 500 comes from your additional overburden 
imposed at the top of the domain. I see a residual that has a depth 
dependence with outliers at the top and bottom. A similar plot for the 
stress in the cells shows a match to rho*g*(y-500).

The stress field within the material has a linear variation with depth. 
The stresses get projected onto the fault via the finite-element 
formulation. The fault vertices in the middle of the domain have cells 
above and below, so they get some average of the stress in the cells 
that are above and below the fault vertices. As a result, they end up 
quite close to -rho*g*(y-500), where y is the coordinate of the fault 
vertex. The vertices at the top and bottom boundary have cells only 
below or above them, so the projection from the finite-element 
formulation comes from only below or above the vertex, so it has a 
greater residual. Decreasing the size of the cells will reduce the residual.

An additional note is that linear triangular cells are constant 
strain/stress cells, so the stresses/strains within a cell are uniform. 
As a result, they cannot represent a linearly varying stress field as 
well as quadrilateral cells or higher order triangular cells. Thus, for 
a similar resolution mesh, quad cells may reduce the residuals with 
respect to rho*g*y. We don't yet have higher order basis functions 
implemented in PyLith (this is high on our TODO list).

Brad


On 11/11/2014 10:26 AM, Romain Jolivet wrote:
> Dear Brad et al.,
>
> I am running a 2-D problem with an elastic rectangle (200km*50km) cut by a fault to simulate a thrust fault.
> Gravity is on in this problem. For a picture, please look at the file geometry.png.
>
> Because gravity is on, I have to prescribe initial stress conditions (isotropic stress tensor everywhere = rho*g*depth).
> For the same reason, I prescribe a Neumann condition on the bottom boundary with sigma_n = rho*g*depth and tau = 0 Pa
>
> The fault is a RateStateAgeing with slip strengthening properties. I push on the right side (dirichlet condition, 1 cm/year) and I want to see it sliding all the way, steady. The other side is fixed.
> Because I do not want to wait for too many time steps, I prescribe an initial shear stress along the fault equal to the normal stress multiplied by the standard friction coefficient (normal stress which should be equal to rho*g*depth).
> Figure initialTraction.png shows the prescribed traction on the fault (x-axis is depth in m, 0 m is the surface).
>
> Still, when I run the problem, it turns out that the pre-step computes a normal traction on the fault that is not equal to rho*g*depth, at least not everywhere. It seems fine for all the nodes inside the material, but the two nodes on the edges are off the linear trend where they should be (see afterElasticPrestep.png). I might be doing something wrong, but I wonder how tractions are estimated on the edges of the fault?
> This case leads to non-constant slip along the fault, especially for the first time step.
>
> I could compute the pre-step and then feed the tractions as a spatialdb, but I feel this problem is too simple to require that and there should be a solution.
>
> I joined a tar ball with the problem configure files, spatialdb, etc, if you want to try it (If you feel like numbers for stress are a bit awkward, it is because I want to do as if the whole thing was buried under 500 m of material. Removing that 500m of excess does not change anything in my case.).
> Let me know what you think,
> Romain
>
> —————————————————————————————————————
> —————————————————————————————————————	
> !!!!!! ADDRESS CHANGED !!!!!
>
> Romain Jolivet
> Postdoctoral Fellow
>
> University of Cambridge
> Department of Earth Sciences
> Bullard Labs
> Madingley Rise
> Madingley Road
> Cambridge CB3 0EZ
> United Kingdom
>
> email: rpj29 at cam.ac.uk
> Phone: +44 1223 748 938
> Mobile: +44 7596 703 148
>
> France: +33 6 52 91 76 39
> US: +1 (626) 560 6356
> —————————————————————————————————————
> —————————————————————————————————————
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>