Skip to content

GJ calculation for composite sections using G_eff / E_eff #590

@giovanniboscu

Description

@giovanniboscu

hello,
In the retrieving "Section Properties" example [1], the torsional rigidity of a composite section is computed as:

ej = sec.get_ej()
sec.calculate_geometric_properties()
gj = sec.get_g_eff() / sec.get_e_eff() * ej

I am not sure this gives the correct G.J for composite sections.
My understanding is that get_ej() is already a torsion result weighted through the material elastic_modulus field. Then get_g_eff() / get_e_eff() applies only a global effective correction afterwards.
I tested this against another approach: rebuilding the same section, but assigning each material:

elastic_modulus_i = G_i

and then reading get_ej() as G.J.
For isotropic materials:

G_i = E_i / (2 * (1 + nu_i))

The two approaches give different results in the attached example.
So my question is:
Is the documented formula intended only as an approximate/effective correction?
And would using elastic_modulus_i = G_i in a torsion-only analysis be a valid way to obtain a directly G-weighted G.J?
Related previous discussion: #76
[1] https://sectionproperties.readthedocs.io/en/latest/examples/results/get_results.html

Material A  E = 3.0000e+04  G = 1.2500e+04  G/E = 0.416667
Material B  E = 2.0000e+05  G = 7.6923e+04  G/E = 0.384615

CASE A  G_eff/E_eff = 0.388796  E.J = 1.2959e+14  G.J = 5.0383e+13
CASE B  G_eff/E_eff = 0.388796  E.J = 1.6079e+14  G.J = 6.2513e+13

CASE A  G.J documented = 5.0383e+13  G.J carrier = 5.2415e+13  diff = -3.88%
CASE B  G.J documented = 6.2513e+13  G.J carrier = 6.3745e+13  diff = -1.93%

import matplotlib.pyplot as plt
from shapely import Polygon
from sectionproperties.analysis import Section
from sectionproperties.pre import CompoundGeometry, Geometry, Material

B, D         = 300.0, 400.0            # section width, height [mm]
STEEL_FRAC   = 0.50                    # Material B fraction of total area
E_C, NU_C    = 30.0e3,  0.20           # Material A: E [N/mm2], Poisson
E_S, NU_S    = 200.0e3, 0.30           # Material B: E [N/mm2], Poisson
G_C          = E_C / (2 * (1 + NU_C)) # Material A: G [N/mm2]
G_S          = E_S / (2 * (1 + NU_S)) # Material B: G [N/mm2]


def make_rect(x1, x2, y1, y2, mat):
    return Geometry(geom=Polygon([(x1,y1),(x2,y1),(x2,y2),(x1,y2)]), material=mat)

def make_geometry(case_id, g_carrier=False):
    ea, eb = (G_C, G_S) if g_carrier else (E_C, E_S)  # E or G as carrier
    mat_a  = Material("Material A", ea, NU_C, 30.0,  2.4e-6,  "lightgrey")
    mat_b  = Material("Material B", eb, NU_S, 500.0, 7.85e-6, "grey")
    x1, x2 = -B/2, B/2
    y1, y2 = -D/2, D/2
    h_s    = D * STEEL_FRAC
    if case_id == "A":
        h     = h_s / 2
        parts = [make_rect(x1,x2, y1,     y1+h,   mat_b),
                 make_rect(x1,x2, y1+h,   y2-h,   mat_a),
                 make_rect(x1,x2, y2-h,   y2,     mat_b)]
    elif case_id == "B":
        parts = [make_rect(x1,x2, y1,     y1+h_s, mat_b),
                 make_rect(x1,x2, y1+h_s, y2,     mat_a)]
    geom = CompoundGeometry(geoms=parts)
    geom.create_mesh(mesh_sizes=1000.0)
    return geom

def analyse(case_id):
    sec = Section(make_geometry(case_id))
    sec.calculate_frame_properties()
    ej = sec.get_ej()
    sec.calculate_geometric_properties()
    ratio = sec.get_g_eff() / sec.get_e_eff()
    return sec, ratio, ej, ratio * ej

def analyse_g_carrier(case_id):
    sec = Section(make_geometry(case_id, g_carrier=True))
    sec.calculate_frame_properties()
    return sec.get_ej()                # E.J read as G.J (E_SP := G_i)

fig, axes = plt.subplots(1, 2, figsize=(10, 5))
results   = {}
for ax, case_id in zip(axes, ("A", "B")):
    sec, ratio, ej, gj = analyse(case_id)
    results[case_id]   = (ratio, ej, gj)
    sec.plot_mesh(materials=True, ax=ax)
    ax.set_title(f"Case {case_id}")
plt.tight_layout()
plt.show()

print(f"Material A  E = {E_C:.4e}  G = {G_C:.4e}  G/E = {G_C/E_C:.6f}")
print(f"Material B  E = {E_S:.4e}  G = {G_S:.4e}  G/E = {G_S/E_S:.6f}")
print()
for case_id in ("A", "B"):
    ratio, ej, gj = results[case_id]
    print(f"CASE {case_id}  G_eff/E_eff = {ratio:.6f}  E.J = {ej:.4e}  G.J = {gj:.4e}")

# --- g-carrier comparison ---
print()
for case_id in ("A", "B"):
    gj_doc     = results[case_id][2]
    gj_carrier = analyse_g_carrier(case_id)
    diff       = (gj_doc / gj_carrier - 1.0) * 100.0
    print(f"CASE {case_id}  G.J documented = {gj_doc:.4e}  G.J carrier = {gj_carrier:.4e}  diff = {diff:+.2f}%")

Metadata

Metadata

Assignees

No one assigned

    Labels

    documentationImprovements or additions to documentation

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions