引言
面内功能梯度材料薄板结构在土木工程、海洋工程等领域应用广泛, 该结构由功能梯度材料所组成(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.
面内变刚度薄板弯曲理论
本文研究变厚度薄板或弹性模量参数在面内变化的薄板的弹性小变形弯曲, 设薄板的厚度函数为
根据Kirchhoff板理论基本假定, 几何方程为
$$ left. {begin{array}{*{20}{l}} {{varepsilon _x} = - zdfrac{{{partial ^2}w(x,y)}}{{partial {x^2}}}} {{varepsilon _y} = - zdfrac{{{partial ^2}w(x,y)}}{{partial {y^2}}}} {{gamma _{xy}} = - zdfrac{{{partial ^2}w(x,y)}}{{partial xpartial y}}} end{array}} ight} $$ | (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}} ight} $$ | (2) |
平衡方程为
$$ frac{{{partial ^2}{M_x}}}{{partial {x^2}}} + 2frac{{{partial ^2}{M_{xy}}}}{{partial xpartial 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} ight)z text{d}z=&qquad -frac{E(x,y)h{(x,y)}^{3}}{12left(1-{nu }^{2} ight)}left[frac{{partial }^{2}w(x,y)}{partial {x}^{2}}+nu frac{{partial }^{2}w(x,y)}{partial {y}^{2}} ight]end{split} $$ | (4) |
记弯曲刚度函数为
$$ D(x,y) = frac{{Eleft( {x,y} ight)h{{left( {x,y} ight)}^3}}}{{12left( {1 - {nu ^2}} ight)}} $$ | (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}}}} ight] $$ | (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}}}} ight] $$ | (7) |
$$ {M_{xy}} = int_{ - h/2}^{h/2} {{tau _{xy}}} z{ m{d}}z = - (1 - nu )D(x,y)frac{{{partial ^2}w(x,y)}}{{partial xpartial y}} $$ | (8) |
对于固支边界条件
$$ {left( {w,frac{{partial w}}{{partial n}}} ight)_s} = 0 $$ | (9) |
对于简支边界条件
$$ {left( {w,{M_n}} ight)_s} = 0 $$ | (10) |
对于自由边界条件
$$ left( {{M_n},F_{{{s}}n}^t} ight) = 0 $$ | (11) |
其中
将式(4) ~ 式(6)代入平衡方程(3)即可得面内变刚度薄板弯曲偏微分控制方程
$$ begin{split}& {nabla ^2}left( {D{nabla ^2}w} ight) - (1 - nu )left( {frac{{{partial ^2}D}}{{partial {y^2}}}frac{{{partial ^2}w}}{{partial {x^2}}}} ight. - 2frac{{{partial ^2}D}}{{partial xpartial y}}frac{{{partial ^2}w}}{{partial xpartial y}} + & qquad left. {frac{{{partial ^2}D}}{{partial {x^2}}}frac{{{partial ^2}w}}{{partial {y^2}}}} ight) + q = 0 end{split} $$ | (12) |
式中
2.
神经网络模型
2.1
神经网络模型的构建
本文方法并非直接设计网络来求解方程(12), 而是采用两个神经网络模型来进行求解, 如图1所示, 将待求解的4阶偏微分控制方程转换为求解4个二阶偏微分方程组, 该解法本质上仍属于强形式的求解方案. 如果仅以挠度作为预测解, 在试函数的构造上对于不同形状的求解域以及简支、自由边界条件的构造会出现困难. 本文采用挠度网络预测薄板挠度
m{w}}}}
ight)$
m{m}}}}
ight) = $
m{m}}}}
ight),{{bar M}_{xy}}left( {x,y;{{{{boldsymbol{theta }}}}_{
m{m}}}}
ight),{{bar M}_y}left( {x,y;{{{{boldsymbol{theta}} }}_{
m{m}}}}
ight)}
ight}$
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-1.jpg'"
class="figure_img
figure_type2 ccc " id="Figure1" />
图
1
本文神经网络模型示意图
Figure
1.
The schematic diagram of neural network model in this paper
下载:
全尺寸图片
幻灯片
2.2
误差函数
根据广义应力?应变关系式(6) ~ 式(8)可求得挠度二阶偏导的弯矩表达式
$$ frac{{{partial ^2}w}}{{partial {x^2}}} = frac{1}{{Dleft( {1 - {nu ^2}} ight)}}left( {nu {M_y} - {M_x}} ight) $$ | (13) |
$$ frac{{{partial ^2}w}}{{partial {y^2}}} = frac{1}{{Dleft( {1 - {nu ^2}} ight)}}left( {nu {M_x} - {M_y}} ight) $$ | (14) |
$$ frac{{{partial ^2}w}}{{partial xpartial y}} = - frac{{{M_{xy}}}}{{D(1 - nu )}} $$ | (15) |
误差函数的构造是神经网络训练的核心, 由于本文方法引入两个网络进行计算, 故在训练中需要考虑两者之间的耦合误差. 若采用无约束优化方案, 本文误差函数主要根据挠度与弯矩网络在边界处的误差、弯矩网络在力平衡方程(3)中的误差、预测的挠度与弯矩通过式(13) ~ 式(15)建立的耦合误差来构造.
采用均方误差(mean square error, MSE)来衡量神经网络的拟合误差, 记挠度网络与弯矩网络的内部参数分别为
m{w}}} $
m{m}}$
$$ begin{split} Cleft( {{{{{boldsymbol{theta}} }}_{ m{w}}},{{{{boldsymbol{theta}} }}_{ m{m}}}} ight) = & MS{E_{text{R}}} + {k_{ m{P}}}MS{E_{ m{P}}} + &{k_1}MS{E_{{varGamma _1}}} + {k_2}MS{E_{{varGamma _2}}} + {k_3}MS{E_{{varGamma _3}}} end{split} $$ | (16) |
如果采用构造试函数的形式使得边界误差强制满足, 则误差函数无需计算边界误差
$$ Cleft( {{{{{boldsymbol{theta}} }}_w},{{{{boldsymbol{theta}} }}_{ m{m}}}} ight) = MS{E_{text{R}}} + {k_{ m{P}}}MS{E_{ m{P}}} $$ | (17) |
其中
$$ begin{split} MS{E_{ m{R}}} =& frac{1}{{{N_varOmega }}}sumlimits_{i = 1}^{{N_varOmega }} {left| {frac{{{partial ^2}{{bar M}_x}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight)}}{{partial {x^2}}}} ight. + 2frac{{{partial ^2}{{bar M}_{xy}}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight)}}{{partial xpartial y}}} + & {left. {frac{{{partial ^2}{{bar M}_y}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight)}}{{partial {y^2}}} + q} ight|^2} end{split} $$ | (18) |
$$ MS{E_{{varGamma _1}}} = frac{1}{{{N_{{varGamma _1}}}}}sumlimits_{i = 1}^{{N_{{varGamma _1}}}} {left[{{{left| {bar wleft( {{{boldsymbol{x}}_{{varGamma _1}}};{{boldsymbol{theta }}_{ m{w}}}} ight)} ight|}^2} + left. {{{left| {frac{{partial bar wleft( {{{boldsymbol{x}}_{{varGamma _1}}};{{boldsymbol{theta }}_{ m{w}}}} ight)}}{{partial n}}} ight|}^2}} ight]} ight.} ;;;;;;$$ | (19) |
$$ MS{E_{{varGamma _2}}} = frac{1}{{{N_{{varGamma _2}}}}}sumlimits_{i = 1}^{{N_{{varGamma _2}}}} {left[ {{{left| {bar wleft( {{{boldsymbol{x}}_{{varGamma _2}}};{{boldsymbol{theta }}_{ m{w}}}} ight)} ight|}^2} + left| {{{left. {{{bar M}_n}left( {{{boldsymbol{x}}_{{varGamma _2}}};{{boldsymbol{theta }}_{ m{m}}}} ight) - {{tilde M}_n}} ight|}^2}} ight.} ight]} $$ | (20) |
$$ begin{split} MS{E_{{varGamma _3}}} = &frac{1}{{{N_{{varGamma _3}}}}}sumlimits_{i = 1}^{{N_{{varGamma _3}}}} {Biggl[{23} {{{left| {left. {{{bar M}_n}left( {{{boldsymbol{x}}_{{varGamma _3}}};{{boldsymbol{theta }}_{ m{m}}}} ight) - {{tilde M}_n}} ight|} ight.}^2}} } hfill + & {{{left| {frac{{partial {{bar M}_{ns}}left( {{{boldsymbol{x}}_{{varGamma _3}}};{{boldsymbol{theta }}_{ m{m}}}} ight)}}{{partial s}} + {Q_n}left( {{{boldsymbol{x}}_{{varGamma _3}}};{{boldsymbol{theta }}_{ m{m}}}} ight)} ight|}^2}} Biggl]{23} hfill end{split};;;;;; ;;;;;;;;;;;;;;;;;;$$ | (21) |
$$ MS{E_{ m{P}}} = MS{E_{{ m{P}}1}} + MS{E_{{ m{P}}2}} + MS{E_{{ m{P}}3}};;;;;;;;;;;;;;;;;;;;;;;;;;; $$ | (22) |
其中, 向量
ight)$
m{m}}}}
ight)$
m{m}}}}
ight)$
m{p}}}$
式(22)中
$$ begin{split} MS{E_{{ m{P}}1}} =& frac{1}{{{N_varOmega }}}sumlimits_{i = 1}^{{N_varOmega }} {left| {frac{{{partial ^2}bar wleft( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{w}}}} ight)}}{{partial {x^2}}}} ight.} hfill & {left. { - frac{1}{{Dleft( {1 - {nu ^2}} ight)}}left[ {nu {{bar M}_y}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight) - {{bar M}_x}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight)} ight]} ight|^2} hfill end{split} $$ | (23) |
$$ begin{split} MS{E_{{ m{P}}2}} = &frac{1}{{{N_varOmega }}}sumlimits_{i = 1}^{{N_varOmega }} {left| {frac{{{partial ^2}bar wleft( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{w}}}} ight)}}{{partial {y^2}}}} ight.} hfill &{left.{- frac{1}{{Dleft( {1 - {nu ^2}} ight)}}left[ {nu {{bar M}_x}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight) - {{bar M}_y}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight) } ight]} ight|^2} hfill end{split} $$ | (24) |
$$ MS{E_{{ m{P}}3}} = frac{1}{{{N_varOmega }}}sumlimits_{i = 1}^{{N_varOmega }} {{{left| {frac{{{partial ^2}bar wleft( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{w}}}} ight)}}{{partial xpartial y}} + frac{{{{bar M}_{xy}}left( {{boldsymbol{x}};{{boldsymbol{theta }}_{ m{m}}}} ight)}}{{D(1 - nu )}}} ight|}^2}} $$ | (25) |
本文的误差函数表达式中包含挠度、弯矩对自变量的二阶偏导数项, 对于这些偏导项, 一方面可以利用神经网络的输出构造差分求解格式来近似求解, 但采用该方案需要较大的计算量才能得到精确的计算结果; 另一方面, 基于计算图的自动微分技术(automatic differentiation, AD)可以高效地处理神经网络对输入变量求导过程, 当前的深度学习框架如Tensoflow, Pytorch, MindsSpore等均支持自动微分. 本文基于Pytorch提供的自动微分接口实现对上述偏导项及误差函数梯度的计算.
在实际计算中, 也可灵活采用混合边界误差的形式进行求解, 如对于部分简单的边界条件构造特解, 而对于复杂的边界采用相应的无约束优化方案. 建立本文的误差函数后, 在每个训练批次(epoch)中均需计算其梯度并结合误差反向传播算法更新网络的内部参数, 关于该过程, Tang和Yang[21]对其进行了详细的讨论.
2.3
学习率选取方案
学习率的选取可直接影响神经网络的训练, 目前神经网络学习率的选取仍带有一定的经验性, 但总体而言, 在训练初期选取较大的学习率可以加快误差收敛速度, 在训练后期, 此时神经网络模型已经学习到相应的特征, 此时往往需要降低学习率, 以便对神经网络内部参数进行微调, 使得误差波动幅度不至于过大. 经过本文的实践, 本文学习率选取方案如下
$$ {alpha _t} = left{ {begin{array}{*{20}{l}} 1times{{{10}^{ - 3}},}&{t leqslant 5000} {5 times {{10}^{ - 4}},}&{5000 < t leqslant 10;000} 1times{{{10}^{ - 4}},}&{10;000 < t leqslant 20;000} {5 times {{10}^{ - 5}},}&{20;000 < t leqslant 30;000} 1times{{{10}^{ - 5}},}&{30;000 < t leqslant 40;000} {5 times {{10}^{ - 6}},}&{40;000 < t leqslant 45;000} 1times{{{10}^{ - 6}},}&{45;000 < t} end{array}} ight. $$ | (26) |
其中,
2.4
算法流程
算法1. 本文算法
输入: 学习率
输出: 挠度网络; 弯矩网络
初始化神经网络模型及其内部参数
m{w}}},{{{{boldsymbol{theta}} }}_{
m{m}}}$
1,
2 while (
{
(1)在求解域中随机生成数据点坐标
m{d}}},{y_{
m{d}}}}
ight) in varOmega$
m{b}}},{y_{
m{b}}}}
ight) in$
(2)将数据点坐标作为参数输入神经网络模型, 通过前向传播算法, 求得
m{w}}}}
ight)$
m{m}}}}
ight)$
m{m}}}}
ight)$
m{m}}}}
ight)$
(3)通过式(16)计算误差函数
(4)计算误差函数
m{w}}},{{{{boldsymbol{theta}} }}_{
m{m}}}$
m{w}}},{{{{boldsymbol{theta}} }}_{
m{m}}}$
(5)
/}
算法中, 本文的数据点的默认生成参数为
m{w}}},{{{{boldsymbol{theta}} }}_{
m{m}}}$
3.
算例
3.1
双向面内变刚度圆形薄板轴对称弯曲
图2所示受横向均布荷载
ho ) = {D_0}{{text{e}}^{ - m
ho /a}}$
ho = sqrt {{x^2} + {y^2}} $
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-2.jpg'"
class="figure_img
figure_type1 bbb " id="Figure2" />
图
2
圆形面内变刚度薄板
Figure
2.
Circular thin plate with in-plane stiffness gradient
下载:
全尺寸图片
幻灯片
本文选取梯度参数m分别为0, 0.5, 1, 2的情况进行计算. 本算例在计算过程中仅需要施加位移边界条件, 设
m{w}}}}
ight)$
m{w}}}}
ight) $
ight)^2}$
本算例的挠度模型以及弯矩模型均采用具有6层隐藏层, 每层隐藏层具备30个神经元的网络结构(记为6 × 30), 采用
m{p}}} = 100$
每个梯度参数下神经网络模型的计算误差如图3所示, 由此可见, m = 0时为刚度恒定的薄板, 此时误差函数收敛较快, 相比其他工况其误差最终的收敛值最小; 随着梯度系数的增大, 训练误差最终的收敛值出现增大的趋势.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-3.jpg'"
class="figure_img
figure_type1 bbb " id="Figure3" />
图
3
神经网络训练误差曲线
Figure
3.
The convergence curve of neural network training error
下载:
全尺寸图片
幻灯片
采用无量纲计算公式
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-4.jpg'"
class="figure_img
figure_type1 bbb " id="Figure4" />
图
4
PINN求解圆形面内变刚度薄板弯曲问题的训练误差收敛曲线(m = 2)
Figure
4.
Training error convergence curve of PINN (m = 2)
下载:
全尺寸图片
幻灯片
表
1
本文方法计算
Table
1.
Comparison of dimensionless
table_type1 ">
m | (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
本文方法计算
Table
2.
Comparison of dimensionless
table_type2 ">
m | (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
|显示表格
为了说明本文方法在求解面内变刚度薄板弯曲问题上的优点, 本算例也利用PINN来求解其四阶偏微分控制方程(12), 采用隐藏层层数为6, 每层隐藏层具有30个神经元的神经网络模型对工况m = 2进行求解, 激活函数为Tanh函数, 训练的数据点由求解域中随机生成, 数据点的产生有两种方案:
方案(1)为在整个求解域中随机生成训练数据点, 此时的数据点可在原点附近生成;
方案(2)为在求解域中随机生成的数据点但离原点较远. 此时两方案在训练过程中的误差收敛情况如图4所示, 可见采用相同的模型, 而生成的数据点不同则会导致模型的训练出现不同的结果, 虽然采用数据点生成方案(1)的模型训练也收敛, 但由于其误差此时收敛于一个较大的值, 得不到正确解. 经过本文分析, 这主要是由于本算例的弯曲刚度函数D的二阶导数在靠近原点区域出现“爆炸”式变化的原因, 即刚度函数的二阶偏导在原点处不收敛, 在靠近原点处
3.2
单向面内变刚度方形薄板受非线性荷载作用
如图5所示边长为
ight)}}$
本算例选取6 × 30的挠度网络模型, 5 × 50的弯矩网络模型, 激活函数选择Tanh函数,
m{p}}} = $
m{w}}}}
ight)$
m{m}}}}
ight) = left{ {{M_x},{M_{xy}},{M_y}}
ight}$
$$ left. begin{array}{l} {w^*} = wleft( {x,y;{{boldsymbol{theta }}_{ m{w}}}} ight){left[left(sqrt 3 ax + y - a ight)left( - sqrt 3 ax + y - a ight) ight]^2} hfill left{ {{M_x}^*,{M_{xy}}^*,{M_y}^*} ight} = left{ {{M_x}(x - 1),{M_{xy}},{M_y}} ight} hfill end{array} ight} $$ |
如果将本算例的应力边界条件也通过挠度网络进行构造, 则挠度表达式将复杂许多.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-5.jpg'"
class="figure_img
figure_type1 bbb " id="Figure5" />
图
5
方形面内变刚度薄板
Figure
5.
Thin square plates with in-plane stiffness gradient
下载:
全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-6.jpg'"
class="figure_img
figure_type1 bbb " id="Figure6" />
图
6
神经网络训练误差曲线
Figure
6.
Neural network training error-curve
下载:
全尺寸图片
幻灯片
将本算例的计算结果与有限元解答对比, 有限元计算中每个单元的弯曲刚度根据单元的形心坐标计算, 采用50
由挠度计算图7可知, 在挠度的求解上, 本文解法与有限元解法一致. 由图8的弯矩对比图可发现, 本文弯矩解与有限元解答基本吻合, 而当梯度系数m = 2时, 虽然本文解与有限元解在部分点上弯矩的相对误差增大, 但整体上解答吻合.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-7.jpg'"
class="figure_img
figure_type1 bbb " id="Figure7" />
图
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)
下载:
全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-8.jpg'"
class="figure_img
figure_type1 bbb " id="Figure8" />
图
8
不同梯度参数下本文弯矩
Figure
8.
Comparison of dimensionless bending moment
下载:
全尺寸图片
幻灯片
3.3
三角形面内变刚度薄板受线性分布荷载作用
工程中地下掩体结构有时可看作三角形薄板结构. 如图9所示, 厚度为
ight)}
ight]$
ight)}}$
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-9.jpg'"
class="figure_img
figure_type1 bbb " id="Figure9" />
图
9
三角形面内变刚度薄板
Figure
9.
Thin triangular plate with in-plane stiffness gradient
下载:
全尺寸图片
幻灯片
本算例的挠度网络与弯矩网络分别采用隐藏层层数为6和5, 每层分别具备30和50个神经元的神经网络结构进行计算, 两个模型均采用
m{p}}} = 100$
本算例的挠度、弯矩试函数为
$$ {begin{array}{*{20}{l}} {{w^*} = wleft( {x,y;{{{{boldsymbol{theta}} }}_{ m{w}}}} ight){{left[left(sqrt 3 ax + y - a ight)left( - sqrt 3 ax + y - a ight) ight]}^2}} left{ {{M_x}^*,{M_{xy}}^*,{M_y}^*} ight} = hfill qquad left{ {{M_x}left( {x,y;{{{{boldsymbol{theta}} }}_{ m{m}}}} ight),{M_{xy}}left( {x,y;{{boldsymbol{theta }}_{ m{m}}}} ight)y,{M_y}left( {x,y;{{boldsymbol{theta }}_{ m{m}}}} ight)y} ight} hfill end{array}} $$ |
此时自由边的应力边界条件并未精确满足, 故在误差函数中引入自由边的剪力误差项
$$ MS{E_{{varGamma _3}}} = frac{1}{{{N_{{varGamma _3}}}}}sumlimits_{i = 1}^{{N_{{varGamma _3}}}} {{{left| {frac{{partial {M_y}left( {x,y;{{boldsymbol{theta}} _{ m{m}}}} ight)}}{{partial y}} + frac{{partial {M_{xy}}left( {x,y;{{boldsymbol{theta}} _{ m{m}}}} ight)}}{{partial x}}} ight|}^2}} $$ |
根据式(16), 误差函数为
$$ Cleft( {{{boldsymbol{theta }}_w},{{boldsymbol{theta }}_{{m}}}} ight) = MS{E_varOmega } + {k_{ m{p}}}MS{E_{ m{p}}} + {k_3}MS{E_{{varGamma _3}}} $$ |
神经网络训练收敛后, 采用无量纲计算公式
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-10.jpg'"
class="figure_img
figure_type1 bbb " id="Figure10" />
图
10
三角形面内变刚度薄板沿轴线x = 0上的挠度
Figure
10.
Dimensionless deflection variation of thin triangular plate with in-plane stiffness gradient along axis x = 0
下载:
全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-11.jpg'"
class="figure_img
figure_type2 ccc " id="Figure11" />
图
11
m = 0.2时三角形面内变刚度薄板弯矩的有限元计算结果
Figure
11.
Finite element calculation of bending moment of thin plate with triangular in-plane variable stiffness (m = 0.2)
下载:
全尺寸图片
幻灯片
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-12.jpg'"
class="figure_img
figure_type2 ccc " id="Figure12" />
图
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.
本文方法计算效率及其有限元对比分析
以
本文的有限元求解程序采用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 paper
table_type2 ">
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
table_type1 ">
Numerical example | m | Number of nodes |
2 | 0 | 441 |
1 | 841 | |
2 | 1681 | |
3 | 0 | 305 |
0.2 | 1331 | |
0.5 | 1604 |
下载:
导出CSV
|显示表格
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-13.jpg'"
class="figure_img
figure_type1 bbb " id="Figure13" />
图
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}_{{
m{p}}} $ 、隐藏层层数与神经元个数、激活函数对网络收敛性的影响分析
m{p}}} $
5.1
网络耦合系数${k_{
m{p}}}$
m{p}}}$
本文模型由于引入两个神经网络模型对面内变刚度薄板弯曲问题进行计算, 与采用单网络的神经网络方法相比, 在训练过程中需考虑网络之间的耦合误差, 本文引入网络耦合系数
m{p}}} $
m{p}}}$
m{p}}}$
m{p}}}$
m{p}}}$
薄板变形能计算公式为
$$ {W_{{ m{int}}}} = frac{1}{2}int_varOmega {{{boldsymbol{k}}^{ m{T}}}} {boldsymbol{M}}{text{d}}varOmega = frac{1}{2}int_varOmega {{{boldsymbol{k}}^{ m{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}} ight] $$ | (28) |
$$ {boldsymbol{k}} = left{ {begin{array}{*{20}{c}} {{k_x}} {{k_y}} {{k_{xy}}} end{array}} ight} = - left{ {begin{array}{*{20}{c}} {dfrac{{{partial ^2}w}}{{partial {x^2}}}} {dfrac{{{partial ^2}w}}{{partial {y^2}}}} {2dfrac{{{partial ^2}w}}{{partial xpartial y}}} end{array}} ight} $$ | (29) |
外力所作实功为
$$ {W}_{text{ext }}=frac{1}{2}{displaystyle {int }_{varOmega }q}w text{d}varOmega $$ | (30) |
其则根据挠度网络的输出值计算.
采用高斯积分, 由上述公式计算神经网络训练时薄板的变形能、应变能的变化情况, 根据挠度网络输出计算外力实功
m{ext}}}}$
m{int}}}}$
m{p}}}$
m{p}}}$
m{p}}}$
m{p}}}$
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-14.jpg'"
class="figure_img
figure_type1 bbb " id="Figure14" />
图
14
神经网络训练时不同
m{p}}}$
Figure
14.
Changes in the deformation energy and the external force work during the training process with different
m{p}}}$
下载:
全尺寸图片
幻灯片
5.2
隐藏层层数与神经元个数
隐藏层层数与每层的神经元个数的选取均会影响到神经网络的训练效率, 不失一般性, 在其他计算参数不变的情况下, 本节选取算例2中m = 2的情况分别讨论隐藏层层数、神经元个数的改变对计算过程中误差函数收敛情况的影响, 计算结果如图15所示.
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-15.jpg'"
class="figure_img
figure_type1 bbb " id="Figure15" />
图
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函数等. 为讨论激活函数对神经网络训练的影响, 其余计算参数不变, 本节选取
onerror="this.onerror=null;this.src='https://lxxb.cstam.org.cn/fileLXXB/journal/article/lxxb/2021/9//21-273-16.jpg'"
class="figure_img
figure_type1 bbb " id="Figure16" />
图
16
不同的激活函数对计算误差的收敛影响对比
Figure
16.
Comparison of the effects of different activation functions on the convergence of loss function
下载:
全尺寸图片
幻灯片
6.
结论与展望
本文基于深度学习技术与强形式的求解方案建立了一种直角坐标下求解面内变刚度薄板弯曲问题的神经网络方法, 通过几个算例分析, 得出以下结论:
(1)本文解答与理论解、有限元解吻合, 证明了本文方法在求解面内变刚度薄板弯曲问题上的正确性, 本文的神经网络模型不需要对弯曲刚度函数求偏导, 其适应性更强, 同时在薄板的位移边界条件、应力边界条件的施加上较为方便.
(2)本文方法属于强形式的数值解法, 其计算所得结果具备连续性与可导性. 理论上, 本文方法可以求解弹性模量以及厚度在面内连续变化的薄板弯曲问题. 弯矩网络的求解受到梯度系数的影响, 在梯度变化较大处弯矩网络的求解精度受到一定的影响, 但对挠度网络的求解精度影响不大.
(3)由于神经网络方法为迭代类解法, 本文方法在薄板线性弯曲问题求解上的收敛速度较有限元慢, 但其计算所需内存较小. 通过本文的模型结构可看出, 神经网络方法具备相当大的灵活性, 根据这一特点, 可进一步发展求解面内变刚度功能梯度薄板非线性弯曲问题的神经网络方法, 神经网络方法在非线性问题的求解中具备潜在优势.
(4)本文模型仍存在优化空间: 一方面在本文模型的训练过程中, 误差函数及其梯度的计算在整个训练过程中占据大部分的时间, 可以考虑优化误差函数的构建过程, 如引入有限元中形函数的思想对算法进行优化; 另一方面为了使得挠度、弯矩网络具备较强的表达能力, 本文模型采用了两个具有独立参数的网络进行计算, 这导致了本文模型的训练参数较多, 为此后续优化中可将本文的两个网络合并为一个网络(2个输入, 4个输出), 对网络结构进行改良, 以减少训练参数.