[CIG-SHORT] Elastoplastic material question

Sat Nov 23 19:22:07 PST 2013

Hi Eric,

The Neumann boundary conditions should work the same no matter what material you’re using.  I’m assuming that plastic yield is occurring prior to the strange stress results.  One thing to try might be varying the dilatation angle.  When it is equal to the yield angle, you can get unrealistic amounts of dilatation.  A value of 0 would mean no dilatation (incompressible).  If you do this, you will probably need to change some solver settings because the plasticity is no longer associated, and your stiffness matrix is no longer symmetric.

Cheers,
Charles

On 23/11/2013, at 11:54 am, Eric Lindsey <elindsey at ucsd.edu> wrote:

> Hi Brad,
> 
> Thanks a lot for looking at this; after following your suggestions I'm
> still having a problem though.
> 
> First, sorry about the strange boundary conditions on the sides -- I
> think I originally did have the correct conditions for simple shear,
> but took them out while testing to see what was causing the issue, and
> just never put them back. I've fixed them up as you suggest, so that
> for the elastic case there is a constant value of each stress
> component throughout the domain at every timestep, and it never goes
> into absolute tension at any location. I've also used the no-fault
> petsc solver settings as suggested.
> 
> However, I still get the same "Infeasible stress state" error when
> running the plastic case, unless I allow for tensile yield. In either
> case the solution looks strange, with significant spatial variations
> in the stress, even before any shear is applied. The problem starts at
> the corners, which are going into absolute tension because the sides
> have extended outwards, despite the normal-traction conditions which
> prevented this in the elastic case.
> 
> Is there some difference in how the Neumann boundary conditions need
> to be applied for the DruckerPrager material? I'm wondering if they
> are being interpreted differently in the two cases -- or perhaps I
> have a problem with the solver, but I don't see any indication of that
> in the log; convergence is fine up until I get the error message.
> 
> Thanks,
> Eric
> 
> On Tue, Nov 19, 2013 at 4:22 PM, Brad Aagaard <baagaard at usgs.gov> wrote:
>> 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
>>> 
>> 
>> _______________________________________________
>> CIG-SHORT mailing list
>> CIG-SHORT at geodynamics.org
>> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
> <log_plastic.txt><pylithapp.cfg><plastic.cfg><elastic.cfg>_______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short

Charles A. Williams
Scientist
GNS Science
1 Fairway Drive, Avalon
PO Box 30368
Lower Hutt  5040
New Zealand
ph (office): 0064-4570-4566
fax (office): 0064-4570-4600
C.Williams at gns.cri.nz

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://geodynamics.org/pipermail/cig-short/attachments/20131124/6c0f841b/attachment.html>