Abaqus - Strain energy of each ply

Hi, there,

I am working on analyzing a shell structure consisting of composite materials and trying to get the strain energy for each element's ply (composite layer). But it seems not working. The following is an example.

The 60 × 40 shell (plate) is subjected to a concentrated force (F=1e4) at the center and pined at four corners. The size of each finite element is 1 × 1, and its type is S4R. The shell consists of two plies, and the thickness of each ply is 0.5. Thus, the volume of each element is 1. Each pile is set as 3 integration points. Based on the Abaqus analyzed result, I can get the stress (including S11, S22, and S12), the strain (including E11, E22, and E12), and the integration point volume of each play in an element. However, I cannot use these data to match the elemental strain energy getting from the 'ESEDEN.' Taking element 1 as an example.

# Equation of shell elemental strain energy

u = (1/2) * (S11*E11 + S22*E22 + S12*E12)

# Elemental strain energy

u_a (elemental strain energy) = 331736.6875

# Stress

s_a1 (stress at integration point p1) = [153491.53125, 99121.953125, 0.0, -2863.34155273438]

s_a2 (stress at p2) = [74377.9609375, 50263.5703125, 0.0, -1431.5791015625]

s_a3 (stress at p3) = [-4735.6015625, 1405.18798828125, 0.0, 0.18339005112648]

s_a4 (stress at p4) = [1949.40014648438, -3915.52783203125, 0.0, -0.18339005112648]

s_a5 (stress at p5) = [-50737.09765625, -73807.6953125, 0.0, -1431.94580078125]

s_a6 (stress at p6) = [-103423.6015625, -143699.875, 0.0, -2863.70825195313]

# Strain

e_a1 (strain at integration point p1) = [7.58951663970947, 3.531658411026, 0.0, -1.7895884513855]

e_a2 (strain at p2)= [3.61816906929016, 1.88301622867584, 0.0, -0.894736886024475]

e_a3 (strain at p3) = [-0.353178441524506, 0.234374091029167, 0.0, 0.000114618786028586]

e_a4 (strain at p4) = [0.234374091029167, -0.353178441524506, 0.0, -0.000114618786028586]

e_a5 (strain at p5) = [-1.41426801681519, -4.32452583312988, 0.0, -0.894966125488281]

e_a6 (strain at p6) = [-3.06291007995605, -8.29587364196777, 0.0, -1.7898176908493]

# Integration point volume

v_a1 (volume at integration point p1) = 0.0833333358168602

v_a2 (volume at p2) = 0.333333343267441

v_a3 (volume at p3) = 0.0833333358168602

v_a4 (volume at p4) = 0.0833333358168602

v_a5 (volume at p5) = 0.333333343267441

v_a6 (volume at p6) = 0.0833333358168602

# Middle surface

When using data at p2 and p5 (the middle surface of each ply):

u2+u5 = 378630.34375 > u_a

# Volume influence

When considering the influence of integration point volumes (e.g., u1 = (1/2)* s_a1*e_al*v_a1):

u1+u2+u3+u4+u5+u6 = 252792.45528262202 < u_a