In [33]:

    

import numpy as np
from sympy.abc import z

from meshless.composite.laminate import read_stack

E1 = 49.627e9
E2 = 15.43e9
nu12 = 0.38
G12 = 4.8e9
G13 = 4.8e9
G23 = 4.8e9
laminaprop = (E1, E2, nu12, G12, G13, G23)

tmap = {
      45: 0.143e-3,
     -45: 0.143e-3,
       0: 1.714e-3
    }
X = 1
angles = [-45, +45, 0, +45, -45, 0]*X + [0, -45, +45, 0, +45, -45]*X
plyts = [tmap[angle] for angle in angles]

lam = read_stack(angles, plyts=plyts, laminaprop=laminaprop)

def calc_scf(self):
    """Calculate shear correction factors

    Reference:

        Vlachoutsis, S. "Shear correction factors for plates and shells",
        Int. Journal for Numerical Methods in Engineering, Vol. 33,
        1537-1552, 1992.
        
        http://onlinelibrary.wiley.com/doi/10.1002/nme.1620330712/full
        

    Using "one shear correction factor" (see reference), assuming:

    - constant G13, G23, E1, E2, nu12, nu21 within each ply
    - g1 calculated using z at the middle of each ply
    - zn1 = Laminate.offset

    Returns
    -------

    k13, k23 : tuple
        Shear correction factors. Also updates attributes: `scf_k13`
        and `scf_k23`.

    """
    D1 = 0
    R1 = 0
    den1 = 0
    
    D2 = 0
    R2 = 0
    den2 = 0
    
    offset = self.offset
    zbot = -self.t/2 + offset
    z1 = zbot
    
    for ply in self.plies:
        z2 = z1 + ply.t
        e1 = (ply.matobj.e1 * np.cos(np.deg2rad(ply.theta)) +
              ply.matobj.e2 * np.sin(np.deg2rad(ply.theta)))
        e2 = (ply.matobj.e2 * np.cos(np.deg2rad(ply.theta)) +
              ply.matobj.e1 * np.sin(np.deg2rad(ply.theta)))
        nu12 = (ply.matobj.nu12 * np.cos(np.deg2rad(ply.theta)) +
              ply.matobj.nu21 * np.sin(np.deg2rad(ply.theta)))
        nu21 = (ply.matobj.nu21 * np.cos(np.deg2rad(ply.theta)) +
          ply.matobj.nu12 * np.sin(np.deg2rad(ply.theta)))
        
        D1 += e1 / (1 - nu12*nu21)
        R1 += D1*((z2 - offset)**3/3 - (z1 - offset)**3/3)
        g13 = ply.matobj.g13
        d1 = g13 * ply.t
        den1 += d1 * (self.t / ply.t) * D1**2*(15*offset*z1**4 + 30*offset*z1**2*zbot*(2*offset - zbot) - 15*offset*z2**4 + 30*offset*z2**2*zbot*(-2*offset + zbot) - 3*z1**5 + 10*z1**3*(-2*offset**2 - 2*offset*zbot + zbot**2) - 15*z1*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2) + 3*z2**5 + 10*z2**3*(2*offset**2 + 2*offset*zbot - zbot**2) + 15*z2*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2))/(60*g13)

        D2 += e2 / (1 - nu12*nu21)
        R2 += D2*((z2 - self.offset)**3/3 - (z1 - self.offset)**3/3)
        g23 = ply.matobj.g23
        d2 = g23 * ply.t
        den2 += d2 * (self.t / ply.t) * D2**2*(15*offset*z1**4 + 30*offset*z1**2*zbot*(2*offset - zbot) - 15*offset*z2**4 + 30*offset*z2**2*zbot*(-2*offset + zbot) - 3*z1**5 + 10*z1**3*(-2*offset**2 - 2*offset*zbot + zbot**2) - 15*z1*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2) + 3*z2**5 + 10*z2**3*(2*offset**2 + 2*offset*zbot - zbot**2) + 15*z2*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2))/(60*g23)
        
        z1 = z2

    self.scf_k13 = R1**2 / den1
    self.scf_k23 = R2**2 / den2

    return self.scf_k13, self.scf_k23

print(calc_scf(lam))


    

    




(0.76833028877157949, 0.77602930040164342)

In [32]:

    

from sympy.abc import z
import sympy

sympy.var('z1, z2, D1, offset, zbot, g13')
g1 = - D1*sympy.integrate(z - offset, (z, zbot, z))
expr = sympy.integrate((g1**2 / g13), (z, z1, z2))
print(expr.simplify())


    

    




D1**2*(15*offset*z1**4 + 30*offset*z1**2*zbot*(2*offset - zbot) - 15*offset*z2**4 + 30*offset*z2**2*zbot*(-2*offset + zbot) - 3*z1**5 + 10*z1**3*(-2*offset**2 - 2*offset*zbot + zbot**2) - 15*z1*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2) + 3*z2**5 + 10*z2**3*(2*offset**2 + 2*offset*zbot - zbot**2) + 15*z2*zbot**2*(4*offset**2 - 4*offset*zbot + zbot**2))/(60*g13)