[CIG-SHORT] Elastoplastic material question

[Brad Aagaard]

Eric,

When I run the elastic case, I see the xx and yy (axial) components of 
the stress tensor becoming tensile near the corners as the +y (top) 
boundary moves to the right (+x direction). The domain is not under pure 
shear because there are no shear tractions on the +x and -x sides of the 
domain (the shear stress is zero on those boundaries); there is 
deformation in the y-direction on the +x and -x sides of the domain.

In the plastic case the areas with the tensile stresses are where the 
plastic strain first appears. The plastic strain soon afterwards shows 
up in three bands across the middle of the domain, where the shear 
strain is largest. These bands are one cell wide, suggesting the width 
is not well resolved. As the shear strain increases, the bands grow. 
This appears consistent with your comments. I don't see any indication 
that there are bugs in the code or that the solver is not working as 
intended.

Also note that in a nearly uniform resolution mesh with slight 
irregularities associated with the triangular discretization, there will 
be small spatial variations in the accuracy of the solution that can 
lead to localization of the yielding.

To generate pure shear, I recommend
u_x = u_y = 0 on -y boundary
u_y = 0 on +y boundary
T_shear = f(t) on -x, +x, and +y boundaries

For PETSc solver settings, when you don't have a fault you should use

[pylithapp.timedependent.formulation]
matrix_type = aij

[pylithapp.petsc]
pc_type = ml

When you have a fault, use

[pylithapp.timedependent.formulation]
split_fields = True
use_custom_constraint_pc = True
matrix_type = aij

[pylithapp.petsc]
fs_pc_type = fieldsplit
fs_pc_use_amat = true
fs_pc_fieldsplit_type = multiplicative
fs_fieldsplit_0_pc_type = ml
fs_fieldsplit_1_pc_type = jacobi
fs_fieldsplit_0_ksp_type = preonly
fs_fieldsplit_1_ksp_type = preonly

For more info on the solver parameters, you may want to review Matt's 
presentation and slides from the June online tutorial 
(http://www.geodynamics.org/cig/community/workinggroups/short/workshops/cdm2013/agenda). 
There is also a comparison of some combinations of solver parameters in 
our JGR paper (http://dx.doi.org/10.1002/jgrb.50217).

Regards,
Brad


On 11/18/2013 05:01 PM, Eric Lindsey wrote:
> I'm having trouble understanding the DruckerPragerPlaneStrain material. I'm
> using a simple homogeneous domain with no faults, and a simple Dirichlet BC
> imposing shear on the top/bottom. I've imposed an initial isotropic
> compression, then I add the shear displacement gradually until it should
> exceed the yield stress. On the sides I have a Neumann condition to
> maintain the normal stress; for the moment I just set the shear tractions
> to zero, but this value doesn't affect the results I'm getting. I expected
> to see uniform plastic shear throughout the domain, but instead get the
> message:
>
> RuntimeError: Infeasible stress state - cannot project back to yield
> surface.
>
> If I include "allow_tensile_yield = True", the model runs, but instead of
> homogeneous strain I'm getting the attached result. This is strange because
> at no point should the material be under absolute tension; the maximum
> shear stress in the elastic case is never larger than the magnitude of
> compressive stress. I think I am misunderstanding the setup of this
> material somehow; hopefully it's a simple error? Input files are attached,
> the elastic case works just fine. Any suggestions would be much appreciated.
>
> Thanks,
> Eric
>
> Relevant lines from the spatial databases:
>
>
>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>