DEFLECTION-BENDING MOMENT COUPLING NEURAL NETWORK METHOD FOR THE BENDING PROBLEM OF THIN PLATES WITH IN-PLANE STIFFNESS GRADIENT
-
摘要: 发展了一种求解面内变刚度功能梯度薄板弯曲问题的神经网络方法. 面内变刚度薄板弯曲问题的偏微分控制方程为一复杂的4阶偏微分方程, 传统的基于强形式的神经网络解法在求解该偏微分方程时可能会遇到难以收敛、边界条件难以处理的情况. 本文基于Kirchhoff薄板弯曲理论, 提出了一种直角坐标系下任意面内变刚度薄板弯曲问题的神经网络解法. 神经网络模型包含挠度网络与弯矩网络, 分别用于预测薄板的挠度与弯矩, 从而将求解4阶偏微分方程转换为求解一系列二阶偏微分方程组, 通过对挠度、弯矩试函数的构造可使得神经网络计算结果严格满足边界条件. 在误差的反向传播中, 根据本文提出的误差函数公式计算训练误差并结合Adam优化算法更新模型的内部参数. 求解了不同边界条件、形状的面内变刚度薄板弯曲问题, 并将所得计算结果与理论解、有限元解进行对比. 研究表明, 本文模型对于求解面内变刚度薄板弯曲问题具备适应性, 虽然模型中的弯矩网络收敛较挠度网络要慢, 但本文方法在试函数的构造上更为简单、适应性更强.Abstract: A neural network method is developed to solve the bending problems of functionally graded thin plates with in-plane stiffness gradient in this paper. The partial differential equation (PDE) of the bending of thin plates with in-plane stiffness gradient is a complex fourth-order PDE. The conventional neural network solution based on strong form, may face the problem of slow convergence and the boundary conditions are difficult to handle when solving the PDE. According to the Kirchhoff thin plate bending theory, a neural network solution to the bending problem of thin plates with any in-plane stiffness gradient in rectangular coordinate system is proposed in this paper. The neural network model includes deflection network and bending moment network, which are used to predict the deflection and bending moment of the thin plate respectively. Thus, the solution of the fourth-order PDE is transformed into a series of second-order PDEs. By constructing trial function of the deflection and bending moment, the results of neural network calculation can be strictly satisfied the boundary conditions. In the back propagation process, training error is calculated according to the error function formula proposed in this paper and combining Adam optimization algorithm to update the internal parameters of the neural networks. In this paper, the bending problems of thin plate with in-plane stiffness gradient with different boundary conditions and shapes are solved, and the calculated results are compared with theoretical solutions or those of finite element solutions. It shows that the proposed method is suitable for solving the bending problem of thin plate with in-plane stiffness gradient. And the convergence of bending moment network is slower than the deflection network. However, it is robust and easier in dealing with boundary conditions.
-
引言
面内功能梯度材料薄板结构在土木工程、海洋工程等领域应用广泛, 该结构由功能梯度材料所组成(functionally graded material, FGM), 其材料特性随着空间位置的变化而表现为梯度性的变化[1]. 在现有研究中, 于天崇等[2]研究了面内变刚度薄板在特定边界下弯曲问题的Levy解. 朱竑祯等[3]研究了周边固支圆形面内变刚度薄板轴对称弯曲问题的级数解答. 何建璋等[4]研究了面内变刚度矩形薄板自由振动问题的辛弹性解. 以上的理论解答仅针对特定的功能梯度函数及特定边界才成立, 一般的情况下难以得出理论解答, 而在数值解法上仍以有限元为主. Santare和Lambros[5]发展了一种针对材料属性为指数分布的梯度有限元求解格式. Kim和Paulino[6]研究了梯度单元以及分层单元在不同荷载下的计算性能. 黄立新等[7]基于分层法思想分析了功能梯度材料的平面应力问题. 田云德和秦世伦[8]采用分层法研究了功能梯度厚板的热应力问题. 对于面内变刚度功能梯度薄板, 采用分层法, 薄板求解域采用有限元网格划分, 每个单元的材料参数为常数, 而其材料参数则根据功能梯度函数由单元内特定点进行计算. 有限元网格划分越密其计算结果越精确, 而在实际计算中, 越精细的网格会导致总体刚度矩阵规模巨大, 需要耗费大量的计算机内存. 无论采用何种数值方法, 其最终目的均是求得面内变刚度薄板弯曲控制偏微分方程的近似解答, 为进一步丰富该类研究, 本文拟结合深度学习技术并发展求解该类问题的新解法.
在早期就有研究[9-10]将人工神经网络作为一类偏微分方程的求解器用于求解偏微分方程, 但由于其对计算机计算能力的要求过高以及优化算法中存在的问题, 这一解法在当时并未得到很好的发展. 而如今, 自深度学习在计算机视觉、语音文字识别取得成功的应用后, 深度学习技术也在各个学科领域加速发展. 在力学领域, Weinan和Yu[11]提出深度Ritz法, 该方法采用变分求解形式对偏微分方程进行求解. Raissi等[12]提出了用于求解高阶非线性偏微分方程的物理驱动的神经网络(physics-informed neural networks, PINNs). Sirignano和Spiliopoulos[13]则提出求解高阶微分方程的深度伽辽金法(deep galerkin method, DGM). Samaniego等[14]建立了深度能量法并将其应用于求解弹性、超弹性等力学问题. 瞿同明等[15]基于深度学习技术, 研究了细观力学中的颗粒本构关系. 谢晨月等[16]发展了一种模拟湍流大涡的神经网络方法. 刘宇翔等[17]基于卷积神经网络研究了无网格方法中影响域的优化问题. 郭宏伟和庄晓莹[18]采用深度配点法以及深度能量法求解了薄板弯曲问题. 陈豪龙和柳占立[19]基于数据驱动的神经网络模型求解了热传导反问题.
在上述研究中, 神经网络解法[11-14]并不像有限元解法一样可以轻松施加边界条件, 早期的研究采取根据边界条件构造满足偏微分方程特解试函数的形式来处理边界条件, 但采用该方法会使得简支边、自由边试函数的表达式变得复杂, 导致程序的实现较为复杂. 近期的研究则采用罚函数的方法将边界处的误差纳入神经网络的训练误差中, 从而将原问题转换为无约束优化问题, 在实际计算中, 也会存在着由于边界误差项难以收敛而影响求解精度的情况[20].
同时由于弯曲刚度函数是面内坐标的连续函数, 面内变刚度薄板弯曲问题的控制方程为一包含了弯曲刚度导数项的复杂4阶偏微分方程, 在实际计算中采用DGM和PINN等方法对其求解时, 会存在由于弯曲刚度偏导数在域内不收敛而导致网络拟合不佳的问题.
基于上述原因, 本文针对薄板弯曲问题求解的特点, 结合前面所述的两种边界处理方案, 建立了一种针对面内变刚度薄板弯曲问题的非全连接前馈神经网络模型, 该模型包含挠度网络与弯矩网络: 挠度网络用于预测薄板的挠度, 弯矩网络用于预测薄板的弯矩, 进而将问题转换为求解4个二阶偏微分方程组. 在边界条件的处理上, 本文仍采用罚函数方法, 不同之处在于本文模型的输出为挠度、弯矩, 因而可根据位移边界条件对挠度网络构造试函数, 根据广义应力边界条件对弯矩网络构造试函数, 这使得本文模型对于常见的边界条件的施加更为简便, 进而减小边界误差项带来的影响, 同时计算效率也得到提高. 本文采用Pytorch深度学习框架编写求解程序, 选取不同边界条件的面内变刚度薄板算例, 在Ubuntu Kylin操作系统上进行计算, 计算机的CPU配置为Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz, 8GB内存, 并将计算所得结果与理论解、有限元解进行对比分析, 以验证本文方法的有效性.
1. 面内变刚度薄板弯曲理论
本文研究变厚度薄板或弹性模量参数在面内变化的薄板的弹性小变形弯曲, 设薄板的厚度函数为
$ h(x,y), $ 材料的泊松比$\nu $ 为常数, 弹性模量函数为$ E(x,y) $ .根据Kirchhoff板理论基本假定, 几何方程为
$$ \left. {\begin{array}{*{20}{l}} {{\varepsilon _x} = - z\dfrac{{{\partial ^2}w(x,y)}}{{\partial {x^2}}}} \\ {{\varepsilon _y} = - z\dfrac{{{\partial ^2}w(x,y)}}{{\partial {y^2}}}} \\ {{\gamma _{xy}} = - z\dfrac{{{\partial ^2}w(x,y)}}{{\partial x\partial y}}} \end{array}} \right\} $$ (1) 物理方程为
$$ \left. {\begin{array}{*{20}{l}} {{\sigma _x} = \dfrac{{E(x,y)}}{{1 - {\nu ^2}}}({\varepsilon _x} + \nu {\varepsilon _y})} \\ {{\sigma _y} = \dfrac{{E(x,y)}}{{1 - {\nu ^2}}}({\varepsilon _y} + \nu {\varepsilon _x})} \\ {{\tau _{xy}} = \dfrac{{E(x,y)}}{{2(1 + \nu )}}{\gamma _{xy}}} \end{array}} \right\} $$ (2) 平衡方程为
$$ \frac{{{\partial ^2}{M_x}}}{{\partial {x^2}}} + 2\frac{{{\partial ^2}{M_{xy}}}}{{\partial x\partial y}} + \frac{{{\partial ^2}{M_y}}}{{\partial {y^2}}} + q = 0 $$ (3) 广义应力−应变关系为
$$ \begin{split} &{M}_{x}={\displaystyle {\int }_{-h/2}^{h/2}\frac{E}{1-{\nu }^{2}}}\left({\epsilon }_{x}+\nu {\epsilon }_{y}\right)z \text{d}z=\\ &\qquad -\frac{E(x,y)h{(x,y)}^{3}}{12\left(1-{\nu }^{2}\right)}\left[\frac{{\partial }^{2}w(x,y)}{\partial {x}^{2}}+\nu \frac{{\partial }^{2}w(x,y)}{\partial {y}^{2}}\right]\end{split} $$ (4) 记弯曲刚度函数为
$$ D(x,y) = \frac{{E\left( {x,y} \right)h{{\left( {x,y} \right)}^3}}}{{12\left( {1 - {\nu ^2}} \right)}} $$ (5) 则
$$ {M_x} = - D(x,y)\left[ {\frac{{{\partial ^2}w(x,y)}}{{\partial {x^2}}} + \nu \frac{{{\partial ^2}w(x,y)}}{{\partial {y^2}}}} \right] $$ (6) 同理
$$ {M_y} = - D(x,y)\left[ {\frac{{{\partial ^2}w(x,y)}}{{\partial {y^2}}} + \nu \frac{{{\partial ^2}w(x,y)}}{{\partial {x^2}}}} \right] $$ (7) $$ {M_{xy}} = \int_{ - h/2}^{h/2} {{\tau _{xy}}} z{\rm{d}}z = - (1 - \nu )D(x,y)\frac{{{\partial ^2}w(x,y)}}{{\partial x\partial y}} $$ (8) 对于固支边界条件
${\varGamma _1}$ $$ {\left( {w,\frac{{\partial w}}{{\partial n}}} \right)_s} = 0 $$ (9) 对于简支边界条件
${\varGamma _2}$ $$ {\left( {w,{M_n}} \right)_s} = 0 $$ (10) 对于自由边界条件
${\varGamma _3}$ $$ \left( {{M_n},F_{{{s}}n}^t} \right) = 0 $$ (11) 其中
$n$ 表示边界的外法线方向,$s$ 表示边界的切线方向.将式(4) ~ 式(6)代入平衡方程(3)即可得面内变刚度薄板弯曲偏微分控制方程
$$ \begin{split} & {\nabla ^2}\left( {D{\nabla ^2}w} \right) - (1 - \nu )\left( {\frac{{{\partial ^2}D}}{{\partial {y^2}}}\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right. - 2\frac{{{\partial ^2}D}}{{\partial x\partial y}}\frac{{{\partial ^2}w}}{{\partial x\partial y}} + \\ & \qquad \left. {\frac{{{\partial ^2}D}}{{\partial {x^2}}}\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right) + q = 0 \end{split} $$ (12) 式中
${\nabla ^2} = \dfrac{{{\partial ^2}}}{{\partial {x^2}}} + \dfrac{{{\partial ^2}}}{{\partial {y^2}}}$ 为Laplace算子.2. 神经网络模型
2.1 神经网络模型的构建
本文方法并非直接设计网络来求解方程(12), 而是采用两个神经网络模型来进行求解, 如图1所示, 将待求解的4阶偏微分控制方程转换为求解4个二阶偏微分方程组, 该解法本质上仍属于强形式的求解方案. 如果仅以挠度作为预测解, 在试函数的构造上对于不同形状的求解域以及简支、自由边界条件的构造会出现困难. 本文采用挠度网络预测薄板挠度
$\bar w\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{w}}}} \right)$ , 弯矩网络预测薄板弯矩,$\bar M\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right) = $ $ \left\{ {{{\bar M}_x}\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right),{{\bar M}_{xy}}\left( {x,y;{{{{\boldsymbol{\theta }}}}_{\rm{m}}}} \right),{{\bar M}_y}\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right)} \right\}$ , 这样的做法可以使得位移边界条件由挠度网络施加, 广义应力边界条件由弯矩网络施加.2.2 误差函数
根据广义应力−应变关系式(6) ~ 式(8)可求得挠度二阶偏导的弯矩表达式
$$ \frac{{{\partial ^2}w}}{{\partial {x^2}}} = \frac{1}{{D\left( {1 - {\nu ^2}} \right)}}\left( {\nu {M_y} - {M_x}} \right) $$ (13) $$ \frac{{{\partial ^2}w}}{{\partial {y^2}}} = \frac{1}{{D\left( {1 - {\nu ^2}} \right)}}\left( {\nu {M_x} - {M_y}} \right) $$ (14) $$ \frac{{{\partial ^2}w}}{{\partial x\partial y}} = - \frac{{{M_{xy}}}}{{D(1 - \nu )}} $$ (15) 误差函数的构造是神经网络训练的核心, 由于本文方法引入两个网络进行计算, 故在训练中需要考虑两者之间的耦合误差. 若采用无约束优化方案, 本文误差函数主要根据挠度与弯矩网络在边界处的误差、弯矩网络在力平衡方程(3)中的误差、预测的挠度与弯矩通过式(13) ~ 式(15)建立的耦合误差来构造.
采用均方误差(mean square error, MSE)来衡量神经网络的拟合误差, 记挠度网络与弯矩网络的内部参数分别为
$ {\boldsymbol{\theta }}_{{\rm{w}}} $ ,${\boldsymbol{\theta }}_{\rm{m}}$ , 本文模型的误差函数可构造为$$ \begin{split} C\left( {{{{{\boldsymbol{\theta}} }}_{\rm{w}}},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right) = & MS{E_{\text{R}}} + {k_{\rm{P}}}MS{E_{\rm{P}}} + \\ &{k_1}MS{E_{{\varGamma _1}}} + {k_2}MS{E_{{\varGamma _2}}} + {k_3}MS{E_{{\varGamma _3}}} \end{split} $$ (16) 如果采用构造试函数的形式使得边界误差强制满足, 则误差函数无需计算边界误差
$$ C\left( {{{{{\boldsymbol{\theta}} }}_w},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right) = MS{E_{\text{R}}} + {k_{\rm{P}}}MS{E_{\rm{P}}} $$ (17) 其中
$$ \begin{split} MS{E_{\rm{R}}} =& \frac{1}{{{N_\varOmega }}}\sum\limits_{i = 1}^{{N_\varOmega }} {\left\| {\frac{{{\partial ^2}{{\bar M}_x}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)}}{{\partial {x^2}}}} \right. + 2\frac{{{\partial ^2}{{\bar M}_{xy}}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)}}{{\partial x\partial y}}} +\\ & {\left. {\frac{{{\partial ^2}{{\bar M}_y}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)}}{{\partial {y^2}}} + q} \right\|^2} \end{split} $$ (18) $$ MS{E_{{\varGamma _1}}} = \frac{1}{{{N_{{\varGamma _1}}}}}\sum\limits_{i = 1}^{{N_{{\varGamma _1}}}} {\left[{{{\left\| {\bar w\left( {{{\boldsymbol{x}}_{{\varGamma _1}}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)} \right\|}^2} + \left. {{{\left\| {\frac{{\partial \bar w\left( {{{\boldsymbol{x}}_{{\varGamma _1}}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)}}{{\partial n}}} \right\|}^2}} \right]} \right.} \;\;\;\;\;\;$$ (19) $$ MS{E_{{\varGamma _2}}} = \frac{1}{{{N_{{\varGamma _2}}}}}\sum\limits_{i = 1}^{{N_{{\varGamma _2}}}} {\left[ {{{\left\| {\bar w\left( {{{\boldsymbol{x}}_{{\varGamma _2}}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)} \right\|}^2} + \left\| {{{\left. {{{\bar M}_n}\left( {{{\boldsymbol{x}}_{{\varGamma _2}}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right) - {{\tilde M}_n}} \right\|}^2}} \right.} \right]} $$ (20) $$ \begin{split} MS{E_{{\varGamma _3}}} = &\frac{1}{{{N_{{\varGamma _3}}}}}\sum\limits_{i = 1}^{{N_{{\varGamma _3}}}} {\Biggl[{23} {{{\left\| {\left. {{{\bar M}_n}\left( {{{\boldsymbol{x}}_{{\varGamma _3}}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right) - {{\tilde M}_n}} \right\|} \right.}^2}} } \hfill +\\ & {{{\left\| {\frac{{\partial {{\bar M}_{ns}}\left( {{{\boldsymbol{x}}_{{\varGamma _3}}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)}}{{\partial s}} + {Q_n}\left( {{{\boldsymbol{x}}_{{\varGamma _3}}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)} \right\|}^2}} \Biggl]{23} \hfill \end{split}\;\;\;\;\;\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ (21) $$ MS{E_{\rm{P}}} = MS{E_{{\rm{P}}1}} + MS{E_{{\rm{P}}2}} + MS{E_{{\rm{P}}3}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; $$ (22) 其中, 向量
${\boldsymbol{x}} = \left( {x,y} \right)$ 表示神经网络的输入; n和s分别为边界的法线、切线方向;$\varOmega $ 表示求解域;$ {{\varGamma }_1} $ ,$ {{\varGamma }_{2}} $ ,$ {{\varGamma }_3} $ 分别为固支、简支、自由边界,$\partial \varOmega = {\varGamma _1} \cup {\varGamma _2} \cup {\varGamma _3}$ ;$ {\tilde M_n} $ 为施加于边界处的弯矩,${M_n}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)$ ,${{{Q}}_n}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)$ 分别为根据弯矩网络输出求得的弯矩、剪力;${k_{\rm{p}}}$ 为网络耦合系数, 取值范围为1 ~ 1000, 该系数的选取会影响挠度网络与弯矩网络之间的耦合效果;${k_1}$ ,${k_2}$ ,${k_3}$ 为边界处的罚系数, 取值范围为1 ~ 10 000.式(22)中
$$ \begin{split} MS{E_{{\rm{P}}1}} =& \frac{1}{{{N_\varOmega }}}\sum\limits_{i = 1}^{{N_\varOmega }} {\left\| {\frac{{{\partial ^2}\bar w\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)}}{{\partial {x^2}}}} \right.} \hfill \\ & {\left. { - \frac{1}{{D\left( {1 - {\nu ^2}} \right)}}\left[ {\nu {{\bar M}_y}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right) - {{\bar M}_x}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)} \right]} \right\|^2} \hfill \end{split} $$ (23) $$ \begin{split} MS{E_{{\rm{P}}2}} = &\frac{1}{{{N_\varOmega }}}\sum\limits_{i = 1}^{{N_\varOmega }} {\left\| {\frac{{{\partial ^2}\bar w\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)}}{{\partial {y^2}}}} \right.} \hfill \\ & {\left.{- \frac{1}{{D\left( {1 - {\nu ^2}} \right)}}\left[ {\nu {{\bar M}_x}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right) - {{\bar M}_y}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right) } \right]} \right\|^2} \hfill \end{split} $$ (24) $$ MS{E_{{\rm{P}}3}} = \frac{1}{{{N_\varOmega }}}\sum\limits_{i = 1}^{{N_\varOmega }} {{{\left\| {\frac{{{\partial ^2}\bar w\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)}}{{\partial x\partial y}} + \frac{{{{\bar M}_{xy}}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)}}{{D(1 - \nu )}}} \right\|}^2}} $$ (25) 本文的误差函数表达式中包含挠度、弯矩对自变量的二阶偏导数项, 对于这些偏导项, 一方面可以利用神经网络的输出构造差分求解格式来近似求解, 但采用该方案需要较大的计算量才能得到精确的计算结果; 另一方面, 基于计算图的自动微分技术(automatic differentiation, AD)可以高效地处理神经网络对输入变量求导过程, 当前的深度学习框架如Tensoflow, Pytorch, MindsSpore等均支持自动微分. 本文基于Pytorch提供的自动微分接口实现对上述偏导项及误差函数梯度的计算.
在实际计算中, 也可灵活采用混合边界误差的形式进行求解, 如对于部分简单的边界条件构造特解, 而对于复杂的边界采用相应的无约束优化方案. 建立本文的误差函数后, 在每个训练批次(epoch)中均需计算其梯度并结合误差反向传播算法更新网络的内部参数, 关于该过程, Tang和Yang[21]对其进行了详细的讨论.
2.3 学习率选取方案
学习率的选取可直接影响神经网络的训练, 目前神经网络学习率的选取仍带有一定的经验性, 但总体而言, 在训练初期选取较大的学习率可以加快误差收敛速度, 在训练后期, 此时神经网络模型已经学习到相应的特征, 此时往往需要降低学习率, 以便对神经网络内部参数进行微调, 使得误差波动幅度不至于过大. 经过本文的实践, 本文学习率选取方案如下
$$ {\alpha _t} = \left\{ {\begin{array}{*{20}{l}} 1\times{{{10}^{ - 3}},}&{t \leqslant 5000} \\ {5 \times {{10}^{ - 4}},}&{5000 < t \leqslant 10\;000} \\ 1\times{{{10}^{ - 4}},}&{10\;000 < t \leqslant 20\;000} \\ {5 \times {{10}^{ - 5}},}&{20\;000 < t \leqslant 30\;000} \\ 1\times{{{10}^{ - 5}},}&{30\;000 < t \leqslant 40\;000} \\ {5 \times {{10}^{ - 6}},}&{40\;000 < t \leqslant 45\;000} \\ 1\times{{{10}^{ - 6}},}&{45\;000 < t} \end{array}} \right. $$ (26) 其中,
$t$ 为训练次数.2.4 算法流程
算法1. 本文算法
输入: 学习率
$\alpha $ ; 最大迭代次数$N$ ; 训练误差下限ϵ; 数据点生成参数${N_\varOmega }$ ,${N_{{\varGamma _1}}}$ ,${N_{{\varGamma _2}}}$ ,${N_{{\varGamma _3}}}$ 输出: 挠度网络; 弯矩网络
初始化神经网络模型及其内部参数
${{{{\boldsymbol{\theta}} }}_{\rm{w}}},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}$ , 训练误差$C$ 1,
$t$ =0 (initialization)2 while (
$C$ > = ϵ) and ($t$ <$N$ ) do{
(1)在求解域中随机生成数据点坐标
$\left( {{x_{\rm{d}}},{y_{\rm{d}}}} \right) \in \varOmega$ ,$\left( {{x_{\rm{b}}},{y_{\rm{b}}}} \right) \in$ $\partial \varOmega$ (2)将数据点坐标作为参数输入神经网络模型, 通过前向传播算法, 求得
$\bar w\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{w}}}} \right)$ ,${\bar M_x}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)$ ,${\bar M_{xy}}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)$ ,${\bar M_y}\left( {{\boldsymbol{x}};{{\boldsymbol{\theta }}_{\rm{m}}}} \right)$ (3)通过式(16)计算误差函数
$C$ (4)计算误差函数
$C$ 关于神经网络参数${{{{\boldsymbol{\theta }}}}_{\rm{w}}},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}$ 的梯度, 结合Adam优化算法更新${{{{\boldsymbol{\theta}} }}_{\rm{w}}},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}$ (5)
$ t \leftarrow t + 1 $ /}
算法中, 本文的数据点的默认生成参数为
${N_\varOmega }$ =500,${N_{{\varGamma _1}}}$ =${N_{{\varGamma _2}}}$ =${N_{{\varGamma _3}}}$ =100, 实际算例分析时会进行相应的调整; 在网络模型初始化中, 采用Xavier均匀分布[22]对参数${{{{\boldsymbol{\theta}} }}_{\rm{w}}},{{{{\boldsymbol{\theta}} }}_{\rm{m}}}$ 初始化.3. 算例
3.1 双向面内变刚度圆形薄板轴对称弯曲
图2所示受横向均布荷载
$q(x,y) = - {q_0}$ 作用的周边固支圆形薄板, 其半径R,$\nu = 0.3$ , 弯曲刚度函数沿半径变化$D(\rho ) = {D_0}{{\text{e}}^{ - m\rho /a}}$ , 其中$\rho = \sqrt {{x^2} + {y^2}} $ ,$m$ 为梯度系数,${q_0}$ ,${D_0}$ 为常数. 该问题存在理论解[14].本文选取梯度参数m分别为0, 0.5, 1, 2的情况进行计算. 本算例在计算过程中仅需要施加位移边界条件, 设
$w\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{w}}}} \right)$ 为挠度网络的输出, 考虑到本算例边界条件较为简单, 构造挠度试函数为${w^*} = w\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{w}}}} \right) $ $ {\left( {{x^2} + {y^2} - {R^2}} \right)^2}$ .本算例的挠度模型以及弯矩模型均采用具有6层隐藏层, 每层隐藏层具备30个神经元的网络结构(记为6 × 30), 采用
${x^2}$ 作为激活函数,${k_{\rm{p}}} = 100$ ,${N_\varOmega } = 640$ ; 采用Adam优化算法, 对各个工况所采用的学习率方案均一致. 本文算例在不同的梯度下上述基本参数不变, 仅更改梯度系数. 为更加详细地显示误差函数的变化情况, 如无特殊说明, 本文均取误差的十进制对数作为等效误差并绘制相应的训练误差曲线.每个梯度参数下神经网络模型的计算误差如图3所示, 由此可见, m = 0时为刚度恒定的薄板, 此时误差函数收敛较快, 相比其他工况其误差最终的收敛值最小; 随着梯度系数的增大, 训练误差最终的收敛值出现增大的趋势.
采用无量纲计算公式
$\bar w = \dfrac{{{{10}^2}{D_0}}}{{{q_0}{R^4}}}{w^*}$ ,$\bar M = \dfrac{{10{M^*}{D_0}}}{{{q_0}{R^2}}}$ 将计算所得结果无量纲化, 其中${w^*}$ 和${M^*}$ 为实际计算所得的挠度、弯矩. 由图3和图4可以看出本文方法计算的结果与理论解答吻合. 由表1的挠度计算结果可见, 本文方法在挠度的求解上达到了相当高的精度, 本算例也选取各个梯度参数下薄板部分点的弯矩${\bar M_x}$ 进行对比分析, 结果如表2所示, 理论解答中弯矩在原点处不存在解答, 故选取离原点较近的点(0.02, 0)求得相应的理论解, 由于本文采用了两个网络模型进行计算, 本文解中,${\bar M_{x1}}$ 为弯矩网络直接输出的结果,${\bar M_{x2}}$ 为根据挠度网络由式(4)求得, 两者相对理论解的误差分别为相对误差1、相对误差2. 可以发现当梯度系数m = 0时, 此时薄板为刚度恒定的薄板,${\bar M_{x1}}$ 和${\bar M_{x2}}$ 均取得较高的精度, 而随着梯度系数的增大, 弯矩网络在靠近原点处所求得的弯矩与理论解的误差也随之增大, 在点(0.02, 0)处, 工况m = 2时${\bar M_{x1}}$ 误差最大, 为4.086%, 而挠度网络计算的结果并未出现较大的差别, 这说明梯度系数的增大会导致本文模型的挠度与弯矩耦合误差在靠近原点处增大, 这也解释了图2所示不同梯度下神经网络误差收敛程度不同的原因, 但这并不影响挠度模型的求解精度.![]()
图 4 PINN求解圆形面内变刚度薄板弯曲问题的训练误差收敛曲线(m = 2)Figure 4. Training error convergence curve of PINN (m = 2)表 1 本文方法计算$\bar w$ 与理论解对比(无量纲)Table 1. Comparison of dimensionless$\bar w$ calculated by neural network method and the theoretical solutionm (x, y) $ \bar {w} $ Theory Relative error/% 0 (0.0, 0) −1.5625 −1.5625 0 (0.4, 0) −1.1024 −1.1025 −0.009 (0.8, 0) −0.2025 −0.2025 0 1 (0.0, 0) −0.8420 −0.8418 0.024 (0.4, 0) −0.5517 −0.5525 −0.145 (0.8, 0) −0.0877 −0.0877 0 2 (0.0, 0) −0.4520 −0.4522 −0.044 (0.4, 0) −0.2725 −0.2725 0 (0.8, 0) −0.0374 −0.0373 0.267 注: 本文误差计算公式为$e = \dfrac{{u - {u^*}}}{u} \times 100\% $, $u$为本文方法计算结果, ${u^*}$为理论解 (Note: The relative error in this paper is calculated by $e = \dfrac{ {u - {u^*} } }{u} \times $$ 100\%$, where $u$ is the calculation results of this paper, ${u^*}$ is the theoretical solution) 表 2 本文方法计算${\bar M_x}$ 与理论解对比(无量纲)Table 2. Comparison of dimensionless${\bar M_x}$ calculated by neural network method and the theoretical solutionm (x, y) ${\bar M_{x1}}$ ${\bar M_{x2}}$ Theory Relative error 1/% Relative error 2/% 0 (0.02, 0) −0.8122 −0.8118 −0.8117 0.062 0.020 (0.4, 0) −0.4828 −0.4824 −0.4825 0.071 −0.027 (0.8, 0) 0.5077 0.5075 0.5075 0.037 −0.004 (1.0, 0) 1.2512 1.2499 1.2500 0.094 −0.008 1 (0.02, 0) −0.6027 −0.5910 −0.5853 2.885 0.956 (0.4, 0) −0.3064 −0.3076 −0.3103 −1.278 −0.886 (0.8, 0) 0.6264 0.6249 0.6282 −0.287 −0.524 (1.0, 0) 1.3536 1.3634 1.3525 0.077 0.797 2 (0.02, 0) −0.4335 −0.4150 −0.4158 4.086 −0.192 (0.4, 0) −0.1652 −0.1651 −0.1695 −2.599 −2.684 (0.8, 0) 0.7273 0.7268 0.7319 −0.637 −0.712 (1.0, 0) 1.4344 1.4543 1.4409 −0.458 0.921 为了说明本文方法在求解面内变刚度薄板弯曲问题上的优点, 本算例也利用PINN来求解其四阶偏微分控制方程(12), 采用隐藏层层数为6, 每层隐藏层具有30个神经元的神经网络模型对工况m = 2进行求解, 激活函数为Tanh函数, 训练的数据点由求解域中随机生成, 数据点的产生有两种方案:
方案(1)为在整个求解域中随机生成训练数据点, 此时的数据点可在原点附近生成;
方案(2)为在求解域中随机生成的数据点但离原点较远. 此时两方案在训练过程中的误差收敛情况如图4所示, 可见采用相同的模型, 而生成的数据点不同则会导致模型的训练出现不同的结果, 虽然采用数据点生成方案(1)的模型训练也收敛, 但由于其误差此时收敛于一个较大的值, 得不到正确解. 经过本文分析, 这主要是由于本算例的弯曲刚度函数D的二阶导数在靠近原点区域出现“爆炸”式变化的原因, 即刚度函数的二阶偏导在原点处不收敛, 在靠近原点处
$\dfrac{{{\partial ^2}D}}{{\partial {x^2}}}$ 等的解答急剧增大, 这会导致PINN采用方案(1)训练时, 遇到靠近原点处的点, 计算所得域内误差突然增大, 进而导致误差训练难以收敛. PINN最初提出时并未考虑求解域内存在奇异点的情况, 对于该情况, 一般情况下可在生成的数据点中排除掉奇异点, 但对于本算例中奇异点处被施予荷载的情况, 如果不能很好地处理则会影响求解的精度. 对此, 本文认为可以弱化相应的偏微分控制方程再利用PINN求解, 也可参考本文思路, 结合神经网络解法的特点, 根据具体问题对原偏微分控制方程等效化处理. 本文方法在求解时并非直接从方程(12)入手, 而是通过求解一系列偏微分方程组来逼近真实解, 避开了对弯曲刚度函数求偏导数, 故其求解仅与域内的刚度值有关, 其适应性更强, 对薄板弯曲问题的求解更具“鲁棒性”.3.2 单向面内变刚度方形薄板受非线性荷载作用
如图5所示边长为
$a$ 的方形薄板, 厚度$h$ ,$\nu = 0.3$ , 1, 2, 3边固支, 4边简支, 其弯曲刚度函数为$D(x) = {D_0}{(x + 1)^m}$ ,${D_0} = \dfrac{{{E_0}h}}{{12\left( {1 - {\nu ^2}} \right)}}$ ,${E_0}$ 为常数, m为梯度系数, 受横向非线性荷载$q(x) = - {{\text{e}}^{\frac{x}{a}}}$ 作用, 利用本文方法求解m分别为0, 1, 2情况下的挠度、内力.本算例选取6 × 30的挠度网络模型, 5 × 50的弯矩网络模型, 激活函数选择Tanh函数,
${k_{\rm{p}}} = $ 100,${N_\varOmega } = 450$ , 采用Adam优化算法. 训练误差曲线如图6所示. 设$w\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{w}}}} \right)$ 为挠度网络的输出, 弯矩网络的输出为$M\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right) = \left\{ {{M_x},{M_{xy}},{M_y}} \right\}$ , 考虑本算例的位移及广义应力边界条件, 挠度及弯矩试函数构造为$$ \left. \begin{array}{l} {w^*} = w\left( {x,y;{{\boldsymbol{\theta }}_{\rm{w}}}} \right){\left[\left(\sqrt 3 ax + y - a\right)\left( - \sqrt 3 ax + y - a\right)\right]^2} \hfill \\ \left\{ {{M_x}^*,{M_{xy}}^*,{M_y}^*} \right\} = \left\{ {{M_x}(x - 1),{M_{xy}},{M_y}} \right\} \hfill \end{array} \right\} $$ 如果将本算例的应力边界条件也通过挠度网络进行构造, 则挠度表达式将复杂许多.
将本算例的计算结果与有限元解答对比, 有限元计算中每个单元的弯曲刚度根据单元的形心坐标计算, 采用50
$ \times $ 50的矩形薄板非协调单元来对求解域进行离散, 离散方案通过小片测试, 对本文的3种工况该网格离散方案均收敛. 将本文与有限元计算结果的无量纲挠度、弯矩进行对比分析, 无量纲计算公式为$ \overline{W}=\dfrac{{10}^{2}{D}_{0}}{{q}_{0}{a}^{4}} {w}^{*} $ ,$\bar M = \dfrac{{100{M^*}{D_0}}}{{{q_0}{a^2}}}$ ,${w^*}$ 和${M^*}$ 为实际计算所得的挠度、弯矩值.由挠度计算图7可知, 在挠度的求解上, 本文解法与有限元解法一致. 由图8的弯矩对比图可发现, 本文弯矩解与有限元解答基本吻合, 而当梯度系数m = 2时, 虽然本文解与有限元解在部分点上弯矩的相对误差增大, 但整体上解答吻合.
![]()
图 7 不同梯度参数下本文挠度计算结果与有限元对比(y = 0.5 m)Figure 7. Comparison of dimensionless deflection calculation results in this paper with FEM when different gradient parameters (y = 0.5 m)![]()
图 8 不同梯度参数下本文弯矩$ {\overline{M}}_{x1} $ 计算结果与有限元对比(y = 0.5 m)Figure 8. Comparison of dimensionless bending moment$ {\overline{M}}_{x1} $ calculation results of this paper with FEM when different gradient parameters (y = 0.5 m)3.3 三角形面内变刚度薄板受线性分布荷载作用
工程中地下掩体结构有时可看作三角形薄板结构. 如图9所示, 厚度为
$h$ 的等边三角形面内变刚度薄板两斜边固支, 水平边作自由边考虑, 弯曲刚度函数为$D(y) = {D_0}\left[ {1 + m\sin \left( {\dfrac{{\text{π}} }{a}y} \right)} \right]$ , 其中${D_0} = \dfrac{{{E_0}h}}{{12\left( {1 - {\nu ^2}} \right)}}$ ,${E_0}$ 为常数,$\nu = $ 0.3. 其受线性分布横向荷载$q(y) = $ $ - \dfrac{{{q_0}y}}{a}$ , 计算m = 0, 0.2, 0.5时薄板的挠度、弯矩.本算例的挠度网络与弯矩网络分别采用隐藏层层数为6和5, 每层分别具备30和50个神经元的神经网络结构进行计算, 两个模型均采用
${x^2}$ 作为激活函数, 采用Adam优化算法,${k_{\rm{p}}} = 100$ ,${N_\varOmega } = 560$ ,${N_{{\varGamma _3}}} = 50$ .本算例的挠度、弯矩试函数为
$$ {\begin{array}{*{20}{l}} {{w^*} = w\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{w}}}} \right){{\left[\left(\sqrt 3 ax + y - a\right)\left( - \sqrt 3 ax + y - a\right)\right]}^2}} \\ \left\{ {{M_x}^*,{M_{xy}}^*,{M_y}^*} \right\} = \hfill \\ \qquad \left\{ {{M_x}\left( {x,y;{{{{\boldsymbol{\theta}} }}_{\rm{m}}}} \right),{M_{xy}}\left( {x,y;{{\boldsymbol{\theta }}_{\rm{m}}}} \right)y,{M_y}\left( {x,y;{{\boldsymbol{\theta }}_{\rm{m}}}} \right)y} \right\} \hfill \end{array}} $$ 此时自由边的应力边界条件并未精确满足, 故在误差函数中引入自由边的剪力误差项
$$ MS{E_{{\varGamma _3}}} = \frac{1}{{{N_{{\varGamma _3}}}}}\sum\limits_{i = 1}^{{N_{{\varGamma _3}}}} {{{\left\| {\frac{{\partial {M_y}\left( {x,y;{{\boldsymbol{\theta}} _{\rm{m}}}} \right)}}{{\partial y}} + \frac{{\partial {M_{xy}}\left( {x,y;{{\boldsymbol{\theta}} _{\rm{m}}}} \right)}}{{\partial x}}} \right\|}^2}} $$ 根据式(16), 误差函数为
$$ C\left( {{{\boldsymbol{\theta }}_w},{{\boldsymbol{\theta }}_{{m}}}} \right) = MS{E_\varOmega } + {k_{\rm{p}}}MS{E_{\rm{p}}} + {k_3}MS{E_{{\varGamma _3}}} $$ 神经网络训练收敛后, 采用无量纲计算公式
$\bar w = \dfrac{{{{10}^3}{D_0}}}{{{q_0}{a^4}}}w$ ,$\bar M = \dfrac{{{{10}^2}M{D_0}}}{{{q_0}{a^2}}}$ , 将挠度网络及弯矩网络的输出无量纲化后与有限元计算结果对比, 有限元采用DKT薄板单元进行求解, 有限元计算收敛时单元总数为3032, 对比结果如图10 ~ 图12所示. 可以看出, 本文方法求解结果与有限元计算结果一致, 同时由于对挠度、弯矩输出进行了构造, 使得大部分边界条件严格满足, 也可看出弯矩在边界与边界交接处的收敛情况良好, 在一定程度上, 本文方法可避免由于边界处误差难以收敛而带来的影响.![]()
图 10 三角形面内变刚度薄板沿轴线x = 0上的挠度$ \bar {w} $ 分布Figure 10. Dimensionless deflection variation of thin triangular plate with in-plane stiffness gradient along axis x = 0![]()
图 11 m = 0.2时三角形面内变刚度薄板弯矩的有限元计算结果Figure 11. Finite element calculation of bending moment of thin plate with triangular in-plane variable stiffness (m = 0.2)![]()
图 12 m = 0.2时三角形面内变刚度薄板弯矩的本文计算结果Figure 12. Neural network method calculation of bending moment of thin plate with triangular in-plane variable stiffness (m = 0.2)4. 本文方法计算效率及其有限元对比分析
以
$T$ 表示算例各个工况下的平均用时,${T^\prime }$ 表示各个工况下误差函数的构建及其梯度计算的平均耗时, 本文各算例求解时迭代所需的数据点数及所需的平均时间、内存如表3所示, 可看出本文误差函数的构建及其梯度的求解占据了神经网络训练总时长的70%左右.本文的有限元求解程序采用python语言进行编写, 选用每个节点有3个自由度的薄板弯曲单元计算, 刚度矩阵以 compressed sparse column (CSC)格式的稀疏矩阵存储, 利用科学计算库scipy中的线性求解器求解刚度方程. 根据有限元解答的最小挠度判断有限元解答是否收敛, 算例2和算例3中各个工况下有限元计算收敛时的节点数目如表4所示, 本文有限元计算所需的节点数与计算所需内存的关系如图13所示.
表 3 本文各算例求解所需的数据点数、内存、时间Table 3. The number of training data points, computational memory and computing time of numerical examples in this paperNumerical example Number of data points Computational memory/MiB $T(s)$ ${T^\prime }(s)$ ${T^\prime }(s)/T(s)$/% 1 560 121.3 551.99 374.66 67.8 2 450 139.4 673.81 475.22 70.5 3 610 169.3 646.14 479.81 74.3 表 4 算例2、算例3的有限元求解收敛所需节点数目Table 4. The number of nodes needed for the convergence of the finite element solution of numerical example 2 and 3Numerical example m Number of nodes 2 0 441 1 841 2 1681 3 0 305 0.2 1331 0.5 1604 ![]()
图 13 有限元计算所需内存与节点数的关系(薄板单元)Figure 13. The relationship between the number of nodes and the memory needed in finite element calculation using thin plate bending element可以发现有限元解答收敛时所需的节点(网格)数目随着梯度系数的增大而增大, 节点数与所需内存并非呈现线性关系, 节点数的增多会导致所需计算内存的急剧增大, 而神经网络方法在梯度系数增大时, 在单次迭代时仍可以较少的数据点进行迭代, 这使得本文方法在求解时所需内存较小. 在时间上, 以算例3中梯度系数m = 0.5为例, 此时有限元求解收敛时, 所需内存为291.5 MiB, 刚度方程从组装到求解用时1.798 s, 可以看出本文解法的求解速度明显慢于有限元, 一方面是由于本文方法求解的是强形式的偏微分控制方程, 与求解弱形式的方程相比, 往往需要更多的迭代次数, 另一方面, 神经网络方法是一类以数据为驱动的解法, 在求解过程中需要往复迭代, 而对于本文研究的线性问题, 有限元仅需求解一次刚度方程.
5. 网络耦合系数
$ {k}_{{\rm{p}}} $ 、隐藏层层数与神经元个数、激活函数对网络收敛性的影响分析5.1 网络耦合系数
${k_{\rm{p}}}$ 本文模型由于引入两个神经网络模型对面内变刚度薄板弯曲问题进行计算, 与采用单网络的神经网络方法相比, 在训练过程中需考虑网络之间的耦合误差, 本文引入网络耦合系数
$ {k}_{{\rm{p}}} $ 以加强网络之间的耦合,${k_{\rm{p}}}$ 的取值过小会导致网络耦合不佳, 过大则可能会导致训练不收敛, 因此需要对${k}_{{\rm{p}}}$ 的取值进行讨论, 选取${k_{\rm{p}}}$ =1, 10, 100, 1000四种情况对算例2进行计算分析(梯度系数m = 2), 以网络训练时薄板的变形能、外力功变化情况衡量不同${k_{\rm{p}}}$ 下本文模型的训练效果.薄板变形能计算公式为
$$ {W_{{\rm{int}}}} = \frac{1}{2}\int_\varOmega {{{\boldsymbol{k}}^{\rm{T}}}} {\boldsymbol{M}}{\text{d}}\varOmega = \frac{1}{2}\int_\varOmega {{{\boldsymbol{k}}^{\rm{T}}}} {\boldsymbol{Dk}}{\text{d}}\varOmega $$ (27) 其中
$$ {\boldsymbol{D}} = D(x,y)\left[ {\begin{array}{*{20}{c}} 1&v&0 \\ v&1&0 \\ 0&0&{(1 - v)/2} \end{array}} \right] $$ (28) $$ {\boldsymbol{k}} = \left\{ {\begin{array}{*{20}{c}} {{k_x}} \\ {{k_y}} \\ {{k_{xy}}} \end{array}} \right\} = - \left\{ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^2}w}}{{\partial {x^2}}}} \\ {\dfrac{{{\partial ^2}w}}{{\partial {y^2}}}} \\ {2\dfrac{{{\partial ^2}w}}{{\partial x\partial y}}} \end{array}} \right\} $$ (29) $ {\boldsymbol{k}} $ 为根据弯矩网络的输出由式(8)和式(9)计算所得.外力所作实功为
$$ {W}_{\text{ext }}=\frac{1}{2}{\displaystyle {\int }_{\varOmega }q}w \text{d}\varOmega $$ (30) 其则根据挠度网络的输出值计算.
采用高斯积分, 由上述公式计算神经网络训练时薄板的变形能、应变能的变化情况, 根据挠度网络输出计算外力实功
${W_{{\rm{ext}}}}$ , 根据弯矩网络输出计算薄板变形能${W_{{\rm{int}}}}$ . 当神经网络解答收敛时, 根据能量原理, 其外力实功应与变形能相等. 此时各工况下的变形能与外力实功在神经网络训练过程中的变化情况如图14所示. 可以看出, 随着训练的进行, 外力实功与薄板的变形能逐渐收敛, 在各个工况下, 变形能的收敛速度相比于外力功要慢, 进而说明挠度网络收敛速度较快. 当${k_{\rm{p}}}$ =1时, 此时薄板变形能与外力功收敛时, 两者相差较大; 而随着${k}_{{\rm{p}}}$ 的增大, 神经网络计算的变形能收敛加快, 此时当变形能与外力功收敛时, 两者之间的差别减小, 而当${k_{\rm{p}}}$ =1000时, 其对训练的影响与${k_{\rm{p}}}$ =100时差别不大, 在实际计算中, 该系数不能过大, 过大会导致神经网络学习不到正确的特征, 进而导致训练不收敛.![]()
图 14 神经网络训练时不同${k_{\rm{p}}}$ 下薄板变形能与外力功的变化情况Figure 14. Changes in the deformation energy and the external force work during the training process with different${k_{\rm{p}}}$ 5.2 隐藏层层数与神经元个数
隐藏层层数与每层的神经元个数的选取均会影响到神经网络的训练效率, 不失一般性, 在其他计算参数不变的情况下, 本节选取算例2中m = 2的情况分别讨论隐藏层层数、神经元个数的改变对计算过程中误差函数收敛情况的影响, 计算结果如图15所示.
![]()
图 15 隐藏层数、每层神经元个数对计算误差的收敛影响对比Figure 15. Comparison of the effects of different number of hidden layers and neurons on the convergence of loss function可以发现, 在一定程度内误差函数的收敛速度随着网络层数、隐藏层每层神经元数的增大而加快, 当两者都增大到一定程度时, 此时误差的收敛速度会趋向于一个“饱和”状态. 适当地增加层数或神经元的个数有利于误差函数的收敛, 目前在利用神经网络方法求解偏微分方程的研究中, 在两者的选取上仍带有一定的经验性. 结合上述的计算结果, 考虑到计算机硬件能力的限制, 本文算例的隐藏层层数在4 ~ 6层间选取, 每层神经元数目在30 ~ 50之间选取.
5.3 激活函数
在神经网络的训练过程中, 有多种非线性激活函数可以选择, 非线性的激活函数是使神经网络具备拟合非线性函数能力的重要原因, 常用的激活函数有Tanh, ReLU, Sigmoid, Swish函数等. 为讨论激活函数对神经网络训练的影响, 其余计算参数不变, 本节选取
${x^2}$ 作为激活函数并与Tanh, Swish函数进行对比分析, 对算例1、算例3进行计算, 训练过程中的误差走向如图16所示. 在误差函数的收敛速度上,${x^2}$ 优于Tanh, Swish函数, 同时由于其函数形式较为简单, 在自动微分计算中其所需计算量较小. 本文经验表明, 在薄板弯曲问题的求解上, 采用多项式函数${x^2}$ 作为激活函数可加快神经网络的收敛.![]()
图 16 不同的激活函数对计算误差的收敛影响对比Figure 16. Comparison of the effects of different activation functions on the convergence of loss function6. 结论与展望
本文基于深度学习技术与强形式的求解方案建立了一种直角坐标下求解面内变刚度薄板弯曲问题的神经网络方法, 通过几个算例分析, 得出以下结论:
(1)本文解答与理论解、有限元解吻合, 证明了本文方法在求解面内变刚度薄板弯曲问题上的正确性, 本文的神经网络模型不需要对弯曲刚度函数求偏导, 其适应性更强, 同时在薄板的位移边界条件、应力边界条件的施加上较为方便.
(2)本文方法属于强形式的数值解法, 其计算所得结果具备连续性与可导性. 理论上, 本文方法可以求解弹性模量以及厚度在面内连续变化的薄板弯曲问题. 弯矩网络的求解受到梯度系数的影响, 在梯度变化较大处弯矩网络的求解精度受到一定的影响, 但对挠度网络的求解精度影响不大.
(3)由于神经网络方法为迭代类解法, 本文方法在薄板线性弯曲问题求解上的收敛速度较有限元慢, 但其计算所需内存较小. 通过本文的模型结构可看出, 神经网络方法具备相当大的灵活性, 根据这一特点, 可进一步发展求解面内变刚度功能梯度薄板非线性弯曲问题的神经网络方法, 神经网络方法在非线性问题的求解中具备潜在优势.
(4)本文模型仍存在优化空间: 一方面在本文模型的训练过程中, 误差函数及其梯度的计算在整个训练过程中占据大部分的时间, 可以考虑优化误差函数的构建过程, 如引入有限元中形函数的思想对算法进行优化; 另一方面为了使得挠度、弯矩网络具备较强的表达能力, 本文模型采用了两个具有独立参数的网络进行计算, 这导致了本文模型的训练参数较多, 为此后续优化中可将本文的两个网络合并为一个网络(2个输入, 4个输出), 对网络结构进行改良, 以减少训练参数.
-
![]()
图 4 PINN求解圆形面内变刚度薄板弯曲问题的训练误差收敛曲线(m = 2)
Figure 4. Training error convergence curve of PINN (m = 2)
![]()
图 7 不同梯度参数下本文挠度计算结果与有限元对比(y = 0.5 m)
Figure 7. Comparison of dimensionless deflection calculation results in this paper with FEM when different gradient parameters (y = 0.5 m)
![]()
图 8 不同梯度参数下本文弯矩
$ {\overline{M}}_{x1} $ 计算结果与有限元对比(y = 0.5 m)Figure 8. Comparison of dimensionless bending moment
$ {\overline{M}}_{x1} $ calculation results of this paper with FEM when different gradient parameters (y = 0.5 m)![]()
图 10 三角形面内变刚度薄板沿轴线x = 0上的挠度
$ \bar {w} $ 分布Figure 10. Dimensionless deflection variation of thin triangular plate with in-plane stiffness gradient along axis x = 0
![]()
图 11 m = 0.2时三角形面内变刚度薄板弯矩的有限元计算结果
Figure 11. Finite element calculation of bending moment of thin plate with triangular in-plane variable stiffness (m = 0.2)
![]()
图 12 m = 0.2时三角形面内变刚度薄板弯矩的本文计算结果
Figure 12. Neural network method calculation of bending moment of thin plate with triangular in-plane variable stiffness (m = 0.2)
![]()
图 13 有限元计算所需内存与节点数的关系(薄板单元)
Figure 13. The relationship between the number of nodes and the memory needed in finite element calculation using thin plate bending element
![]()
图 14 神经网络训练时不同
${k_{\rm{p}}}$ 下薄板变形能与外力功的变化情况Figure 14. Changes in the deformation energy and the external force work during the training process with different
${k_{\rm{p}}}$ ![]()
图 15 隐藏层数、每层神经元个数对计算误差的收敛影响对比
Figure 15. Comparison of the effects of different number of hidden layers and neurons on the convergence of loss function
![]()
图 16 不同的激活函数对计算误差的收敛影响对比
Figure 16. Comparison of the effects of different activation functions on the convergence of loss function
表 1 本文方法计算
$\bar w$ 与理论解对比(无量纲)Table 1 Comparison of dimensionless
$\bar w$ calculated by neural network method and the theoretical solutionm (x, y) $ \bar {w} $ Theory Relative error/% 0 (0.0, 0) −1.5625 −1.5625 0 (0.4, 0) −1.1024 −1.1025 −0.009 (0.8, 0) −0.2025 −0.2025 0 1 (0.0, 0) −0.8420 −0.8418 0.024 (0.4, 0) −0.5517 −0.5525 −0.145 (0.8, 0) −0.0877 −0.0877 0 2 (0.0, 0) −0.4520 −0.4522 −0.044 (0.4, 0) −0.2725 −0.2725 0 (0.8, 0) −0.0374 −0.0373 0.267 注: 本文误差计算公式为$e = \dfrac{{u - {u^*}}}{u} \times 100\% $, $u$为本文方法计算结果, ${u^*}$为理论解 (Note: The relative error in this paper is calculated by $e = \dfrac{ {u - {u^*} } }{u} \times $$ 100\%$, where $u$ is the calculation results of this paper, ${u^*}$ is the theoretical solution) 
下载: 导出CSV
表 2 本文方法计算
${\bar M_x}$ 与理论解对比(无量纲)Table 2 Comparison of dimensionless
${\bar M_x}$ calculated by neural network method and the theoretical solutionm (x, y) ${\bar M_{x1}}$ ${\bar M_{x2}}$ Theory Relative error 1/% Relative error 2/% 0 (0.02, 0) −0.8122 −0.8118 −0.8117 0.062 0.020 (0.4, 0) −0.4828 −0.4824 −0.4825 0.071 −0.027 (0.8, 0) 0.5077 0.5075 0.5075 0.037 −0.004 (1.0, 0) 1.2512 1.2499 1.2500 0.094 −0.008 1 (0.02, 0) −0.6027 −0.5910 −0.5853 2.885 0.956 (0.4, 0) −0.3064 −0.3076 −0.3103 −1.278 −0.886 (0.8, 0) 0.6264 0.6249 0.6282 −0.287 −0.524 (1.0, 0) 1.3536 1.3634 1.3525 0.077 0.797 2 (0.02, 0) −0.4335 −0.4150 −0.4158 4.086 −0.192 (0.4, 0) −0.1652 −0.1651 −0.1695 −2.599 −2.684 (0.8, 0) 0.7273 0.7268 0.7319 −0.637 −0.712 (1.0, 0) 1.4344 1.4543 1.4409 −0.458 0.921
下载: 导出CSV
表 3 本文各算例求解所需的数据点数、内存、时间
Table 3 The number of training data points, computational memory and computing time of numerical examples in this paper
Numerical example Number of data points Computational memory/MiB $T(s)$ ${T^\prime }(s)$ ${T^\prime }(s)/T(s)$/% 1 560 121.3 551.99 374.66 67.8 2 450 139.4 673.81 475.22 70.5 3 610 169.3 646.14 479.81 74.3
下载: 导出CSV
表 4 算例2、算例3的有限元求解收敛所需节点数目
Table 4 The number of nodes needed for the convergence of the finite element solution of numerical example 2 and 3
Numerical example m Number of nodes 2 0 441 1 841 2 1681 3 0 305 0.2 1331 0.5 1604
下载: 导出CSV
-
[1] 徐芝纶. 弹性力学简明教程. 北京: 高等教育出版社, 2013 (Xu Zhilun. A Concise Course in Elasticity. Beijing: Higher Education Press, 2013 (in Chinese))
[2] 于天崇, 聂国隽, 仲政. 变刚度矩形板弯曲问题的Levy解. 力学季刊, 2012, 33(1): 53-59 (Yu Tianchong, Nie Guojun, Zhong Zheng. Levy-type solution for the bending of rectangular plates with variable stiffness. Chinese Quarterly of Mechanics, 2012, 33(1): 53-59 (in Chinese) doi: 10.3969/j.issn.0254-0053.2012.01.008 [3] 朱竑祯, 王纬波, 高存法等. 面内变刚度功能梯度圆形薄板的轴对称弯曲. 船舶力学, 2018, 22(11): 1364-1375 (Zhu Hongzhen, Wang Weibo, Gao Cunfa, et al. Axisymmetric bending of functionally graded thin circular plates with in-plane stiffness gradient. Journal of Ship Mechanics, 2018, 22(11): 1364-1375 (in Chinese) doi: 10.3969/j.issn.1007-7294.2018.11.006 [4] 何建璋, 褚福运, 仲政等. 面内变刚度矩形薄板自由振动问题的辛弹性分析. 同济大学学报(自然科学版), 2013, 41(9): 1310-1317 (He Jianzhang, Chu Fuyun, Zhong Zheng, et al. Symplectic elasticity approach for free vibration of a rectangular plate with in-plane material inhomogeneity. Journal of Tongji University (Natural Science) , 2013, 41(9): 1310-1317 (in Chinese) doi: 10.3969/j.issn.0253-374x.2013.09.006 [5] Santare MH, Lambros J. Use of graded finite elements to model the behavior of nonhomogeneous materials. Journal of Applied Mechanics, 2000, 67(4): 819-822 doi: 10.1115/1.1328089
[6] Kim JH, Paulino GH. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials. Journal of Applied Mechanics, 2002, 69(4): 502-514
[7] 黄立新, 姚祺, 张晓磊等. 基于分层法的功能梯度材料有限元分析. 玻璃钢/复合材料, 2013, 2: 43-48 (Huang Lixin, Yao Qi, Zhang Xiaolei, et al. Finite element analysis of functionally graded materials based on layering method. Composites Science and Engineering, 2013, 2: 43-48 (in Chinese) [8] 田云德, 秦世伦. 梯度厚板热应力分层计算方法的改进. 计算力学学报, 2009, 26(6): 961-965 (Tian Yunde, Qin Shilun. The improvement of stratified method about thermal stress of functionally gradient material thick plates. Chinese Journal of Computational Mechanics, 2009, 26(6): 961-965 (in Chinese) [9] Lee H, Kang IS. Neural algorithm for solving differential equations. Journal of Computational Physics, 1990, 91(1): 110-131 doi: 10.1016/0021-9991(90)90007-N
[10] Lagaris IE, Likas A, Fotiadis DI. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 1998, 9(5): 987-1000
[11] Eukinan W, Yu B. The deep ritz method: a deep learning-based numerical algorithm for solving variational problems. Communications in Mathematics and Statistics, 2018, 6(1): 1-12 doi: 10.1007/s40304-018-0127-z
[12] Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 2019, 378: 686-707 doi: 10.1016/j.jcp.2018.10.045
[13] Sirignano JA, Spiliopoulos K. DGM: A deep learning algorithm for solving partial differential equations. Journal of Computational Physics, 2018, 375: 1339-1364 doi: 10.1016/j.jcp.2018.08.029
[14] Samaniego E, Anitescu C, Goswami S, et al. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020: 362
[15] 瞿同明, 冯云田, 王孟琦等. 基于深度学习和细观力学的颗粒材料本构关系研究. 力学学报, 2021, 53(7): 1-12 (Qu Tongming, Feng Yuntian, Wang Mengqi, et al. Constitutive relations of granular materials by integrating micromechanical knowledge with deep learning. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1-12 (in Chinese) [16] 谢晨月, 袁泽龙, 王建春等. 基于人工神经网络的湍流大涡模拟方法. 力学学报, 2021, 53(1): 1-16 (Xie Chenyue, Yuan Zelong, Wang Jianchun, et al. Artificial neural network-based subgrid-scale models for large-eddy simulation of turbulence. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(1): 1-16 (in Chinese) [17] 刘宇翔, 王东东, 樊礼恒等. 基于卷积神经网络的无网格形函数影响域优化研究. 固体力学学报, 2021, 42(3): 302-319 (Liu Yuxiang, Wang Dongdong, Fan Liheng, et al. A CNN-based approach for optimizing support selection of meshfree methods. Chinese Journal of Solid Mechanics, 2021, 42(3): 302-319 (in Chinese) [18] 郭宏伟, 庄晓莹. 采用两步优化器的深度配点法与深度能量法求解薄板弯曲问题. 固体力学学报, 2021, 42(3): 249-266 (Guo Hongwei, Zhuang Xiaoying. The application of deep collocation method and deep energy method with a two-step optimizer in the bending analysis of kirchhoff thin plate. Chinese Journal of Solid Mechanics, 2021, 42(3): 249-266 (in Chinese) [19] 陈豪龙, 柳占立. 基于数据驱动模型求解热传导反问题. 计算力学学报, 2021, 38(3): 272-279 (Chen Haolong, Liu Zhanli. Solving the inverse heat conduction problem based on data driven model. Chinese Journal of Computational Mechanics, 2021, 38(3): 272-279 (in Chinese) [20] Cai Z, Chen J, Liu M, et al. Deep least-squares methods:an unsupervised learning-based numerical method for solving elliptic PDEs. Journal of Computational Physics, 2020: 420
[21] Tang SQ, Yang Y. Why neural networks apply to scientific computing? Theoretical and Applied Mechanics Letters, 2021 (in press)
[22] Glorot X, Bengio Y. Understanding the difficulty of training deep feedforward neural networks. Journal of Machine Learning Research, 2010, 9: 249-256
-
期刊类型引用(6)
1. 郭远，傅卓佳，闵建，刘肖廷，赵海涛. 课程-迁移学习物理信息神经网络用于长时间非线性波传播模拟. 力学学报. 2024(03): 763-773 . 
本站查看
2. 郭万金，吴广，朱恒，陈柱宏，詹子立，曹雏清，赵立军，郭磊，朱文锋. 考虑廓形特征的机器人异型盘磨削接触力建模. 表面技术. 2024(12): 167-180 . 百度学术
3. 方卫华，徐孟启，王润英. 基于物理信息神经网络的薄板反问题研究. 固体力学学报. 2023(04): 483-496 . 百度学术
4. 陈俊，邓立克，袁凯，张灿辉，王东东. 无网格法驱动径向基函数前置卷积神经网络. 力学季刊. 2023(03): 512-524 . 百度学术
5. 欧阳晔，江巍，吴怡，冯强，郑宏. 广义乘子法求解构造变分问题的神经网络方法. 工程力学. 2023(11): 11-20 . 百度学术
6. 陈健 ，王东东 ，刘宇翔 ，陈俊 . 无网格动力分析的循环卷积神经网络代理模型. 力学学报. 2022(03): 732-745 . 本站查看
其他类型引用(5)
计量
- 文章访问数: 1194
- HTML全文浏览量: 344
- PDF下载量: 189
- 被引次数: 11
目录
引言
1. 面内变刚度薄板弯曲理论
2. 神经网络模型
2.1 神经网络模型的构建
2.2 误差函数
2.3 学习率选取方案
2.4 算法流程
3. 算例
3.1 双向面内变刚度圆形薄板轴对称弯曲
3.2 单向面内变刚度方形薄板受非线性荷载作用
3.3 三角形面内变刚度薄板受线性分布荷载作用
4. 本文方法计算效率及其有限元对比分析
5. 网络耦合系数
$ {k}_{{\rm{p}}} $ 、隐藏层层数与神经元个数、激活函数对网络收敛性的影响分析5.1 网络耦合系数
${k_{\rm{p}}}$ 5.2 隐藏层层数与神经元个数
5.3 激活函数
6. 结论与展望
