[CIG-SHORT] Elastoplastic material question
Brad Aagaard
baagaard at usgs.gov
Mon Nov 25 07:57:06 PST 2013
Eric,
The zz component of the initial stress is given via a spatial database
just like all other potentially spatially varying parameters. In this
case we include it in the spatial database for the initial state.
Here are example settings for specifying spatial databases for the
initial stress and initial state.
[pylithapp.timedependent.materials.elastic]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress
db_initial_stress.iohandler.filename = initial_stress.spatialdb
db_initial_stress.query_type = linear
db_initial_state = spatialdata.spatialdb.SimpleDB
db_initial_state.label = Initial state variables
db_initial_state.iohandler.filename = initial_state.spatialdb
db_initial_state.query_type = linear
A spatial database for the initial state:
#SPATIAL.ascii 1
SimpleDB {
num-values = 5
value-names = stress-zz-initial plastic-strain-xx
plastic-strain-yy plastic-strain-xy plastic-strain-zz
value-units = MPa None None None None
num-locs = 5
data-dim = 1
space-dim = 2
cs-data = cartesian {
to-meters = 1.0e+3
space-dim = 2
}
}
0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 -11.950 -134.34391 0.0 0.0 0.0 0.0
0.0 -11.951 -166.72940 0.0 0.0 0.0 0.0
0.0 -11.952 -199.12032 0.0 0.0 0.0 0.0
0.0 -80.000 -1332.8 0.0 0.0 0.0 0.0
Brad
On 11/24/2013 03:26 PM, Eric Lindsey wrote:
> I have an idea where the problem is now - it seems I need to specify the
> initial stress in the ZZ direction. This isn't necessary for the elastic
> case, because it can be solved for explicitly, but in the plastic case I'm
> getting significant strain in the ZZ direction (third component of the
> 4-component "plastic_strain" property, as I understand), even if no shear
> is applied. I'm guessing that the default assumption is a zero sigma_ZZ
> stress component -- in combination with the X and Y compression I've
> applied in the initial stresses, this leads to problems. I didn't think of
> checking this earlier, I assumed "plane strain" meant it would constrain
> this component of strain to be zero!
>
> There is an option for the 2D viscoelastic materials called
> "stress_zz_initial" but it's not documented anywhere except a brief mention
> in the manual; it's not used in any examples I can find. Is this also an
> option for the DruckerPragerPlaneStrain material, and how do I implement
> it? I've tried:
>
> [pylithapp.timedependent.materials.crust]
> stress_zz_initial = -0.2e+06*Pa
>
> but this gives the error:
>
> -- unrecognized property
> 'timedependent.materials.druckerpragerplanestrain.stress_zz_initial'
>
> Am I on the right track?
>
> Thanks,
> Eric
>
>
> On Sat, Nov 23, 2013 at 7:22 PM, Charles Williams <willic3 at gmail.com> wrote:
>
>> 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
>>
>>
>> _______________________________________________
>> 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
>
More information about the CIG-SHORT
mailing list