Topic: MITC4 stresses

Consider the membrane patch test by Macneal, Richard H., and Robert L. Harder. "A proposed standard set of problems to test finite element accuracy." Finite elements in Analysis and Design 1.1 (1985): 3-20 [link] for the MITC4 shell.

INPUT FILE:

patch_0.out
"Patch test of MITC4Shell elements, shear"
StaticStructural nsteps 1
domain 3dshell
OutputManager tstep_all dofman_all element_all
ndofman 8 nelem 5 ncrosssect 1 nmat 1 nbc 4 nic 0 nltf 1 nset 5
node 1 coords 3 0.04 0.02 0.00
node 2 coords 3 0.18 0.03 0.00
node 3 coords 3 0.16 0.08 0.00
node 4 coords 3 0.08 0.08 0.00
node 5 coords 3 0.00 0.00 0.00
node 6 coords 3 0.24 0.00 0.00
node 7 coords 3 0.24 0.12 0.00
node 8 coords 3 0.00 0.12 0.00
mitc4shell 1 nodes 4 5 6 2 1
mitc4shell 2 nodes 4 6 7 3 2
mitc4shell 3 nodes 4 7 8 4 3
mitc4shell 4 nodes 4 8 5 1 4
mitc4shell 5 nodes 4 1 2 3 4
SimpleCS 1 thick 0.001 material 1 set 1 drillType 1 relDrillStiffness 0.01
IsoLE 1 d 2.  E 1e6  n 0.25  tAlpha 0.000012
BoundaryCondition 1 loadTimeFunction 1 dofs 3 1 2 3 values 3 0       0       0       set 2
BoundaryCondition 2 loadTimeFunction 1 dofs 2 1 2   values 2 0.00024 0.00012         set 3
BoundaryCondition 3 loadTimeFunction 1 dofs 2 1 2   values 2 0.00030 0.00024         set 4
BoundaryCondition 4 loadTimeFunction 1 dofs 2 1 2   values 2 0.00006 0.00012         set 5
ConstantFunction 1 f(t) 1.0
Set 1 elementranges {(1 5)}
Set 2 nodes 1 5
Set 3 nodes 1 6
Set 4 nodes 1 7
Set 5 nodes 1 8

OUTPUT (selected):

element 1 (       1):
  GP 1 :  forces      1.3333e+00 1.3333e+00 0.0000e+00 0.0000e+00 0.0000e+00 4.0000e-01
          moments     0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
          strains     1.0000e-03 1.0000e-03 0.0000e+00 0.0000e+00 0.0000e+00 1.0000e-03
          curvatures  0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
          GP 1.1 :    strains     1.0000e-03 1.0000e-03 0.0000e+00 0.0000e+00 0.0000e+00 1.0000e-03
                      stresses    1.3333e+03 1.3333e+03 0.0000e+00 0.0000e+00 4.0000e+02 0.0000e+00
          GP 1.2 :    strains     1.0000e-03 1.0000e-03 0.0000e+00 0.0000e+00 0.0000e+00 1.0000e-03
                      stresses    1.3333e+03 1.3333e+03 0.0000e+00 0.0000e+00 4.0000e+02 0.0000e+00

ISSUE:
Looks like that σxy is messed with σzx.

Re: MITC4 stresses

Hi,
thanks for the notice. We've fixed the wrong ordering of stress components in printing the output, current version should be working correctly.

Edita

Re: MITC4 stresses

Hi Edita,
thanks for the fix.
However I still have some strange findings from this element regarding the stress recovery. Please consider the following case:

3 o ---------- o 4
  |            |
  |            |
  |            | Ly = 1
  |            |
1 o ---------- o 2
    Lx = 1

Nodes 1, 3 are assumed fixed. Vertical load is applied on 2, 4.
The corresponding input file reads:

bending_plate.out
"bending plate"
StaticStructural nsteps 1
domain 3dshell
OutputManager tstep_all dofman_all element_all
ndofman 4 nelem 1 ncrosssect 1 nmat 1 nbc 4 nic 0 nltf 1 nset 5
node 1 coords 3       0.0000       0.0000       0.0000
node 2 coords 3       1.0000       0.0000       0.0000
node 3 coords 3       0.0000       1.0000       0.0000
node 4 coords 3       1.0000       1.0000       0.0000
mitc4shell    1 nodes 4    1    2    4    3
SimpleCS 1 thick 1.000000 material 1 set 1 drillType 1 relDrillStiffness 0.0
IsoLE 1 d 2.  E 1000000.000000  n 0.000000  tAlpha 0.000012
BoundaryCondition 1 loadTimeFunction 1 dofs 6 1 2 3 4 5 6 values 6   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000 set 2
BoundaryCondition 2 loadTimeFunction 1 dofs 6 1 2 3 4 5 6 values 6   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000   0.00000000 set 3
NodalLoad 3 loadTimeFunction 1 dofs 1 3 components 1  -0.50000000 set 4
NodalLoad 4 loadTimeFunction 1 dofs 1 3 components 1  -0.50000000 set 5
ConstantFunction 1 f(t) 1.0
Set 1 elementranges {(1 1)}
Set 2 nodes 1 1
Set 3 nodes 1 3
Set 4 nodes 1 2
Set 5 nodes 1 4

And the corresponding output is:

Element output:
---------------
element 1 (       1):
  GP 1 :  forces      0.0000e+00 0.0000e+00 0.0000e+00 -4.1752e-17 -1.0000e+00 0.0000e+00
          moments     5.0000e-01 2.7835e-17 0.0000e+00 0.0000e+00 0.0000e+00 3.0565e-17
          strains     0.0000e+00 0.0000e+00 0.0000e+00 -8.3504e-23 -2.0000e-06 0.0000e+00
          curvatures  6.0000e-06 3.3402e-22 0.0000e+00 0.0000e+00 0.0000e+00 -3.0126e-22
          GP 1.1 :    strains     -1.7321e-06 -9.6422e-23 0.0000e+00 -8.3504e-23 -2.0000e-06 -2.1176e-22
                      stresses    -1.7321e+00 -9.6422e-17 0.0000e+00 -4.1752e-17 -1.0000e+00 -1.0588e-16
          GP 1.2 :    strains     1.7321e-06 9.6422e-23 0.0000e+00 -8.3504e-23 -2.0000e-06 2.1176e-22
                      stresses    1.7321e+00 9.6422e-17 0.0000e+00 -4.1752e-17 -1.0000e+00 1.0588e-16
  GP 2 :  forces      0.0000e+00 0.0000e+00 0.0000e+00 -4.1752e-17 -1.0000e+00 0.0000e+00
          moments     5.0000e-01 2.7835e-17 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
          strains     0.0000e+00 0.0000e+00 0.0000e+00 -8.3504e-23 -2.0000e-06 0.0000e+00
          curvatures  6.0000e-06 3.3402e-22 0.0000e+00 0.0000e+00 0.0000e+00 -5.4578e-22
          GP 2.1 :    strains     -1.7321e-06 -9.6422e-23 0.0000e+00 -8.3504e-23 -2.0000e-06 0.0000e+00
                      stresses    -1.7321e+00 -9.6422e-17 0.0000e+00 -4.1752e-17 -1.0000e+00 0.0000e+00
          GP 2.2 :    strains     1.7321e-06 9.6422e-23 0.0000e+00 -8.3504e-23 -2.0000e-06 0.0000e+00
                      stresses    1.7321e+00 9.6422e-17 0.0000e+00 -4.1752e-17 -1.0000e+00 0.0000e+00
  GP 3 :  forces      0.0000e+00 0.0000e+00 0.0000e+00 -1.1187e-17 -1.0000e+00 0.0000e+00
          moments     5.0000e-01 7.4583e-18 0.0000e+00 0.0000e+00 0.0000e+00 -7.6412e-18
          strains     0.0000e+00 0.0000e+00 0.0000e+00 -2.2375e-23 -2.0000e-06 0.0000e+00
          curvatures  6.0000e-06 8.9500e-23 0.0000e+00 0.0000e+00 0.0000e+00 -8.9500e-23
          GP 3.1 :    strains     -1.7321e-06 -2.5836e-23 0.0000e+00 -2.2375e-23 -2.0000e-06 5.2940e-23
                      stresses    -1.7321e+00 -2.5836e-17 0.0000e+00 -1.1187e-17 -1.0000e+00 2.6470e-17
          GP 3.2 :    strains     1.7321e-06 2.5836e-23 0.0000e+00 -2.2375e-23 -2.0000e-06 -5.2940e-23
                      stresses    1.7321e+00 2.5836e-17 0.0000e+00 -1.1187e-17 -1.0000e+00 -2.6470e-17
  GP 4 :  forces      0.0000e+00 0.0000e+00 0.0000e+00 -1.1187e-17 -1.0000e+00 0.0000e+00
          moments     5.0000e-01 7.4583e-18 0.0000e+00 0.0000e+00 0.0000e+00 -1.5282e-17
          strains     0.0000e+00 0.0000e+00 0.0000e+00 -2.2375e-23 -2.0000e-06 0.0000e+00
          curvatures  6.0000e-06 8.9500e-23 0.0000e+00 0.0000e+00 0.0000e+00 -3.3402e-22
          GP 4.1 :    strains     -1.7321e-06 -2.5836e-23 0.0000e+00 -2.2375e-23 -2.0000e-06 1.0588e-22
                      stresses    -1.7321e+00 -2.5836e-17 0.0000e+00 -1.1187e-17 -1.0000e+00 5.2940e-17
          GP 4.2 :    strains     1.7321e-06 2.5836e-23 0.0000e+00 -2.2375e-23 -2.0000e-06 -1.0588e-22
                      stresses    1.7321e+00 2.5836e-17 0.0000e+00 -1.1187e-17 -1.0000e+00 -5.2940e-17

Assuming tension is positive, X.1 Gausspoint is in the bottom and X.2 in the top, it looks to me that the strains/stresses have opposite signs.

Re: MITC4 stresses

I looked a bit closer at this, but it looks right to me.
You push down (-z), which means we will have tension (positive normal force) on the top layer (+z) and compression (negative normal force) on the bottom layer (-z), which matches up with what you get here?

                         | F = 0.5N
    _____________________v
    |       <---> X.2    |
1/2 o--------------------o 2/4
    |_______>---<_X.1____|

Re: MITC4 stresses

I guess you are right, I have uploaded a mismatching input file.
Please try that again, but reverse the element connectivity as:

mitc4shell    1 nodes 4    1    3    4    2

The above results in:

  GP 1 :  forces      0.0000e+00 0.0000e+00 0.0000e+00 4.1752e-17 -1.0000e+00 0.0000e+00
          moments     -5.0000e-01 6.2628e-17 0.0000e+00 0.0000e+00 0.0000e+00 1.5096e-16
          strains     0.0000e+00 0.0000e+00 0.0000e+00 8.3504e-23 -2.0000e-06 0.0000e+00
          curvatures  -6.0000e-06 7.5154e-22 0.0000e+00 0.0000e+00 0.0000e+00 3.0809e-21
          GP 1.1 :    strains     1.7321e-06 -2.1695e-22 0.0000e+00 8.3504e-23 -2.0000e-06 -1.0459e-21
                      stresses    1.7321e+00 -2.1695e-16 0.0000e+00 4.1752e-17 -1.0000e+00 -5.2294e-16
          GP 1.2 :    strains     -1.7321e-06 2.1695e-22 0.0000e+00 8.3504e-23 -2.0000e-06 1.0459e-21
                      stresses    -1.7321e+00 2.1695e-16 0.0000e+00 4.1752e-17 -1.0000e+00 5.2294e-16

I.e. the stresses change signs, which is not supposed to happen as they are reported in the global coordinate system, right?
Therefore, this does not look as a sign problem as I have initially suggested, but probably a stress transformation problem.

Re: MITC4 stresses

Hi,
I think the problem here is in the GP numbering rather than in transformations.  When you reverse the connectivity you basically switch the local z-coordinate upside down and this influences the GP numbering. So with connectivity 1-3-4-2 GP1.1 is located in what seems to us as a top layer (positive global z, however negative local z) and GP1.2 in bottom layer.

Re: MITC4 stresses

It's definitely the case that by turning the element will make switch the order of GPs, but there is still something quite wrong here.
We should get the moments in global c.s either way, so moments should stay identical.

Using the simplest cases (this is exactly the local c.s. used

mitc4shell    1 nodes 4    4    3    1    2
# debug printouts when computing moment:
z = -2.886751e-01
localStress (6): 
-1.732e+00   2.652e-16   0.000e+00   4.235e-16  -1.000e+00   0.000e+00  
z = 2.886751e-01
localStress (6): 
 1.732e+00  -2.652e-16   0.000e+00   4.235e-16  -1.000e+00   0.000e+00  
mLocal (6): 
 5.000e-01  -7.655e-17   0.000e+00   0.000e+00   0.000e+00   6.113e-17  
GtoL (3 x 3): 
 1.000e+00   0.000e+00   0.000e+00  
 0.000e+00   1.000e+00   0.000e+00  
 0.000e+00   0.000e+00   1.000e+00  
answer (6): 
 5.000e-01  -7.655e-17   0.000e+00   0.000e+00   0.000e+00   6.113e-17  

should give exactly the same moments as the one flipped around

mitc4shell    1 nodes 4    2    1    3    4
# debug printouts when computing moments:
z = -2.886751e-01  (local)
localStress (6): 
 1.732e+00  -3.375e-16   0.000e+00   4.235e-16   1.000e+00   1.059e-16  
z = 2.886751e-01  (local)
localStress (6): 
-1.732e+00   3.375e-16   0.000e+00   4.235e-16   1.000e+00  -1.059e-16  
mLocal (6): 
-5.000e-01   9.742e-17   0.000e+00   0.000e+00   0.000e+00  -3.056e-17  
GtoL (3 x 3): 
 1.000e+00   0.000e+00   0.000e+00  
 0.000e+00  -1.000e+00   0.000e+00  
 0.000e+00   0.000e+00  -1.000e+00  
answer (6): 
-5.000e-01   9.742e-17   0.000e+00   0.000e+00   0.000e+00   3.056e-17  

But! The flaw is not with the transformation, it's that mLocal is completely incorrect.
Looking at the first component of M_xx (which happens to be 0.5Nm in the output above), the code currently does:

M_xx = integral sigma_xx *z dz

But this is wrong, it must be

M_xy = integral sigma_xx * z dz

In particular, it's this line of code

mLocal.add(w, localStress);

is incorrect for all the components. Many of the things computed here makes no sense.


So, the moment vector must follow the exact same form as the stress vector, else transformStressVectorTo will not work correctly as all.
M = [M_xx, M_yy, M_zz, M_yz, M_xz, M_xy]
and similarly, the curvature tensor in voigt form should probably follow the strain Voigt form (and must therefor also use transformStrainVectorTo)
K =  [K_xx, K_yy, K_zz, 2*K_yz, 2*K_xz, 2*K_xy]

So, the moment vector, in this case, should be

moments = 0.0   0.0   0.0   0.0   0.0   0.5

So, quite a bit of code needs to change in giveMidplaneIPValue

Re: MITC4 stresses

Hmm, I looked at implementing a fix here, and, well, not I'm unsure about everything I've said so far.
The current code is definitely wrong somewhere though, as it doesn't give the same results when rotating the element, despite computing global values.
At the very least, the StructuralMaterial::transformStressVectorTo method isn't adapted for the components of the shell moment tensor.


I compare with a few other plates, and, most of them have something  like

        answer.at(1) = help.at(1); // mx
        answer.at(2) = help.at(2); // my
        answer.at(3) = 0.0;        // mz
        answer.at(4) = 0.0;        // myz
        answer.at(5) = 0.0;        // mxz
        answer.at(6) = help.at(3); // mxy

And, well, most of them don't  ever bother computing global values, so, if we never transform this tensor, it's just a list of values, printed as they are.
So, you could put them in any order, and it would work fine.


If we wanted to properly store the shell moment as a tensor, wouldn't be unsymmetric with M_xy and M_yx be the primary bending modes of a plate? But, i'm honestly unsure about everything here now other than that global values should remain constant in these examples.