Wrong projection of stress values from Gauss point to the nodes

[mmoti]

Dear Feap community,

I have a problem to plot correctly my stress values at nodes.
Here is what I have already done:

I have implemented my own element for a 3D shell problem (Kirchhoff-Love with quadratic 8 and 9 nodes element) with only 3 displacement DOFs.
In a case where I apply a symmetric boundary condition I expect a symmetric result in displacement as well as  in stresses contour plot.
The applied boundary condition is the simplest test of in-homogeneous on a flat membrane shape. (homogeneous tests shows correct plots)
The displacement is perfectly symmetric as undoubtedly expected!! (plot,count,1  2  3 )
However, the stress plot is not. (plot,stre,1  2  3 = sigma_11 or _22 or _12 )  I have checked the stress values at the symmetric Gauss points in  symmetric positioned elements. The stress values at Gauss points are equal and for shear stresses (sigma_12) have different signs (- or +) which is also expected.
I have checked the same boundary value problem with the shell3d of Feap and there I got the symmetric correct plots for stresses.
I have checked the stress projection subroutine of Feap and I use it in a similar way.
the "pstyp" as it was suggested in the feap programmer manual, I chose 2, however I also tried with pstyp 3.

Long story short, when I am using quadratic shape functions (8 or 9 nodes quadrilateral elements) for 3D shell problem with 3 DOF (only displacement) do I need any special adjustment in the stress nodal projection subroutine?

Thank you.
I look forward to hearing your helpful advices.
M.Moti

[Prof. R.L. Taylor]

How are you projecting the stress values?  What does stress, node, 1, numnp give?   In some states feap expects the stress order to be sig_11, sig_22, sig_33, sig_12 for 2-d states (to use principal stress routines built in).  If you post your routines we can tell more.

[mmoti]

Thanks for your reply.

I attached my subroutine  for nodal stress projecting  as well as the one in the shell3d elmt of Feap. ( should I attach my complete  element subroutine? )
I also put some pictures to show the boundary problem and stress result.

*******************************************************************
here is the result of stress, node, 1, numnp:

 feap * * LINEar-ELAStic.....3D.....QUADrilateral-ELEMent                     

   N o d a l   P r o j e c t i o n s

   Node  1-Pr.Value  2-Pr.Value  3-Pr.Value  1-Pr.Angle
               I_1 Value   J_2 Value   J_3 Value
              1 Value     2 Value     3 Value

       1     2.3129E+02  5.5310E+01  3.6246E+01  0.0000E+00
              1.0762E+02  1.8625E+02  3.2012E+05
              2.3129E+02  5.5310E+01  3.6246E+01

as well as node 17:

   N o d a l   P r o j e c t i o n s

   Node  1-Pr.Value  2-Pr.Value  3-Pr.Value  1-Pr.Angle
               I_1 Value   J_2 Value   J_3 Value
               1 Value     2 Value     3 Value

      17    2.0024E+02  4.6883E+01 -3.4369E+00  0.0000E+00
              8.1230E+01  1.8376E+02  2.3997E+05
              2.0024E+02  4.6883E+01 -3.4369E+00
*******************************************************************
Due to plane stress model,  sig_33=0 and my order in sigma plotting is:
xx-stress   yy-stress   xy-stress   xx-moment    yy-moment    xy-moment ...

Thank you for your attention.

[mmoti]

I thought, perhaps this addition data helps:

a comparison of the stress values at the Gauss Points (9 GP) for an 8 noded element at two symmetric elements, element 1 and 3.
Order of numbering of the Gauss Points is also mentioned.
GP 1  in element 1 is considered  at the symmetric position to the GP 4 of element 3.

at each GP the value of the "  dvol = [det(J) * weight factor]   "  is showed.   

[Prof. R.L. Taylor]

You are only projecting one quadrature point.  You need to see the routines in feap which loop over the number of quadrature points when doing the projection (usually a loop: do l = 1,lint).

[mmoti]

Thank you very much.
Done