四周固定的圆板在均布载荷下大挠度弯曲的计算【2】

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()