simple shear results

[krikorik]

Hi everybody
I would like to ask if anybody knows why of the difference in results on this simple test.

It is a 2D elastic analysis of a square of side of length 1 solicited to simple shear.

The material is linear elastic and isotropic with Young=200E5, and Poisson=0.3

The final deformation gradient is F = [1, 1, 0; 0, 1, 0; 0, 0, 1]

I used either numerical and analytical methods to compute the component of stress and Von Mises equivalent stress

For the numerical I used abaqus implict and explicit. Also you can enter the material constant as linear elastic isotropic or as anisotropic entering the corresponding 21 constants of the isotropic material. In teh last case you are forced to use the *orientation keyword. When orientation was used I set the local axis coincident with the global axis, that is 1,0,0 ; 0,1,0 orientation

Also due to the large deformation I used *Nlgeom=yes

The element is a four node reduced integration CPE4R that can eb used either in explicit or implicit analysis

I run in double precision and make sure that the inertia effects were negligible wrt the elastic response.

************************************************** ***********
Here are the results I got for 11 different ways to present the same problem.


Case| Mises| S11| S22| S33| S12
1| 1.278e7| -3.539e6| 3.539e6| 0| 6.472e6
2| 1.278e7| -3.539e6| 3.539e6| 0| 6.472e6
3| 1.278e7| 3.535e6| -3.535e6| 0| 6.475e6
4| 1.332e7| -6.388e6| +6.388e6| 0| 4.285e6
5| 1.332e7| -6.388e6| +6.388e6| 0| 4.285e6
6| 1.332e7| 0.00000| 0.00000| 0| 7.692e6
7| 1.282e7| 1.324e7| 3.311e6| 0| 9.933e6
8| 1.282e7| -3.310e6| 3.310e6| 0| 6.622e6
9| 1.282e7| -7.284e6| 7.284e6| 0| 1.324e6
10| 1.289e7| 2.578e6| -2.578e6| 0| 6.978e6
11| 1.289e7| 2.578e6| -2.578e6| 0| 6.978e6

Case

1) Abaqus explicit
elastic isotropic behavior entered as a fully anisotropic stiffness matrix
orientation option used in 1,0,0 ; 0,1,0

2) Abaqus explicit
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
orientation option used in 1,0,0 ; 0,1,0

3) Abaqus explicit
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
NO orientation option was used

4) Abaqus implicit
elastic isotropic behavior entered as a fully anisotropic stiffness matrix
orientation option used in 1,0,0 0,1,0

5) Abaqus implicit
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
orientation option used in 1,0,0 0,1,0

6) Abaqus implicit
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
NO orientation option was used

7) Hand calculation of Cauchy using hyperelastic formulation and Hencky or natural strain
I worked with F, got C, got U got Elog, then Spk2 then Cauchy


8) Hand calculation of Spk2 using hyperelastic formulation
same formulation as above up to Spk2

9) “Spk2” or push backwards of Cauchy by using
“Spk2”=Rt Cauchy R


10) VUMAT abaqus manual example
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
orientation option used in 1,0,0 0,1,0


11) VUMAT abaqus manual
elastic isotropic behavior entered as two constants E=200E5; mu=0.3
NO orientation option was used

************************************************** **********

First thing which is good is that Von Mises is approximately the same for all the cases meaning that the theories behind them are all suitable for the this amount of strain and material behavior.

Second good thing is that sigma33 is 0 in all the cases which is also expected for simple shear.

The problem comes with the stress components.

************************************************** ***********
My questions are

1) Why case 1 and 4 give different results? is implicit and explicit analysis using different frames or stresses?

2) Why case 3 and 6 give different results? (same as 1)

3) Why cases 1-6 does not give 7 which is Cauchy stress? I remember reading in abaqus manula that the Cauchy stress is given as a result

4) In case of using UMAT which one of the stresses 1-9 (or another) should the subroutiner give back?

5) In case of using VUMAT which one of the stresses 1-9 (or another) should the subroutine give back? Note that the results you give back in this subroutine are shown rotated in the viewer (try something different than tensile and compression)

6) In case of using UMAT which one of the stresses 1-9 (or another) should you see in the viewer?

7) In case of using VUMAT which one of the stresses 1-9 (or another) should you see in the viewer?

8) Why cases 10 and 11 differ from the rest?

9) Why cases 10 and 11 give you the same result? should not case 11 give the result in the local coordinate like the difference seen in 1,2 wrt 3 and 4,5 wrt 6?

************************************************** **********
Does anybody has really clear in which system each stress tensor is given for the different abaqus results? Also does anyone has clear which stress measure is given in each case?

************************************************** ***********
Appendix

Here I paste an example code for one of the cases if anyone wants to try


*HEADING
1 ELEMENT TEST
**
*NODE
1,0,0,0
2,1,0,0
3,1,1,0
4,0,1,0
*ELEMENT, TYPE=CPE4R
1,1,2,3,4
*ELSET, ELSET=EBLOCK
1
*NSET, NSET=NBLOCK
1,2,3,4
*NSET, NSET=NTOP
3,4
*NSET, NSET=NBOTTOM
1,2
*NSET, NSET=NBOTTOMRIGHT
2
*NSET, NSET=NBOTTOMLEFT
1
*NSET, NSET=NTOPLEFT
4
*NSET, NSET=NLEFT
1,4
*NSET, NSET=NRIGTH
2,3
**
*SOLID SECTION,ELSET=EBLOCK,MATERIAL=MGRAIN1, ORIENTATION=OGRAINS
**
*ORIENTATION, NAME=OGRAINS, SYSTEM=RECTANGULAR, DEFINITION=COORDINATES
1,0,0, 0,1,0
**
*MATERIAL, NAME=MGRAIN1
**
*ELASTIC, TYPE=ANISOTROPIC
+2.69230769E+007,+1.15384615E+007,+2.69230769E+007 ,+1.15384615E+007,+1.15384615E+007,+2.69230769E+00 7,+0.00000000E+000,+0.00000000E+000,
+0.00000000E+000,+7.69230769E+006,+0.00000000E+000 ,+0.00000000E+000,+0.00000000E+000,+0.00000000E+00 0,+7.69230769E+006,+0.00000000E+000,
+0.00000000E+000,+0.00000000E+000,+0.00000000E+000 ,+0.00000000E+000,+7.69230769E+006,
**
**** Load Step 1 -------------------------------------------------------
*STEP,AMPLITUDE=RAMP, NLGEOM=YES
*STATIC
**
*BOUNDARY, TYPE=DISPLACEMENT
NTOP, 1,1,1
NTOP, 2,2,0.0
NBOTTOM, 1,2,0.0
*END STEP

Regards
Gregorio

[Jorgen]

Here are a few comments:

:arrow: It is not a good idea to use linear elasticity for large deformations, use hyperelasticity.

:arrow: You need to be careful with the boundary conditions in simple shear.

I tried the following input file:

** Uniaxial tension + one brick element + displacement control
*Parameter

# (input variables)
maxStrain = 0.01
mu = 1.0
kappa = 50.0

# (derived variables)
lam = exp(maxStrain)
uMax = lam - 1
tMax = 1.0
N = 500
dtMax = tMax / N

*Heading
*Node, nset=alln
1, 0, 0, 0
2, 0, 0, 1
3, 1, 0, 1
4, 1, 0, 0
5, 0, 1, 0
6, 0, 1, 1
7, 1, 1, 1
8, 1, 1, 0
*Nset, nset=nLeft
1, 2, 5, 6
*Nset, nset=nBot
1, 2, 3, 4
*Nset, nset=nBack
1, 4, 5, 8
*Nset, nset=nTop
5, 6, 7, 8
*Nset, nset=nRight
3, 4, 7, 8
*Element, type=C3D8R, elset=allE
1, 1, 2, 3, 4, 5, 6, 7, 8
*Solid section, elset=allE, material=mat
1.0
*Material, name=mat
*User material, constants=2
** NH vumat
** mu, kappa
1.0, 1.23
*Density
1000e-12
**Depvar
**0
*Boundary
nLeft, 1
nBot, 2
nBack, 3
*Amplitude, name=lamp, value=absolute
0.0, 0.0, <tMax>, <uMax>
*Step, nlgeom=yes
*Dynamic, explicit
, <tMax>
*Variable mass scaling, freq=1, dt=<dtMax>, elset=allE, type=below min
*Boundary, amplitude=lamp
nRight, 1, 1, <uMax>
*Output, field, op=NEW, number interval=100
*Node Output
u,
*Element Output
s, le, sdv
*Output, history, op=NEW, time interval=0.01
*Element Output, elset=alle
s, le
*End step


Then I ran the following octave-file:

F = [1 1 0; 0 1 0; 0 0 1]
B = F * F'
V = sqrtm(B)
E = logm(V)

## Neo-Hookean: T = mu * dev[Bstar]
Bstar = B - trace(B)/3 * eye(3)


The output from this octave file is as follows:

F =

1 1 0
0 1 0
0 0 1

B =

2 1 0
1 1 0
0 0 1

V =

1.34164 0.44721 0.00000
0.44721 0.89443 0.00000
0.00000 0.00000 1.00000

E =

0.21520 0.43041 0.00000
0.43041 -0.21520 0.00000
0.00000 0.00000 0.00000

Bstar =

0.66667 1.00000 0.00000
1.00000 -0.33333 0.00000
0.00000 0.00000 -0.33333


Note, since the shear modulus is equal to 1, the stress is the same as Bstar. When I then compared the output from ABAQUS (using the inputfile shown above) I get exactly the stress state indicated in Bstar.

The results seem fine to me.

[krikorik]

Dear Jorgen
Thank you for the answer.

:arrow: I agree that it is not a good idea to use elasticity for large deformations, nevertheless the von mises equivalent stress gave the almost same for all cases. The point is that if you don't use large deformations (I would say large rotations and infinitesimal deformations if you use a natural strain measure) all the stress measure components give you the same in either eulerian or lagrangian coordinates, implicit or explicit. I just want to dig in which differences arise from solving the same problem in different ways using linear elastic material behavior.

:arrow: The boundary conditions represent the case of a deck of cards being slided. I think this represents simple shear in difference of the pure shear case.

I run the code you posted and look at the nh model. You are using the left polar decomposition while I am using the right one. In your case that defines V in the current configuration while I define U in the initial one. This is the main difference in the strain computation.

I also used plain stress conditions and you are also using plane strain conditions but I specifically choose simple shear since in both cases S33 and E33 are 0 in linear elastic isotropic materials and no difference should come from that point.

I see that you get the stress from the deviatoric part of the left cauchy deformation tensor while I use the hydrostatic components as well.
That means we are not comparing apples with apples since your results come from a nh model.
What come next is my interpretation, not your words: I assume that you are saying that you can compute the stress for linear elastic isotropic case from
B = F * F'
V = sqrtm(B)
E = logm(V)
and use Sigma=C:E,
and feed it back to abaqus in the vumat subroutine

If C is isotropic as in this case, this stress tensor will have the same principal directions as V, B and the stress tensor you have obtained, meaning that these stresses are different in value but no other angular correction is necessary to pass them to abaqus.
If you don't think so please give me an explanation.


************************************************** **********
I wrote a code that pass to abaqus the stress computed as above for the same 2D system and got

3311721.43752585 1 1sigma
6621674.29865396 1 2sigma
0.000000000000000E+000 1 3sigma
6621674.29865396 2 1sigma
-3309960.64238762 2 2sigma
0.000000000000000E+000 2 3sigma
0.000000000000000E+000 3 1sigma
0.000000000000000E+000 3 2sigma
528.238568936009 3 3sigma

when passed to abaqus it will show in the viewer
S11= +7.285e6
S22= -7.285e6
S33= 0.0
S12= +1.324e6
Mises=1282e7

It still keeps the same mises but the components are different from all previous ones computed...

************************************************** *********
I was searching to see where those differences migth come from and found at least the solution for question 9 at the abaqus webpage. It is a bug


8O v63_4564
ABAQUS/Explicit writes incorrect values of the stress field
output in a global basis to the output database
file (.odb) for the solid elements with
hyperfoam/hyperelastic/vumat material, when an
orientation is associated with the section
definition (*SOLID SECTION, ORIENTATION). The
history output in the output database file is not
affected.

---------------------------------------------------------------
I was looking also at the values shown in my first post and I saw that the results with or without orientation are in an angular difference given by R which in this case is approx 26.6 degrees.

----------------------------------------------------------------
I still have not clear which reference system implicit and explicit are using.

---------------------------------------------------------------
One more weird thing:
try to look at the principal values of the stress in implicit and explicit for this shear case. You'll see that in the explicit case these values are rotated as expected but they are not in the implicit...

Regards
Gregorio