贺 云 杜 娟 李海滨
*(内蒙古农业大学水利与土木建筑工程学院，呼和浩特 010018)
†(内蒙古工业大学理学院，呼和浩特 010051)
**(内蒙古财经大学统计与数学学院，呼和浩特 010070)
轴对称结构裂纹弹性力学分析是工程实践中重要而基础的问题。传统的解析求解方法有按位移求解、按应力求解和混合求解[1]。这三种方法都要从静力学方面、几何学方面和物理学方面考虑。根据弹性理论可知，弹性力学问题可以归结为具有边界条件的偏微分方程求解问题。按应力求解该问题时，传统的求解方法有逆解法和半逆解法[2]。此类方法都需要针对所要求解的问题，假设应力函数为参数未知的某种函数形式，将其代入应力函数表示的应力分量表达式中。根据相容方程及边界条件来考察假设的应力函数，得到假定应力函数的待求参数，从而得到应力函数及应力分量。但由于弹性体裂纹形状及边界条件的复杂性，一般很难得到弹性力学问题的解析解。众所周知，己知描述物理问题的微分方程模型中，只有极少数部分是可以求出解析解，而绝大多数只能借助近似数值方法求解，弹性力学模型亦不例外，因此寻找有效数值解法就显得尤为重要，而有限单元法是迄今为止在科学与工程计算领域最活跃、应用最广泛的数值计算方法[3-4]。
基于基本的有限单元法，后来一些学者又提出了一些发展、改进的方法，如稳定混合有限单元法[5]、h型混合变阶有限单元法[6]、无闭锁非协调有限单元法[7]，自适应有限单元法[8]、间断有限单元法[9]等。但在实际应用中，有许多因素使得弹性力学问题产生强奇性解，如物理区域是一个非光滑的几何区域（凹角或裂缝等）、材料常数不连续、应力集中等。虽然这些特殊弹性力学问题在进行有限单元法求解时可以通过加密网格来克服求解困难，但是网格加密会引起计算量急剧增加。这是因为由网格加密所带来的自由度成指数增长，导致计算量成指数增加。
由于轴对称结构裂纹弹性分析可以归结为求解具有边界条件的偏微分方程，神经网络方法在此方面也有广泛的应用。Lagaris等[10]研究了通过神经网络求解含有初始边界条件的偏微分方程，初始边界条件的神经网络不含可调参数，偏微分方程为一个前馈型神经网络，将此方法的求解结果与有限单元法进行对比，在得到良好结果的同时，执行速度有了较大的提升。针对含初始边界条件的偏微分方程与神经网络结构及训练集无法有机结合在一起的问题，Aarts等[11]提出了一种基于遗传算法的具有特殊结构神经网络的求解偏微分方程方法，得到了较好的结果。张永恒等[12]通过基于Hopflied神经网络的优化算法来实现打靶法的过程，从而求解了具有初始边界条件的微分方程，此方法具有对初值选择要求较低，而计算精度较高等优点。柴俊霞等[13]设计了以多元二次径向基神经网络为求解单元的偏微分方程求解方法，并通过算例验证了方法的有效性。Mosleh等[14]运用模糊神经网络求解具有初始值问题的模糊微分方程，并在一个工程问题中加以运用。Guner等[15]运用exp函数方法求解分数阶微分方程，同样得到了较好的结果。然而，传统的神经网络[16-17]求解偏微分方程的难点在于，单个神经网络无法同时满足偏微分方程及其边界条件，而训练多个神经网络又难以体现神经网络之间的联系，且传统神经网络求解偏微分方程问题工作效率低、精度无法保证。
为此，本文在弹性体模型解析求解方法的基础上，拟提出利用多神经网络（multi-neural network，MNN）联合训练方法实现轴对称结构裂纹弹性力学分析。首先，对轴对称结构裂纹进行弹性力学描述，由于按应力求解弹性力学应力问题可以简化为应力函数表示的相容方程及其边界条件的求解问题，所以针对该问题，随后将应力函数假设为统一的神经网络形式，构造应力函数表示的相容方程及应力边界条件的神经网络结构，从而组成神经网络系统。运用MNN联合训练方法训练该神经网络系统，从而得到应力分量的神经网络数值解。最后采用两个典型算例的计算和对比分析验证所提方法的有效性。
外部圆形裂纹，裂纹半径为c，无限弹性体任意半径为ρ，转过任意角度为φ，如图1所示。本文将在弹性力学理论基础上，用神经网络方法研究无限弹性体中含有一个外部圆形裂纹的轴对称力学问题。对于上述问题，当裂纹上下表面作用的外应力自相平衡时，可化归为具有应力边界的平面轴对称结构问题，其精确解可由弹性力学理论求得。对于受非对称载荷作用的复杂边界情况，拟给出裂纹面的神经网络求解方法。
图1 外部圆形裂纹Fig.1 External circular crack
根据轴对称结构弹性理论，可以归结为求解平衡微分方程、几何方程和物理方程的方程组。弹性体的平衡微分方程可以表示为
式中，σρ和σφ为极坐标系下的正应力分量，τρφ为极坐标系下的切应力分量，fρ和fφ为极坐标系下的体力分量。
几何方程为
式中，ερ和εφ为极坐标系下的正应变分量；
γρφ为极坐标系下的切应变；
uρ和uφ为极坐标系下的位移分量。
对于应力问题，物理方程为
式中，E为弹性模量，µ为泊松比。
至此，弹性力学应力问题可以转化成求解满足具有边界条件的相容方程。该相容方程由应力函数表示，通过应力函数的表达式，从而可以计算得到应力分量。
对于弹性体模型应力问题，根据弹性理论，用径向坐标ρ和环向坐标φ表示的极坐标系下相容方程为
应力边界条件为
当不计体力时，通过应力函数φ可以求得极坐标系下应力分量为
对于包含未知函数及其偏导数的等式，一般说来，如果x1,x2,...,xn为变量，以u˜ 为未知函数，对任意变量xi最高阶偏导数不超过Ai的线性偏微分方程的一般形式可表达为
在偏微分方程问题中往往还带有边界条件，可表示为
由神经网络相关理论可知，一个多层神经网络可以逼近任意非线性函数，这对于无法用解析方法求解及数值求解较困难的偏微分方程是一种思路。因此，可以将偏微分方程的通解u˜ 写为一个参数待定的带权数的隐式多项式
根据如上理论，轴对称结构应力函数φ(ρ,φ)表示为变量ρ和φ的带权数的隐式多项式，此式结构和神经网络表达式结构类似，即
网络隐层单元激活函数f(x)=ex， 将式（15）代入式（12）～式（14），得到各应力分量表达式为
将式（15）代入式（9），并将式（16）～式（18）代入式（10）和式（11），得到极坐标下偏微分方程及其边界条件的神经网络表达式为
式（19）～式（21）就构成了具有边界条件偏微分方程的神经网络求解表达式。比较式（19）与式（20）、式（21）对应的神经网络可知，各网络在结构上均为两输入、单输出、n个隐层单元的三层前馈型神经网络。网络结构如图2～图4所示。
图2 式（19）对应神经网络结构Fig.2 The neural network structure of Eq.(19)
图4 式（21）对应神经网络结构Fig.4 The neural network structure of Eq.(21)
根据式（19）～式（21）构造极坐标系下的神经网络系统。网络输入均为ρ和φ，网络输出分别为h2 ，f¯ρ，f¯φ。此外，各网络的输入层到隐层单元连接权值为kij，其中j=1,2 ，隐层单元阈值为bi。式（19）～式（21）对应的神经网络与式（15）对应的网络隐层到输出层单元连接权值之间存在的关系为
欲通过学习使上述三个神经网络同时收敛，则三个神经网络构成网络系统的总性能函数误差收敛，其中总性能函数可定义为式（25）。
神经网络系统中ki1，ki2，bi，wi为独立变量，令x=[ki1,ki2,bi,wi]T为网络系统待调节列向量。系统性能函数写为
系统性能函数E˜ 为描述神经网络系统整体误差的函数，其值越小，神经网络计算精度越高，隐层神经元激活函数为f(x)=ex，文献[18]验证了隐层单元激活函数为指数函数时，较logsig和tansig等传统激活函数，同样具有较好的计算精度和泛化能力。
图3 式（20）对应神经网络结构Fig.3 The neural network structure of Eq.(20)
由于运用梯度下降法训练网络时，稳定性要求学习率很小，所以使得训练很慢。动量梯度法速度虽然有所提高，但还是无法达到实际应用的要求。牛顿法收敛速度快，但是在每一次迭代中，需要求出Hessian矩阵，其中包括性能函数的二次导数，这就使得计算量变得很大。还有一些神经网络的改进训练算法[19]，如动量梯度下降法、牛顿法、共轭梯度法、 Levenberg-Marquardt（L-M）算法等。
如果神经网络的性能函数为平方和的形式，就可以采用L-M训练算法。此算法是神经网络计算中较为常用的算法，不需要计算Hessian阵，Hessian阵可以用矩阵H=JTJ来近似替换。其中，J为网络性能函数E˜ 对权值、阈值一阶导数的雅可比矩阵，它是关于权值、阈值的函数。该矩阵的求解相比Hessian矩阵求解要简单很多。
L-M训练算法中网络参数通过式xr+1=xr-[JTJ+µI]-1JTE˜ 进行学习，式中I为单位矩阵，µ为学习率。
L-M训练算法是牛顿法与梯度下降法的结合，当学习率µ很小时，相当于牛顿法的迭代步长。当µ很大时，相当于梯度下降法的迭代步长。在迭代过程中，如果训练成功，就减小µ的值；
如果训练失败，就增加µ的值，这样使性能函数的值逐渐减小。
因此，本文方法结合L-M训练算法，运用上述MNN联合训练方法训练神经网络系统，得到神经网络参数。将网络参数代入式（16）～式（18），从而得到极坐标系下应力分量的数值解。
3.1 算例1
圆环分别受到均布外压力qb= 2 kN/m，均布内压力qa= 1 kN/m。外径R= 3 cm，内径r= 2 cm，如图5，求模型中各点的应力分量。
图5 圆环结构Fig.5 Ring structure
对于圆环轴对称问题，相容方程可简化为
边界条件为
对于该问题，由于应力函数φ只是半径ρ的函数，所以φ的神经网络形式可简化为
将式（29）代入相容性方程及边界条件，构造神经网络结构及训练样本。由于边界条件式（27）自动满足，所以训练样本如表1所示。
表1 神经网络系统训练样本Table 1 Training samples of neural network system
隐层单元激活函数 ，经过调节，算例中隐层单元个数n= 10。运用MNN联合训练方法，联合训练神经网络系统20步，网络系统收敛，误差收敛曲线如图6所示。
图6 算例1误差收敛曲线Fig.6 Error convergence curve in example 1
根据简化后的应力函数神经网络形式式（29），得到径向应力分量和环向应力分量的表达式为
圆环半径在[2,3]范围内均匀取11个点作为测试样本，将其代入式（30）和式（31），得到测试样本点处应力分量的神经网络数值解，将径向应力的数值列于表2中。
根据弹性力学解析求解方法得到径向应力的表达式为
通过式（32）得到该平面问题的弹性力学解析解作为理论解。为了进一步对比结果，由于圆环的对称性，取圆环的1/4进行分析，通过Ansys软件求得该问题的有限单元法解。图7为计算精度和时间随网格数量的变化曲线，曲线1表示结构中计算精度随网格数量收敛的一般曲线，曲线2表示计算时间随网格数量的变化。从图中可看出网格较少时，增加网格数量可以使计算精度明显提高，而计算时间不会有大的增加。当网格数量增加到一定程度后，再继续增加网格时，精度提高甚微，而计算时间却有大幅度增加。所以为了增加网格的经济性，算例中取计算精度大致收敛时的网格数量，P= 160个。两种方法的计算结果一同列于表2。径向应力的神经网络方法解与理论解对比如图8所示，有限单元法解的径向应力云图如图9所示。
图7 有限单元法计算精度和时间曲线Fig.7 Calculation accuracy and time curves of finite element method
图9 有限单元法的径向应力云图Fig.9 Radial stress cloud figure of finite element method
表2 径向应力计算结果Table 2 Calculation results of radial stress
由表2结果可以看出，本文方法总体的计算精度高于有限单元法。并且从图8中可以看出，在越接近边界处，本文方法的精度越高，而有限单元法在边界处的计算精度却较差。
图8 径向应力对比图Fig.8 Comparition figure of radial stress
3.2 算例 2
矩形薄板长a= 8 cm，高b= 7 cm，在离开边界较远处有半径为r=0.5cm 的小圆孔规则裂纹。上下及左右边界分别受到均布压力和拉力，大小均为q= 1 kN/m，如图10所示，求图中半径为R的大圆内各点的应力分量。
图10 带圆孔的矩形薄板结构Fig.10 Structure of rectangular thin plate with round hole
该问题为孔口应力集中问题，具有局部性，一般孔口的应力集中区域约在距孔边1.5倍圆孔直径的范围内。为了不失一般性，本算例取距孔边为2.5倍圆孔直径的大圆（虚线所示）范围内进行分析。在大圆外部，应力情况与无孔时基本相同，本算例将不作计算。该问题的相容方程为式（9）。
边界条件为
运用本文方法构造网络结构及训练样本，式（9）、式（33）、式（34）分别对应的神经网络训练样本如表3～表5所示。
表3 式（9）对应的网络训练样本Table 3 Training samples of Eq.(9) network
表4 式（33）对应的网络训练样本Table 4 Training samples of Eq.(33) network
表5 式（34）对应的网络训练样本Table 5 Training samples of Eq.(34) network
隐层单元激活函数f=ex，经过调节，算例中隐层单元个数n= 50。运用MNN联合训练方法，联合训练神经网络系统20步，网络系统收敛，误差收敛曲线如图11所示。
图11 算例2误差收敛曲线Fig.11 Error convergence curve in example 2
将测试样本代入式（16）～式（18）， 得到测试样本点处的径向应力σρ，环向应力σφ和切应力τρφ，如表6～表8所示。
表6 测试样本点处径向应力 σρ (单位：kPa)Table 6 Radial stress σ ρ of test sample points (unit: kPa)
表8 测试样本点处切应力 (单位：kPa)Table 8 Shear stress of test sample points (unit: kPa)
表8 测试样本点处切应力 (单位：kPa)Table 8 Shear stress of test sample points (unit: kPa)
images/BZ_168_343_2726_422_2759.png images/BZ_168_1338_2693_1430_2727.png0 0.157 1 0.314 2 …… 1.413 7 1.570 8 0.500 0 -0.027 0 -0.021 0 -0.067 6 …… -0.068 6 -0.011 2 0.750 0 -0.005 5 -0.359 4 -0.703 1 …… -0.390 9 -0.118 5 1.000 0 -0.007 2 -0.398 8 -0.711 0 …… -0.417 7 -0.121 7…… …… …… …… …… …… ……2.750 0 -0.001 9 -0.360 2 -0.607 0 …… -0.271 2 -0.113 5 3.000 0 -0.003 9 -0.328 5 -0.605 5 …… -0.310 1 -0.128 2
根据本文方法得到的径向应力σρ，环向应力σφ和切应力τρφ在测试样本点处的大小分布如图12所示。
图12 算例2本文方法的应力分布图Fig.12 Stress distribution figure of the proposed method in example 2
为了便于进行对比，利用有限单元法求解该问题，根据图7的精度收敛曲线，本算例取网格数量为P= 135个，并在应力集中处进行网格加密处理。在Ansys软件中的应力云图如图13所示。
表7 测试样本点处环向应力 (单位：kPa)Table 7 Hoop stress of test sample points (unit: kPa )
表7 测试样本点处环向应力 (单位：kPa)Table 7 Hoop stress of test sample points (unit: kPa )
images/BZ_168_346_2032_425_2066.png images/BZ_168_1341_2000_1433_2033.png0 0.157 1 0.314 2 …… 1.413 7 1.570 8 0.500 0 -4.031 2 -3.744 3 -3.313 7 …… 3.758 5 4.023 0 0.750 0 -1.621 4 -1.441 0 -1.252 2 …… 1.468 7 1.646 9 1.000 0 -1.166 9 -1.188 0 -1.007 6 …… 1.053 7 1.157 5…… …… …… …… …… …… ……2.750 0 -1.061 4 -0.903 2 -0.822 8 …… 0.904 8 0.990 9 3.000 0 -0.945 4 -0.954 8 -0.760 3 …… 0.920 8 0.930 2
根据弹性力学半逆解法，得到应力分量的求解表达式为
通过式（35）～式（37）可以求得该平面问题应力分量的弹性力学解析解作为理论解。有限单元法、本文方法求得的数值解与理论解相比的平均绝对误差列于表9。
表9 测试样本点处应力的平均绝对误差Table 9 Average absolute error of test sample points stress
由表9结果可以看出，对于径向应力和切应力，有限单元法和本文方法的计算误差相差不大。而对于应力集中较为严重的环向应力，本文方法的计算精度高于有限单元法。与传统有限单元法求解时，需要对物理模型进行网格划分相比，本文方法从弹性力学的基本方程出发，引入边界条件方程可以直接求解。
3.3 算例 3
无限大弹性体，在上边界处有半径为a=0.2cm的半圆孔裂纹。上边界受到相反的两个力，大小均为F= 10 kN，如图14所示，分析图中半径为b=0.4cm的半圆内各点的应力分量。
图14 带半圆孔裂纹的弹性体结构Fig.14 Structure of elastomer with semi-circular hole crack
根据相容方程式（9）和式（38），式（39）的边界条件，构造MNN结构，取对称结构的1/2进行分析，运用本文方法计算得到应力分量，如图15所示。
图15 算例3本文方法的应力分布图Fig.15 Stress distribution figure of the proposed method in example 3
利用有限单元法求解该问题，根据精度收敛曲线，本算例取网格数量为P= 1 400个，网格划分如图16所示。
图16 有限元网格划分图Fig.16 The figure of finite element mesh division
应力分量的理论解求解表达式为
有限单元法、本文方法求得的数值解与理论解相比的平均绝对误差列于表10。
由表10结果可以看出：（1）相比算例2，网格划分加密后，有限单元法计算精度略有提高；
（2）对于理论解表达式非线性较强时，本文方法精度略有降低，但是对于径向应力和切应力，两种方法的计算误差相差不大；
（3）对于应力集中较为严重的环向应力，本文方法在计算精度方面的优越性明显高于有限单元法。
表10 测试样本点处应力的平均绝对误差Table 10 Average absolute error of test sample points stress
3.4 算例4
载人潜水器的耐压球壳直径设计为2.1 m，壳体上有一个小端直径为200 mm的主观察窗，两个小端直径为120 mm的侧观察窗对称分布于主观察窗两侧。主侧观察窗的厚度为234 mm、均采用45o角锥台形结构，结构尺寸如图17所示，外部边界处的压强为162 MPa。
图17 观察窗结构图（单位：mm）Fig.17 Structure of observation window (unit: mm)
此轴对称结构观察窗和窗座的接触面可看作为一个外部圆形裂纹，由于此处容易出现应力集中，从而导致观察窗的结构破坏，该算例将对此处进行力学分析。根据相容方程和边界条件，构造MNN结构，取对称结构的1/2进行分析，运用本文方法计算得到应力分量。作为对比，利用有限单元法求解该问题，根据精度收敛曲线，本算例取网格数量为P= 400个，应力云图如图18所示，本文方法与有限单元法在测试点处的结果对比如图19所示。
图18 观察窗有限单元法应力云图Fig.18 Stress cloud of finite element method in observation window
测试样本点处有限单元法、本文方法求得的数值解及最大相对误差列于表11。
由表11结果可以看出：本文方法相比有限单元法，在样本点处径向应力、环向应力的最大相对误差均小于4%，从而说明对于复杂结构，本文方法有较好的适用性。从图19中可以看出：结构边界处两种方法相差较大，这是由于有限单元法在边界处的精度较差，而本文方法更接近于理论值。
表11 测试样本点处的应力Table 11 The stress of test sample points
图19 本文方法和有限元方法应力对比图Fig.19 Stress comparison between the proposed method and finite element method
本文提出了轴对称结构外部圆形裂纹弹性力学分析的神经网络方法。假设应力函数为神经网络形式，通过相容方程及边界条件与应力函数的导数关系，分别构造参数相关联的多个神经网络，通过联合训练方法训练网络系统，实现神经网络的高精度数值求解。数值算例表明，本文方法对包含应力集中在内的弹性力学分析有较好的适用性，在计算精度方面，相比有限单元法的计算结果有所提高。由于本文介绍的神经网络方法直接从理论求解弹性力学问题的偏微分方程出发，所以相比传统有限单元法的求解过程，本文方法在数值求解时无需对弹性体物理模型进行网格划分。
本文贡献在于：（1）将应力函数写成统一的神经网络形式，无需根据具体问题假设应力函数，从而克服了应力函数的求解困难；
（2）通过神经网络方法的求解，使轴对称结构外部圆形裂纹弹性力学分析的数值解精度更高，对复杂形状及复杂边界的弹性体力学分析提供了一种新的思路。另外，本文方法也可以拓展到按位移求解、混合求解的弹性力学分析中。