From fb067c64e7bb7d7958016975bdb47ebd68c4b83d Mon Sep 17 00:00:00 2001 From: Morten Engen Date: Tue, 24 Mar 2026 21:29:32 +0100 Subject: [PATCH] Correct the relation between the in-plane shear strain and the twisting curvature --- .../section_integrators/_shell_integrator.py | 23 +++++++++++-------- tests/test_sections/test_shell_section.py | 12 +++++----- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/structuralcodes/sections/section_integrators/_shell_integrator.py b/structuralcodes/sections/section_integrators/_shell_integrator.py index b619bc30..3f5ad1fd 100644 --- a/structuralcodes/sections/section_integrators/_shell_integrator.py +++ b/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( diff --git a/tests/test_sections/test_shell_section.py b/tests/test_sections/test_shell_section.py index 89ecdfae..2b88dac5 100644 --- a/tests/test_sections/test_shell_section.py +++ b/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)