对于任何使用复杂方程组的项目(如复合力学分析或推力/流体流动分析),这个问题已经出现了几个月。不过,简单的计算、迭代方法和Numpy矩阵本身似乎可以很好地工作。当然,我已经尝试过完全卸载和重新安装Python和Pycharm,但同样的问题依然存在。
在过去的几个月里,任何包含更复杂数学的代码(如上述代码)都将输出不同的值,尽管所有输入值和计算都是常量。不同的代码没有变化或随机性。我注意到这些不正确的输出通常发生在数组/矩阵中,我可以看出输出值不正确,因为预期的数字反而大得离谱。
例如,矩阵中的一个元素预期为5.197e+7,但相反,代码为同一元素输出3.322e+257。然而,在第二次运行代码时,代码会为完全相同的元素生成2.822e+204的输出。只有在第三次运行时,代码才会输出该元素的正确值,即5.197e+7。为了输出正确的值,这个过程可以对相同的、不变的代码进行1到15次单独的运行。
另一个有趣的方面是,我编写的计算通常需要对所述代码进行多次迭代。但是,即使我正在重置临时保存所有值(不再影响计算的最终值除外)的数组,任何导致这些“错误”的因素都会在代码中进行,直到迭代结束。据我所知,情况不应该是这样,因为在迭代结束时,代码会将所有值(已知初始值除外)设置为0并重新计算。这意味着代码会不断地产生相同的错误。
以下是该项目预期产出与实际产出的一些例子。
第一次迭代的预期输出
NM = [[ 2.57939977e+04]
[ 3.03926820e+04]
[-3.55271368e-13]
[ 1.00000000e+00]
[ 1.00000000e+00]
[ 2.60208521e-16]]
[A] =
[[1.50155575e+08 4.45004838e+07 0.00000000e+00]
[4.44974703e+07 1.07288531e+08 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 5.19662175e+07]]
[D] =
[[60.771 7.663 4.019]
[ 7.659 12.322 4.019]
[ 4.019 4.019 9.567]]
第一次尝试:第一次迭代的实际输出
NM =
[[2.57939977e+004]
[7.26576487e+225]
[2.35904846e+253]
[7.25895469e+242]
[1.18107381e+291]
[1.61312569e+291]]
[A] =
[[1.50155575e+08 4.45004838e+07 0.00000000e+00]
[4.44974703e+07 1.07288531e+08 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 5.19662175e+07]]
[D] =
[[60.771 7.663 4.019]
[ 7.659 12.322 4.019]
[ 4.019 4.019 9.567]]
第二次尝试:第一次迭代
NM = [[2.18479897e+158]
[3.03926820e+004]
[1.62552246e+034]
[1.00000121e+000]
[1.07935186e+000]
[2.60208521e-016]]
[A] =
[[1.50155575e+08 4.45004838e+07 0.00000000e+00]
[4.44974703e+07 1.07288531e+08 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 5.19662175e+07]]
[D] =
[[60.771 7.663 4.019]
[ 7.659 12.322 4.019]
[ 4.019 4.019 9.567]]
第三次尝试:第一次迭代
NM = [[ 2.57939977e+004]
[ 3.03926820e+004]
[-3.55271368e-013]
[ 4.54183603e+225]
[ 6.28847407e+094]
[ 2.60208521e-016]]
[A] =
[[1.50155575e+08 4.45004838e+07 0.00000000e+00]
[4.44974703e+07 1.07288531e+08 0.00000000e+00]
[0.00000000e+00 0.00000000e+00 5.19662175e+07]]
[D] =
[[60.771 7.663 4.019]
[ 7.659 12.322 4.019]
[ 4.019 4.019 9.567]]
第四次尝试:第一次迭代
NM = [[ 2.57939977e+04]
[ 3.03926820e+04]
[-3.55271368e-13]
[ 1.00000000e+00]
[ 1.00000000e+00]
[ 2.60208521e-16]]
[A] =
[[1.50155575e+008 4.45004838e+007 0.00000000e+000]
[4.44974703e+007 1.07288531e+008 0.00000000e+000]
[0.00000000e+000 0.00000000e+000 2.39203287e+198]]
[D] =
[[60.771 7.663 4.019]
[ 7.659 12.322 4.019]
[ 4.019 4.019 9.567]]
可以看出,第四次迭代最终生成了正确的NM矩阵;然而,[A]矩阵从正确的值变成了错误的A[3][3]元素。然而,在多次尝试重新运行程序后,程序最终输出所有三个正确的矩阵。如前所述,这可能需要1到15次尝试。
在这一点上,我不知道是什么导致了这些问题。我把我的代码给了其他人,让他们在自己的电脑上运行,这个程序似乎没有问题。任何关于如何解决这个问题的建议都将不胜感激。查看第一次迭代的完整代码和说明如下。如果还需要什么,请告诉我。
重新创造
-
运行代码,确保C.stress_critical_loads()在main方法中被注释掉,而C.tsai_wu()没有被注释掉。
-
运行后,滚动到输出的顶部,然后向下滚动,直到出现NM、[A]、[B]和[D]的第一个实例。
-
如果在输出代码中间看到错误消息,那是因为NM矩阵是错误的。当NM矩阵正确时,错误信息消失。尽管有错误消息,该程序也会运行。
再次感谢任何人提供的建议。
import numpy as np
import math
from tabulate import *
def stress_transformation_matrix(x):
# This function calculates the stress transformation
# matrix given an input value for the ply angle
return np.array([
[math.cos(x)**2, math.sin(x)**2, 2*math.sin(x)*math.cos(x)],
[math.sin(x)**2, math.cos(x)**2, -2*math.sin(x)*math.cos(x)],
[-math.sin(x)*math.cos(x), math.sin(x)*math.cos(x), math.cos(x)**2 - math.sin(x)**2]
])
def strain_transformation_matrix(x):
# This function calculates the strain transformation
# matrix given an input value for the ply angle
return np.array([
[math.cos(x)**2, math.sin(x)**2, math.sin(x)*math.cos(x)],
[math.sin(x)**2, math.cos(x)**2, -math.sin(x)*math.cos(x)],
[-2*math.sin(x)*math.cos(x), 2*math.sin(x)*math.cos(x), math.cos(x)**2-math.sin(x)**2]
])
def local_stiffness_matrix(E11, E22, G12, v12, v21):
# This function calculates the local stiffness
# given the material properties at each ply
return np.array([
[E11 / (1 - v12 * v21), v21 * E11 / (1 - v12 * v21), 0],
[v12 * E22 / (1 - v12 * v21), E22 / (1 - v12 * v21), 0],
[0, 0, G12]
])
class Properties:
def __init__(self):
# Material Properties
self.M_E11 = np.array([181E9, 204E9, 138E9, 38.6E9, 76E9])
self.M_E22 = np.array([10.3E9, 18.5E9, 9E9, 8.3E9, 5.5E9])
self.M_G12 = np.array([7.17E9, 5.59E9, 7.1E9, 4.14E9, 2.3E9])
self.M_v12 = np.array([0.280, 0.230, 0.300, 0.260, 0.340])
self.M_v21 = np.array([0.016, 0.021, 0.019, 0.056, 0.33])
self.E11 = np.array([])
self.E22 = np.array([])
self.G12 = np.array([])
self.v12 = np.array([])
self.v21 = np.array([])
# Ultimate Failure Strengths
self.P_SLT = np.array([1500E6, 1260E6, 1447E6, 1062E6, 1400E6])
self.P_SLc = np.array([1500E6, 2499E6, 1447E6, 610E6, 235E6])
self.P_STt = np.array([40E6, 61E6, 52E6, 31E6, 12E6])
self.P_STc = np.array([246E6, 202E6, 206E6, 118E6, 53E6])
self.P_SLTs = np.array([68E6, 67E6, 93E6, 72E6, 34E6])
self.SLT = np.array([])
self.SLc = np.array([])
self.STt = np.array([])
self.STc = np.array([])
self.SLTs = np.array([])
# Ultimate Failure Strains
self.P_epsLtf = np.array([1.087E-2, 0, 1.380E-2, 2.807E-2, 0])
self.P_epsLcf = np.array([0.652E-2, 0, 1.175E-2, 1.754E-2, 0])
self.P_epsCtf = np.array([0.245E-2, 0, 0.436E-2, 0.456E-2, 0])
self.P_epsCcf = np.array([1.818E-2, 0, 2E-2, 1.2E-2, 0])
self.P_epsLTs = np.array([4E-2, 0, 2E-2, 4E-2, 0])
self.epsLtf = np.array([])
self.epsLcf = np.array([])
self.epsCtf = np.array([])
self.epsCcf = np.array([])
self.epsLTs = np.array([])
# Strength Ratios for stress and strain
self.R11_s = np.array([])
self.R22_s = np.array([])
self.R12_s = np.array([])
self.R11_e = np.array([])
self.R22_e = np.array([])
self.R12_e = np.array([])
self.R_crit_s = np.array([])
self.R_crit_e = np.array([])
# Tsai-Wu Coefficients
self.F11 = np.array([])
self.F22 = np.array([])
self.F12 = np.array([])
self.F66 = np.array([])
self.F1 = np.array([])
self.F2 = np.array([])
# Stiffness and Transformation Matrices
self.Q = np.array([])
self.Q_hat = np.array([])
self.T = np.array([])
self.T_hat = np.array([])
# ABD Matrices
self.A = np.empty(shape=(3, 3))
self.B = np.empty(shape=(3, 3))
self.D = np.empty(shape=(3, 3))
self.ABD = np.empty(shape=(6, 6))
# Laminate Loads
self.Nxx = 200
self.Nyy = 200
self.Nxy = 0
self.Mxx = 1
self.Myy = 1
self.Mxy = 0
self.N_mech = np.array([])
self.M_mech = np.array([])
self.NT = np.empty(shape=(3, 1))
self.MT = np.empty(shape=(3, 1))
self.NM = np.array([])
self.Nxx_cs = np.array([])
self.Mxx_cs = np.array([])
self.Nxx_ce = np.array([])
self.Mxx_ce = np.array([])
# Mid-plane Strains and Curvatures
self.e0_k = np.array([])
self.e0 = np.array([])
self.k0 = np.array([])
# Stresses and Strains
self.sg = np.array([]) # global stress
self.sl = np.array([]) # local stress
self.eg = np.array([]) # global strain
self.el = np.array([]) # local strain
# Ply Orientations
self.lam = None # total laminate thickness
self.z_lam = np.array([]) # laminate thickness from to pto bottom
self.z_m = np.array([]) # mid-plane laminate thickness
# Table and Material Properties
self.mat = np.array([])
self.mat_list = np.array([])
# Total Ply Failure
self.Beta = np.array([])
# Angle at Each Ply (deg)
self.angle = np.array([0, 0, 45, 45, -45, -45, 90, 90, -45, -45, 45, 45, 0, 0])
self.a = np.array([])
# Thickness at Each Ply (m)
self.z = np.array([0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3,
0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3, 0.125E-3])
# Number of Ply
self.n = self.angle.size
# Material at Each Ply
self.m = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
# 0: T300/5208
# 1: B(4) 5505
# 2: AS4/3501
# 3: Scotchply 1002
# 4: Kevlar49/Epoxy
# Thermal Properties
self.a1 = -0.018E-6 # longitudinal thermal expansion, deg. C
self.a2 = 24.3E-6 # transverse thermal expansion, deg. C
self.a3 = 0 # shear thermal expansion, deg. C
self.dT = 100 # change in temperature, deg. C
def ply_angle(self):
a = np.array([])
for i in range(self.n):
a = np.append(a, math.radians(self.angle[i]))
self.a = a
def ply(self):
# set total laminate ply failure array
self.Beta = np.full(self.n, 1)
# Determine the thickness of each ply
self.lam = sum(self.z)
print('Total Thickness =', self.lam)
for i in range(self.n):
self.z_lam = np.append(self.z_lam, -self.lam / 2 + (i + 1) * self.z[i])
print('Top-Bottom =', self.z_lam)
for i in range(self.n):
self.z_m = np.append(self.z_m, -self.lam / 2 - self.z[i] / 2 + (i + 1) * self.z[i])
print('Mid-plane =', self.z_m)
def material(self):
for i in range(self.n):
m = self.m[i]
# Material Properties
self.E11 = np.append(self.E11, self.M_E11[m])
self.E22 = np.append(self.E22, self.M_E22[m])
self.G12 = np.append(self.G12, self.M_G12[m])
self.v12 = np.append(self.v12, self.M_v12[m])
self.v21 = np.append(self.v21, self.M_v21[m])
self.mat_list = np.append(self.mat_list, m)
# Ultimate Failure Strengths
self.SLT = np.append(self.SLT, self.P_SLT[m])
self.SLc = np.append(self.STc, self.P_SLc[m])
self.STt = np.append(self.STt, self.P_STt[m])
self.STc = np.append(self.STc, self.P_STc[m])
self.SLTs = np.append(self.SLTs, self.P_SLTs[m])
# Ultimate Failure Strains
self.epsLtf = np.append(self.epsLtf, self.P_epsLtf[m])
self.epsLcf = np.append(self.epsLcf, self.P_epsLcf[m])
self.epsCtf = np.append(self.epsCtf, self.P_epsCtf[m])
self.epsCcf = np.append(self.epsCcf, self.P_epsCcf[m])
self.epsLTs = np.append(self.epsLTs, self.P_epsLTs[m])
# Tsai-Wu Coefficients
self.F11 = np.append(self.F11, 1 / (self.P_SLT[m] * self.P_SLc[m]))
self.F22 = np.append(self.F22, 1 / (self.P_STt[m] * self.P_STc[m]))
self.F12 = np.append(self.F12, -math.sqrt(self.F11[i] * self.F22[i]) / 2)
self.F66 = np.append(self.F66, 1 / self.P_SLTs[m] ** 2)
self.F1 = np.append(self.F1, 1 / self.P_SLT[m] - 1 / self.P_SLc[m])
self.F2 = np.append(self.F2, 1 / self.P_STt[m] - 1 / self.P_STc[m])
if m == 0: self.mat = np.append(self.mat, 'T300/5208')
elif m == 1: self.mat = np.append(self.mat, 'B(4)/5505')
elif m == 2: self.mat = np.append(self.mat, 'AS4/3501')
elif m == 3: self.mat = np.append(self.mat, 'Scotchply 1002')
else: self.mat = np.append(self.mat, 'Kevlar49/Epoxy')
print('F11 = ', self.F11)
print('F22 =', self.F22)
print('F12 =', self.F12)
print('F66 =', self.F66)
print('F1 =', self.F1)
print('F2 =', self.F2)
def mechanical_loads(self):
self.N_mech = np.array([
[self.Nxx],
[self.Nyy],
[self.Nxy]
])
self.M_mech = np.array([
[self.Mxx],
[self.Myy],
[self.Mxy]
])
def local_thermal_coef(self):
al = np.array([
[self.a1],
[self.a2],
[self.a3]
])
return al
def material_properties(self):
# Creates a table for the material properties of each different ply
mat_list = list(map(int, set(self.mat_list)))
for i in mat_list:
data = np.array([[self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i]]])
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
print('\nMaterial =', self.mat[i], '(in SI units)')
print(tabulate(data, headers=['E11', 'E22', 'G12', 'v12', 'v21']))
print('\nMaterial =', self.mat[i], '(in Pa)')
print('[Q] =\n', np.round(Q, 3))
print('')
def table(self):
# Creates a table of ply number, thickness, orientation and material
table = np.array([])
for i in range(self.n):
data = np.array([[i+1, self.mat[i], self.z[i], math.degrees(self.a[i])]])
if i == 0:
table = np.append(table, data)
else:
table = np.vstack([table, data])
print(tabulate(table, headers=['Lamina', 'Material', 'Thickness', 'Orientation']))
class Calculations(Properties):
def __init__(self):
super().__init__()
self.NM = np.array([]) # Loads and Moments
self.e0_k = np.array([]) # Mid-plane strains and curvatures
self.fail_order = np.array([]) # Ply order of failure
def a_matrix(self):
for i in range(self.n):
T = stress_transformation_matrix(self.a[i])
T_hat = strain_transformation_matrix(self.a[i])
if self.Beta[i] == 1:
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
else:
Q = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]
])
Q_hat = np.matmul(np.linalg.inv(T), Q)
Q_hat = np.matmul(Q_hat, T_hat)
A = Q_hat * (self.z_lam[i] - (self.z_lam[i] - self.z[i]))
self.A = np.add(self.A, A)
# print('\nLamina', i, '\t\u03B8 =', math.degrees(self.a[i]), 't =', self.z[i], 'Material =', self.mat[i-1])
# print('[Q\u0302] =\n', np.round(Q_hat, 3))
print('\nLaminate Stiffness Matrices\n[A] =\n', np.round(self.A, 3))
def b_matrix(self):
for i in range(self.n):
T = stress_transformation_matrix(self.a[i])
T_hat = strain_transformation_matrix(self.a[i])
if self.Beta[i] == 1:
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
else:
Q = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]
])
Q_hat = np.matmul(np.linalg.inv(T), Q)
Q_hat = np.matmul(Q_hat, T_hat)
B = Q_hat * (self.z_lam[i] ** 2 - (self.z_lam[i] - self.z[i]) ** 2)
self.B = np.add(self.B, B)
self.B /= 2
print('\n[B] =\n', np.round(self.B, 3))
def d_matrix(self):
for i in range(self.n):
T = stress_transformation_matrix(self.a[i])
T_hat = strain_transformation_matrix(self.a[i])
if self.Beta[i] == 1:
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
else:
Q = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]
])
Q_hat = np.matmul(np.linalg.inv(T), Q)
Q_hat = np.matmul(Q_hat, T_hat)
D = Q_hat * (self.z_lam[i] ** 3 - (self.z_lam[i] - self.z[i]) ** 3)
self.D = np.add(self.D, D)
self.D /= 3
print('\n[D] =\n', np.round(self.D, 3))
def matrices_combined(self):
self.a_matrix() # call A matrix method
self.b_matrix() # call B matrix method
self.d_matrix() # call D matrix method
self.ABD = np.array([
[np.array(self.A[0][0]), np.array(self.A[0][1]), np.array(self.A[0][2]), np.array(self.B[0][0]), np.array(self.B[0][1]), np.array(self.B[0][2])],
[np.array(self.A[1][0]), np.array(self.A[1][1]), np.array(self.A[1][2]), np.array(self.B[1][0]), np.array(self.B[1][1]), np.array(self.B[1][2])],
[np.array(self.A[2][0]), np.array(self.A[2][1]), np.array(self.A[2][2]), np.array(self.B[2][0]), np.array(self.B[2][1]), np.array(self.B[2][2])],
[np.array(self.B[0][0]), np.array(self.B[0][1]), np.array(self.B[0][2]), np.array(self.D[0][0]), np.array(self.D[0][1]), np.array(self.D[0][2])],
[np.array(self.B[1][0]), np.array(self.B[1][1]), np.array(self.B[1][2]), np.array(self.D[1][0]), np.array(self.D[1][1]), np.array(self.D[1][2])],
[np.array(self.B[2][0]), np.array(self.B[2][1]), np.array(self.B[2][2]), np.array(self.D[2][0]), np.array(self.D[2][1]), np.array(self.D[2][2])],
])
def calculated_loads(self):
self.NM = np.matmul(self.ABD, self.e0_k)
print('\nResulting Loads and Moments\n', np.round(self.NM, 3))
def thermal_loads(self):
al = self.local_thermal_coef()
for i in range(self.n):
T = stress_transformation_matrix(self.a[i])
T_hat = strain_transformation_matrix(self.a[i])
if self.Beta[i] == 1:
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
else:
Q = np.array([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]
])
Q_hat = np.matmul(np.linalg.inv(T), Q)
Q_hat = np.matmul(Q_hat, T_hat)
ag = np.matmul(np.linalg.inv(T_hat), al)
NT = np.matmul(Q_hat, ag) * (self.z_lam[i] - (self.z_lam[i] - self.z[i]))
self.NT = np.add(self.NT, NT)
MT = np.matmul(Q_hat, ag) * (self.z_lam[i] ** 2 - (self.z_lam[i] - self.z[i]) ** 2)
self.MT = np.add(self.MT, MT)
self.NT *= self.dT
self.MT *= (self.dT / 2)
def combined_loads(self):
# Combine the mechanical and thermal loads
self.mechanical_loads()
self.thermal_loads()
N = self.N_mech + self.NT
M = self.M_mech + self.MT
self.NM = np.vstack((N, M))
print('\nNM =', self.NM)
def calculated_midplane_strains_curvatures(self):
self.matrices_combined() # call ABD matrix method
# self.loads()
self.e0_k = np.matmul(np.linalg.inv(self.ABD), self.NM)
# print('\nResulting Mid-Plane Strains and Curvatures\n', self.e0_k)
self.e0 = np.array([
self.e0_k[0],
self.e0_k[1],
self.e0_k[2]
])
self.k0 = np.array([
self.e0_k[3],
self.e0_k[4],
self.e0_k[5]
])
def calculated_stresses_and_strains(self):
self.combined_loads() # call calculated loads method
self.calculated_midplane_strains_curvatures() # call mid-plane strains and curvature method
for i in range(self.n):
eg = np.array([self.e0 + self.z_m[i] * self.k0])
if i == 0:
self.eg = np.array(eg)
else:
self.eg = np.append(self.eg, eg, axis=0)
for i in range(self.n):
T_hat = strain_transformation_matrix(self.a[i])
el = np.matrix(np.matmul(T_hat, self.eg[i]))
if i == 0:
self.el = np.array(el)
else:
self.el = np.append(self.el, el, axis=1)
# All values of local strains are 1x3 arrays, so
# each time the strain is needed, the specific
# array must be selected and then transposed
self.el = np.transpose(self.el)
# print('[Local Strain] =', self.el)
for i in range(self.n):
Q = local_stiffness_matrix(self.E11[i], self.E22[i], self.G12[i], self.v12[i], self.v21[i])
sl = np.matrix(np.matmul(Q, np.transpose(np.matrix(self.el[i]))))
if i == 0:
self.sl = np.array(sl)
else:
self.sl = np.append(self.sl, sl, axis=1)
# All values of local stress are 1x3 arrays, so
# each time the stress is needed, the specific
# array must be selected and then transposed
self.sl = np.transpose(self.sl)
print('\n[Local Stress] =\n', self.sl)
def stress_strength_ratio_11(self):
for i in range(self.n):
sl = np.transpose(self.sl[i])
if float(sl[0]) >= 0:
self.R11_s = np.append(self.R11_s, self.SLT[i] / float(sl[0]))
else:
self.R11_s = np.append(self.R11_s, self.SLc[i] / abs(float(sl[0])))
def stress_strength_ratio_22(self):
for i in range(self.n):
sl = np.transpose(self.sl[i])
if float(sl[1]) >= 0:
self.R22_s = np.append(self.R22_s, self.STt[i] / float(sl[1]))
else:
self.R22_s = np.append(self.R22_s, self.STc[i] / abs(float(sl[1])))
def stress_strength_ratio_12(self):
for i in range(self.n):
sl = np.transpose(self.sl[i])
self.R12_s = np.append(self.R12_s, self.SLTs[i] / abs(float(sl[2])))
def stress_critical_loads(self):
self.calculated_stresses_and_strains()
self.stress_strength_ratio_11() # call R11 method
self.stress_strength_ratio_22() # call R22 method
self.stress_strength_ratio_12() # call R12 method
for i in range(self.n):
R = np.array([
[self.R11_s[i], self.R22_s[i], self.R12_s[i]]
])
if self.Beta[i] == 1:
R_crit = np.min(R)
else:
R_crit = 9999
self.R_crit_s = np.append(self.R_crit_s, R_crit)
self.Nxx_cs = np.append(self.Nxx_cs, R_crit * self.Nxx)
self.Mxx_cs = np.append(self.Mxx_cs, R_crit * self.Mxx)
Nxx_c = np.min(self.Nxx_cs)
N_pos = np.where(self.Nxx_cs == Nxx_c)[0]
Mxx_c = np.min(self.Mxx_cs)
self.fail_order = np.append(self.fail_order, N_pos + 1)
print('\nPly Failure in ply', N_pos + 1)
print('R_crit =', np.min(self.R_crit_s))
print('Critical in-plane loading at ply', N_pos + 1, '=', Nxx_c, '')
print('Critical bending loading at ply', N_pos + 1, '=', Mxx_c, '')
print('Stress: Nxx_c =', Nxx_c)
print('Stress: Mxx_c =', Mxx_c)
Beta = np.full(self.n, 0)
if not np.array_equal(self.Beta, Beta):
self.Beta[N_pos] = 0
print(self.Beta)
# Reset all values
self.A = np.zeros(shape=(3, 3))
self.B = np.zeros(shape=(3, 3))
self.D = np.zeros(shape=(3, 3))
self.ABD = np.zeros(shape=(6, 6))
self.eg = np.zeros([])
self.el = np.zeros([])
self.sl = np.zeros([])
self.NT = np.zeros(shape=(3, 1))
self.MT = np.zeros(shape=(3, 1))
self.NM = np.zeros(shape=(6, 1))
self.R11_s = np.zeros([])
self.R22_s = np.zeros([])
self.R12_s = np.zeros([])
self.R_crit_s = np.zeros([])
self.Nxx_cs = np.zeros([])
self.Mxx_cs = np.zeros([])
# re-run calculations
if not np.array_equal(self.Beta, Beta):
self.stress_critical_loads()
else:
print('Ply Fail Order =', self.fail_order)
def tsai_wu(self):
self.calculated_stresses_and_strains()
a = np.array([])
b = np.array([])
c = np.array([])
R1 = np.array([])
R2 = np.array([])
for i in range(self.n):
sl = np.transpose(self.sl[i])
a = np.append(a, self.F11[i] * sl[0] ** 2 + 2 * self.F12[i] * sl[0] * sl[1]
+ self.F22[i] * sl[1] ** 2 + self.F66[i] * sl[2] ** 2)
b = np.append(b, self.F1[i] * sl[0] + self.F2[i] * sl[1])
c = np.append(c, -1)
for i in range(self.n):
if self.Beta[i] == 1:
R1 = np.append(R1, (-b[i] + math.sqrt(b[i]**2 - 4 * a[i] * c[i])) / (2 * a[i]))
R2 = np.append(R2, (-b[i] - math.sqrt(b[i]**2 - 4 * a[i] * c[i])) / (2 * a[i]))
else:
R1 = np.append(R1, 9999)
R2 = np.append(R2, -9999)
R_crit = np.min(R1)
R_crit_pos = np.where(R1 == R_crit)[0]
# print('R1 =', R1)
# print('R2 =', R2)
print('\nR_cr at ply', R_crit_pos + 1, '=', R_crit, '')
self.fail_order = np.append(self.fail_order, R_crit_pos + 1)
NTW_xxc = R_crit * self.Nxx
MTW_xxc = R_crit * self.Mxx
print('N_xx,cr for ply', R_crit_pos + 1, '=', NTW_xxc)
print('M_xx,cr for ply', R_crit_pos + 1, '=', MTW_xxc)
Beta = np.full(self.n, 0)
if not np.array_equal(self.Beta, Beta):
self.Beta[R_crit_pos] = 0
print('\nFailure Progression =', self.Beta)
# Reset all values
self.A = np.zeros(shape=(3, 3))
self.B = np.zeros(shape=(3, 3))
self.D = np.zeros(shape=(3, 3))
self.ABD = np.zeros(shape=(6, 6))
self.eg = np.zeros([])
self.el = np.zeros([])
self.sl = np.zeros([])
self.NT = np.zeros(shape=(3, 1))
self.MT = np.zeros(shape=(3, 1))
self.NM = np.zeros(shape=(6, 1))
# re-run all values
if not np.array_equal(self.Beta, Beta):
self.tsai_wu()
else:
print('Ply Fail Order =', self.fail_order)
def main():
C = Calculations()
C.ply_angle()
C.ply()
C.material()
C.material_properties()
C.table()
# Comment out for Tsai-Wu failure calculations
# print('\nCritical Stress Criterion')
# C.stress_critical_loads()
# Comment out for Maximum Stress failure calculations
print('\nTsai-Wu Criterion\n')
C.tsai_wu()
if __name__ == '__main__':
main()
