本篇代码功能和上一篇一样，只是优化了挠度函数，
到
代码语言：javascript
复制
from matplotlib import pyplot as pltimport numpy as np
import sympy
from scipy.optimize import fsolve
class RoundPlate(object):
"""四周固定圆板在均布载荷下大挠度弯曲求解"""
def __init__(self, E: float, nu: float, a: float, t: float):
"""弹性模量，泊松比，半径，板厚"""
self.E = E # 弹性模量
self.nu = nu # 泊松比
self.D = E * t ** 3 / (12.0 * (1 - nu ** 2)) # 板的抗弯刚度
self.g_nu = (-2791.0 * nu * nu + 4250.0 * nu + 7505.0) / 17640.0 # 仅与泊松比相关的系数
self.a = a # 圆板的半径
self.t = t # 板的厚度
self.k = E * t / (180180.0 * a ** 2 * (1 - nu ** 2))
self.k2 = 16.0 * self.D / (15.0 * self.a ** 2)
# 下面的参数与载荷 q 有关
self.c0 = 0.0
self.c1 = 0.0
self.a0 = 0.0
self.a1 = 0.0
def w0_linear(self, q: float) -> float: #
"""求线性条件下的板最大位移（圆心处位移），q 为均布载荷（压强）"""
return (self.a ** 4) * q / 64.0 / self.D
def f_w0(self, w0: float, q: float) -> float:
"""求非线性条件下的板最大位移w0（圆心处位移）的方程 f_w0(w0)=0"""
return w0 + self.g_nu / self.t ** 2 * w0 ** 3 - self.w0_linear(q)
def d_f_w0(self, w0: float) -> float:
"""d f_w0(w0)/d(w0), f_w0的一阶导数，for牛顿迭代法"""
return 3.0 * self.g_nu / self.t ** 2 * w0 * w0 + 1.0
def solve(self, q: float, error: float = 1e-10) -> float:
"""牛顿迭代法求解一元方程。q 为均布载荷（压强），error 为求解精度"""
x = self.w0_linear(q) # 用线性条件下的板最大位移 做为迭代的初值
i = 0
while True:
delta_x = self.f_w0(x, q) / self.d_f_w0(x)
x -= delta_x
if delta_x < error * self.t: # 改用相对板厚的误差
break
i += 1
# print(f"// {i} 次迭代后收敛 //")
return x
def cal_eq3_coefficients(self, q0: float):
"计算二元 (C0,C1)3次方程 eq3 的各向系数"
a = (-3192904.0 * self.nu ** 2 + 4862000.0 * self.nu + 8585720.0) / 63.0
a *= self.k
b = (-1119768.0 * self.nu ** 2 + 1489800.0 * self.nu + 2701920.0) / 7.0
b *= self.k
c = (-13344156.0 * self.nu ** 2 + 15530280.0 * self.nu + 30130932.0) / 77.0
c *= self.k
d = (-4941324.0 * self.nu ** 2 + 5017572.0 * self.nu + 10773432.0) / 77.0
d *= self.k
e = 20.0 * self.k2
f = 15.0 * self.k2
g = -q0 * self.a ** 2 / 3.0
return a, b, c, d, e, f, g
def cal_eq4_coefficients(self, q0: float):
"计算二元 (C0,C1)3次方程 eq4 的各向系数"
a = (-373256.0 * self.nu ** 2 + 496600.0 * self.nu + 900640.0) / 7.0
a *= self.k
b = (-13344156.0 * self.nu ** 2 + 15530280.0 * self.nu + 30130932.0) / 77.0
b *= self.k
c = (-14823972.0 * self.nu ** 2 + 15052716.0 * self.nu + 32320296.0) / 77.0
c *= self.k
d = (-73019070 * self.nu ** 2 + 64726452.0 * self.nu + 157413618.0) / 1001.0
d *= self.k
e = 15.0 * self.k2
f = 18.0 * self.k2
g = -0.25 * q0 * self.a ** 2
return a, b, c, d, e, f, g
def solve2(self, q0, c0_init, c1_init=0):
"用 scipy.optimize.fsolve 求解位移方程的系数 C0, C1"
a1, b1, c1, d1, e1, f1, g1 = self.cal_eq3_coefficients(q0)
a2, b2, c2, d2, e2, f2, g2 = self.cal_eq4_coefficients(q0)
def solve_function(init_values):
"""建立方程"""
x, y = init_values
return [
a1 * x ** 3 + b1 * x * x * y + c1 * x * y * y + d1 * y ** 3 + e1 * x + f1 * y + g1,
a2 * x ** 3 + b2 * x * x * y + c2 * x * y * y + d2 * y ** 3 + e2 * x + f2 * y + g2,
]
init_values = c0_init, c1_init
x, y = fsolve(solve_function, init_values)
return x, y
def w(self, rho: float, q0: float) -> float:
"""计算板中面任一点的位移"""
c0_init = self.solve(q0)
c1_init = 0.0
self.c0, self.c1 = self.solve2(q0, c0_init, c1_init)
return self.c0 * (1 - (rho / self.a) ** 2) ** 2 + self.c1 * (1 - (rho / self.a) ** 2) ** 3
def cal_a0_a1(self):
"""计算水平径向位移函数的系数"""
self.a0 = -(25454.0 * self.c0 **2 * self.nu + 61308.0 * self.c0 * self.c1 * self.nu
+ 37503.0 * self.c1**2 * self.nu - 51194.0 * self.c0**2 - 93600.0 * self.c0 * self.c1
- 44685.0 * self.c1**2) / (36036.0 * self.a)
self.a1 = (3718.0 * self.c0**2 * self.nu + 12636.0*self.c0*self.c1*self.nu + 9099.0 * self.c1**2 * self.nu
- 22594.0 * self.c0**2 - 46800.0 * self.c0 * self.c1 - 24273.0 * self.c1**2) / (12012.0 * self.a)
def u(self, rho: float) -> float:
"""计算水平径向位移分量"""
self.cal_a0_a1()
return (1.0 - rho/self.a) * (self.a0 * (rho/self.a) + self.a1 * (rho/self.a)**2)
def epsilon_rho_mid(self, rho: float) -> float: # 中面径向应变
# 先求中面上的（z=0）的应变
ex = (1.0 - rho/self.a)*(self.a0/self.a + 2*self.a1*rho/(self.a**2)) - (self.a0*rho/self.a**2 + self.a1 * rho**2 / self.a**3)
temp = 1.0 - (rho/self.a)**2
dw_drho = -(4 * self.c0 * temp + 6 * self.c1 * temp**2) * rho / self.a**2
ex += 0.5 * dw_drho**2
return ex
def epsilon_phi_mid(self, rho: float) -> float: # 中面周向应变 = u(rho)/rho
return self.u(rho) / rho
# def gama_xy_mid(self, rho: float) -> float: # 中面切应变, 取d W(rho)/ d rho
# temp1 = rho / self.a**2
# temp2 = 1.0 - (rho/self.a)**2
# return -4.0*self.c0 * temp1 * temp2 - 6.0 * self.c1 * temp1 * temp2**2
#
# def epsilon_rho(self, rho: float, z: float) -> float: # 径向应变
# ex = self.epsilon_rho_mid(rho)
# # w 对 rho 的二阶导数
# dw_dz_2 = -30 * self.c1 * rho**4/self.a**6 + 12.0 * rho**2 * (self.c0 + 3*self.c1)/self.a**4 - (4.0 * self.c0 + 6.0 * self.c1)/self.a**2
# ex += -dw_dz_2 * z
# return ex
#
# def sigma_top(self,rho:float):
# temp1 = self.E / (1.0 - self.nu**2)
# sigma_rho = temp1 * (self.epsilon_rho(rho, -0.5 * self.t) + self.nu * self.epsilon_phi_mid(rho))
# sigma_phi = temp1 * (self.epsilon_phi_mid(rho) + self.nu * self.epsilon_rho(rho, -0.5 * self.t))
# gama = self.E/(2.0* (1.0+self.nu)) * self.gama_xy_mid(rho)
# return sigma_rho, sigma_phi, gama
#
# def sigma_bottom(self,rho: float):
# temp1 = self.E / (1.0 - self.nu**2)
# sigma_rho = temp1 * (self.epsilon_rho(rho, +0.5 * self.t) + self.nu * self.epsilon_phi_mid(rho))
# sigma_phi = temp1 * (self.epsilon_phi_mid(rho) + self.nu * self.epsilon_rho(rho, +0.5 * self.t))
# gama = self.E/(2.0*(1.0+self.nu)) * self.gama_xy_mid(rho)
# return sigma_rho, sigma_phi, gama
if name == "main":
# 模型输入
E = 200e9 # 弹性模量
nu = 0.3 # 泊松比。范围 -1< nu<= 0.5
a = 1.0e-3 # 圆板的半径
t = 1.0e-6 # 板的厚度
rp = RoundPlate(E, nu, a, t)
q0 = 10.0 # 分布载荷
print(f"抗弯刚度 D = {rp.D}")
print(f"线性条件下的板最大位移 w0_L = {rp.w0_linear(q0)}")
c0_init = rp.solve(q0)
print(f"非线性(w是rho的4次函数）条件下板最大位移 w0 = {c0_init}")
c0, c1 = rp.solve2(q0, c0_init)
print(f"C0 ={c0}, C1={c1}")
print(f"非线性((w是rho的6次函数)条件下板最大位移 w0 = {rp.w(0, q0)}")
r = a/2
ex = rp.epsilon_rho_mid(r)
print(f"ex = {ex} at rho= {r}")
r = 0
ex = rp.epsilon_rho_mid(r)
print(f"ex = {ex} at rho= {r}")
r = a
ex = rp.epsilon_rho_mid(r)
print(f"ex = {ex} at rho= {r}")
Q = np.linspace(0, q0, 101) # 载荷
# print(Q)
W0_L = np.frompyfunc(rp.w0_linear, 1, 1)(Q) # 线性解
# print(W0_L)
W0_4 = np.frompyfunc(rp.solve, 2, 1)(Q, 1e-10) # 非线性(w是rho的4次函数）
# print(W0_4)
W0_6 = np.frompyfunc(rp.w, 2, 1)(0, Q) # 非线性((w是rho的6次函数)
# print(W0_6)
# 读取专业 FEM 软件在相同 几何、材料和载荷下的轴对称模型的计算结果
FEMdata = np.loadtxt(r"e:\disp.txt", dtype=float, comments='#', delimiter=" ", converters=None, skiprows=1, usecols=None)
print(FEMdata)
plt.title("四周固定的薄圆板在均布载荷下大挠度弯曲最大位移计算")
plt.plot(Q, W0_L, c='r', linestyle="--", label="线性条件下最大位移")
plt.scatter(Q, W0_L, c='r')
plt.plot(Q, W0_4, c='b', linestyle="-", label="非线性(w是rho的4次函数）条件下最大位移")
plt.scatter(Q, W0_4, c='b')
plt.plot(Q, W0_6, c='g', linestyle="-", label="非线性(w是rho的6次函数，更精确）条件下最大位移")
plt.plot(FEMdata[:, 0], FEMdata[:, 1]/1e6, c='cyan', linestyle="-", label="专业FEM软件(开启几何非线性)算得的最大位移") # /1e6 , 将微米转为米
plt.scatter(Q, W0_6, c='g')
plt.xlabel("均布载荷q /Pa")
plt.ylabel("最大位移（中心位移）/m")
plt.legend(loc="lower right")
plt.show()
# print(f"中心点位移：{rp.w(0, q0)}")
# print(f"边缘位移：{rp.w(rp.a, q0)}")
# print(Q.max())
R = np.linspace(-rp.a, rp.a, 101)
W = -np.frompyfunc(rp.w, 2, 1)(np.abs(R), Q.max())
plt.title("四周固定的圆板在均布载荷下大挠度弯曲 挠曲面")
plt.plot(R, W, c='b')
plt.xlabel("坐标 r")
plt.ylabel("垂向位移")
plt.show()
U = np.frompyfunc(rp.u, 1, 1)(np.abs(R))
plt.title("四周固定的圆板在均布载荷下大挠度弯曲 中面径向位移")
plt.plot(R, U, c='b')
plt.xlabel("坐标 r")
plt.ylabel("中面径向位移")
plt.show()
# E_rho_mid = np.frompyfunc(rp.epsilon_rho_mid,1,1)(np.abs(R)) # 中面径向应变
# Epsilon_rho_top = np.frompyfunc(rp.epsilon_rho,2,1)(np.abs(R), -t/2) # 顶面径向应变
# Epsilon_rho_bottom = np.frompyfunc(rp.epsilon_rho,2,1)(np.abs(R), +t/2) # 底面（z向向下）径向应变
# Gama_xy_mid = np.frompyfunc(rp.gama_xy_mid,1,1)(np.abs(R)) # 中面切应变
# plt.title("四周固定的圆板在均布载荷下大挠度弯曲 表面径向应变")
# plt.plot(R, E_rho_mid, c='g', label="epsilon_rho_rho mid")
# plt.plot(R, Epsilon_rho_top, c='r', label="epsilon_rho_rho top")
# plt.plot(R, Epsilon_rho_bottom, c='b', label="epsilon_rho_rho bottom")
# # plt.plot(R, Gama_xy_mid, c='purple', label="Gama_mid")
# plt.xlabel("坐标 r")
# plt.ylabel("表面应变")
# plt.legend(loc="lower right")
# plt.show()
#
#
# E_phi_mid = np.frompyfunc(rp.epsilon_phi_mid, 1, 1)(np.abs(R)) # 中面周向应变
# plt.title("四周固定的圆板在均布载荷下大挠度弯曲 中面周向应变")
# plt.plot(R, E_phi_mid, c='b')
# plt.xlabel("坐标 r")
# plt.ylabel("中面周向应变")
# plt.show()
#
# Sigma_rho_bot, Sigma_phi_bot, Gama_bot = np.frompyfunc(rp.sigma_bottom, 1, 3)(np.abs(R)) # 中面周向应变
# plt.title("四周固定的圆板在均布载荷下大挠度弯曲 底面应力")
# plt.plot(R, Sigma_rho_bot, c='r', label="Sigma_rho")
# plt.plot(R, Sigma_phi_bot, c='g', label="Sigma_phi")
# # plt.plot(R, Gama_bot, c='b', label="Gama_rho_phi")
# plt.xlabel("坐标 r")
# plt.ylabel("底面应力")
# plt.legend(loc="lower right")
# plt.show()
#
# Sigma_rho_top, Sigma_phi_top, Gama_top = np.frompyfunc(rp.sigma_top, 1, 3)(np.abs(R)) # 中面周向应变
# plt.title("四周固定的圆板在均布载荷下大挠度弯曲 顶面应力")
# plt.plot(R, Sigma_rho_top, c='r', label="Sigma_rho")
# plt.plot(R, Sigma_phi_top, c='g', label="Sigma_phi")
# # plt.plot(R, Gama_bot, c='b', label="Gama_rho_phi")
# plt.xlabel("坐标 r")
# plt.ylabel("顶面应力")
# plt.legend(loc="lower right")
# plt.show()</code></pre></div></div><div class="rno-markdown-code"><div class="rno-markdown-code-toolbar"><div class="rno-markdown-code-toolbar-info"><div class="rno-markdown-code-toolbar-item is-type"><span class="is-m-hidden">代码语言：</span>javascript</div></div><div class="rno-markdown-code-toolbar-opt"><div class="rno-markdown-code-toolbar-copy"><i class="icon-copy"></i><span class="is-m-hidden">复制</span></div></div></div><div class="developer-code-block"><pre class="prism-token token line-numbers language-javascript"><code class="language-javascript" style="margin-left:0"></code></pre></div></div>
