fix: correct the relation between the in-plane shear strain and the twisting curvature

structuralcodes/sections/section_integrators/_shell_integrator.py

@@ -47,7 +47,7 @@ def prepare_input(
         Returns:
             Tuple: (prepared_input, z_coords)
         """
-        strain = np.atleast_1d(strain)
+        strain = np.asarray(np.atleast_1d(strain), dtype='float64')
 
         layers = kwargs.get('layers')
         z_coords, dz = (None, None) if layers is None else layers
@@ -68,8 +68,11 @@ def prepare_input(
 
         # A positive in-plane strain gives tension
         # A positive curvature gives a negative strain for a positive z-value
+        # Note the factor 2 for in-plane shear strain due to twisting curvature
         for z in z_coords:
-            fiber_strain = strain[:3] - z * strain[3:]
+            fiber_strain = strain[:3].copy()
+            fiber_strain[:2] -= z * strain[3:5]
+            fiber_strain[2] -= 2 * z * strain[5]
             if integrate == 'stress':
                 integrand = material.constitutive_law.get_stress(fiber_strain)
             elif integrate == 'modulus':
@@ -81,23 +84,22 @@ def prepare_input(
             z_list.append(z)
 
         for r in geo.reinforcement:
-            z_r = r.z
-            material = r.material
-
-            fiber_strain = strain[:3] - z_r * strain[3:]
+            fiber_strain = strain[:3].copy()
+            fiber_strain[:2] -= r.z * strain[3:5]
+            fiber_strain[2] -= 2 * r.z * strain[5]
             eps_sj = r.T @ fiber_strain
 
             if integrate == 'stress':
-                sig_sj = material.constitutive_law.get_stress(eps_sj[0])
+                sig_sj = r.material.constitutive_law.get_stress(eps_sj[0])
                 integrand = r.local_stress_to_global_force * sig_sj
             elif integrate == 'modulus':
-                mod = material.constitutive_law.get_secant(eps_sj[0])
+                mod = r.material.constitutive_law.get_secant(eps_sj[0])
                 integrand = r.local_modulus_to_global_stiffness * mod
             else:
                 raise ValueError(f'Unknown integrate type: {integrate}')
 
             IA.append(integrand)
-            z_list.append(z_r)
+            z_list.append(r.z)
 
         MA = np.stack(IA, axis=0)
         prepared_input = [(z_list, MA)]
@@ -151,6 +153,9 @@ def integrate_modulus(
             B -= z_i * C_layer
             D += z_i**2 * C_layer
 
+        # This is necessary because the twisting moment is assumed proportional
+        # to the double twisting curvature
+        D[-1, -1] *= 2
         return np.block([[A, B], [B, D]])
 
     def integrate_strain_response_on_geometry(

tests/test_sections/test_shell_section.py

@@ -77,7 +77,7 @@ def test_integrate_strain_profile(
     Nxy = A_membrane_shear * eps_xy
     Mx = A_bending * (chi_x + nu * chi_y)
     My = A_bending * (chi_y + nu * chi_x)
-    Mxy = A_bending_shear * chi_xy
+    Mxy = A_bending_shear * 2 * chi_xy
 
     expected = np.array([Nx, Ny, Nxy, Mx, My, Mxy])
 
@@ -101,7 +101,7 @@ def test_integrate_strain_profile_stiffness_matrix(nu):
     A33 = E * t / (2 * (1 + nu))
     A44 = E * t**3 / 12 / (1 - nu**2)
     A14 = nu * A44
-    A55 = E * t**3 / 12 / (2 * (1 + nu))
+    A55 = E * t**3 / 12 / (2 * (1 + nu)) * 2
 
     A = np.array(
         [
@@ -159,10 +159,10 @@ def test_elastic_strain_profile(Ec, nu, t, nx, ny, nxy, mx, my, mxy):
     eps_x = (sig_x - nu * sig_y) / Ec
     eps_y = (sig_y - nu * sig_x) / Ec
     eps_xy = 2 * (1 + nu) * tau_xy / Ec
-    D = Ec * t**3 / 12
-    chi_x = (mx - nu * my) / D
-    chi_y = (my - nu * mx) / D
-    chi_xy = 2 * (1 + nu) * mxy / D
+    D = Ec * t**3 / (12 * (1 - nu**2))
+    chi_x = (mx - nu * my) / (D * (1 - nu**2))
+    chi_y = (my - nu * mx) / (D * (1 - nu**2))
+    chi_xy = mxy / (D * (1 - nu))
     expected = np.array([eps_x, eps_y, eps_xy, chi_x, chi_y, chi_xy])
 
     assert np.allclose(eps, expected)