import numpy as np # Material: m1={'E1':133860, 'E2':7706, 'v12':0.301, 'G12':4306} laminateThickness = 10 thi = laminateThickness n = 10 layup1 = [ {'mat':m1 , 'ori': 0 , 'thi':thi/2} , {'mat':m1 , 'ori': 90 , 'thi':thi/2}] layup2 = [ {'mat':m1 , 'ori': 0 , 'thi':thi/n} , {'mat':m1 , 'ori': 90 , 'thi':thi/n}, {'mat':m1 , 'ori': 0 , 'thi':thi/n} , {'mat':m1 , 'ori': 90 , 'thi':thi/n}, {'mat':m1 , 'ori': 0 , 'thi':thi/n} , {'mat':m1 , 'ori': 90 , 'thi':thi/n}, {'mat':m1 , 'ori': 0 , 'thi':thi/n} , {'mat':m1 , 'ori': 90 , 'thi':thi/n}, {'mat':m1 , 'ori': 0 , 'thi':thi/n} , {'mat':m1 , 'ori': 90 , 'thi':thi/n}] # The code below is taken from: # https://folk.ntnu.no/nilspv/TMM4175/computational-procedures.html def S2D(m): return np.array([[ 1/m['E1'], -m['v12']/m['E1'], 0], [-m['v12']/m['E1'], 1/m['E2'], 0], [ 0, 0, 1/m['G12']]], float) def T2Ds(a): c,s = np.cos(np.radians(a)), np.sin(np.radians(a)) return np.array([[ c*c , s*s , 2*c*s], [ s*s , c*c , -2*c*s], [-c*s, c*s , c*c-s*s]], float) def T2De(a): c,s = np.cos(np.radians(a)), np.sin(np.radians(a)) return np.array([[ c*c, s*s, c*s ], [ s*s, c*c, -c*s ], [-2*c*s, 2*c*s, c*c-s*s ]], float) def Q2D(m): S=S2D(m) return np.linalg.inv(S) def Q2Dtransform(Q,a): return np.dot(np.linalg.inv( T2Ds(a) ) , np.dot(Q,T2De(a)) ) def computeB(layup): B=np.zeros((3,3),float) hbot = -laminateThickness/2 # bottom of first layer for layer in layup: m = layer['mat'] Q = Q2D(m) a = layer['ori'] Qt = Q2Dtransform(Q,a) htop = hbot + layer['thi'] # top of current layer B = B + (1/2)*Qt*(htop**2-hbot**2) hbot=htop # for the next layer return B B1 = computeB(layup1) B2 = computeB(layup2) print(B1[0,0]) print(B2[0,0])