[CIG-SHORT] Elastoplastic material question
Eric Lindsey
elindsey at ucsd.edu
Mon Nov 25 12:56:05 PST 2013
Thanks for the help!
Eric
On Mon, Nov 25, 2013 at 7:57 AM, Brad Aagaard <baagaard at usgs.gov> wrote:
> 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
>>
>>
> _______________________________________________
> CIG-SHORT mailing list
> CIG-SHORT at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-short
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://geodynamics.org/pipermail/cig-short/attachments/20131125/d4a8c4f4/attachment-0001.html>
More information about the CIG-SHORT
mailing list