有限单元法是采用计算机求解数学物理问题的一种数值计算方法，在各种结构的静动力分析、稳定性计算、场分析等领域中发挥了巨大的作用[4-5]。目前，基于冯·诺伊曼串行计算机的有限元算法往往都要花费很长的计算机时间，对于大型复杂结构则更是如此。考虑塑性、损伤、疲劳、蠕变等多种力学行为的全耦合求解问题，计算过程复杂，求解时间长，因此必须解决计算速度问题; 特别是对于一些时变系统要求实时分析和预测的场合更是要求计算具有实时性[6-7]。
从有限元计算步骤可以看出，对于大型复杂结构的有限元计算，整体刚度方程的求解运算在整个计算时间里占有很大的比重[4]。如何实现总体刚度矩阵的快速求解是缩短整个计算时间的关键之一。而卷积神经网络具有强大的映射、实时计算和并行计算的能力，是实现快速求解大型复杂结构问题的关键。总体刚度矩阵是由单元刚矩阵组集而成，快速求解单元刚度矩阵便成为了整个问题的核心，为此本文分析了单元刚度矩阵的特性和组成，提出了基于卷积神经网络的单元刚度矩阵高精度实时计算的方法。
本文方法并不是为了替代传统经典的有限元求解单元刚度矩阵的方法，而是为了探索一种将深度学习与传统结构分析结合的有效途径，为后续使用深度学习解决更复杂的结构问题，和将深度学习与传统有限元进行更深度融合提供一种方法参考。本文选择了典型的平面四边形单元为分析对象。
1 卷积神经网络及模型建立 卷积神经网络是一种前馈神经网络，与普通的神经网络区别在于：卷积神经网络包含了一个由卷积层和子采样层构成的特征提取器。在卷积神经网络的卷积层中，一个神经元只与部分邻层神经元相连。卷积神经网络有两大特点：①有至少一个卷积层，用来提取特征；②卷积层通过权值共享的方式进行工作，大大减少了权重的数量，使得在训练中达到同样的学习效果时，收敛速度明显快于全连接神经网络[8-9]。
卷积神经网络的基本训练学习过程如下：
步骤1卷积神经网络权值和阈值的初始化。阈值均初始化为0，使用Xavier初始化算法[10]对卷积神经网络各层的权值进行初始化。Xavier初始化算法是由输入和输出神经元的数目自动确定权值矩阵的初始化大小。卷积神经网络的权值参数初始化为0均值，并满足如下均匀分布：
 | (1) |
式中:ni为第i层神经元的数目；ni+1为第i+1层神经元的数目。
步骤2给定输入 x 和目标标签?。
步骤3前向传播，计算实际输出 y 。用各层权值矩阵 w l和阈值向量 b l，计算每层神经元的输出 a l，y 为卷积神经网络最后一层神经元的输出。
 | (2) |
 | (3) |
式中:上标l为卷积神经网络的层号；“*”代表卷积操作；a l-1为l-1层神经元的输出；w l为l-1层神经元与l层神经元之间的卷积核的权值矩阵；z l为l层神经元经过卷积操作之后的中间结果； b l为l层神经元的阈值; σ为神经元的激活函数。
本文卷积神经网络用到的激活函数为LeakyReLU[11]，其数学表达式为
 | (4) |
式中：α为(0, 1)区间内固定的参数，本文中设置α=0.1。函数图像如图 1所示。
 |
| 图 1 LeakyReLU激活函数 Fig. 1 LeakyReLU activation function |
| 图选项 |
步骤4损失函数。经过卷积神经网络的前向传播，可以得到输入样本所对应的预测值 y ，用损失函数来衡量样本预测值 y 与真实标签值?之间的差异。本文求解刚度矩阵为回归问题，要达到的目标为将误差控制在小数点后第三位的数字，平方误差会随着误差的减小而导致下降越来越慢，绝对值误差的下降速度随着误差的减小几乎保持不变，但是会出现导数在原点处不连续的情况。为了加快卷积神经网络的训练速度，并且保证导数连续，采用绝对值误差的改进型作为损失函数。
 | (5) |
式中:ε为(0, 1)之间非常小的数，是为了确保损失函数在0连续可导而添加的常量，本文取ε=10－12。
训练深度学习网络往往需要大量的数据，采用传统的梯度下降需要遍历所有数据才能求解出目标下降所需的梯度，对于数据量大时就会花费大量的时间。为了加快梯度的求解速度，采用了小批量梯度下降(Mini-Batch Gradient Descent, MBGD)法[12]。其次为了防止过拟合，在损失函数中加入了L2正则化[13]，将损失函数修正为
 | (6) |
式中:m0为小批量样本的个数；λ为正则化系数；
 | (7) |
步骤5反向传播，更新梯度，本文采用的是Adam[14]优化算法。
步骤6达到误差精度或迭代次数限制，则停止迭代，否则回到步骤2。
2 有限元总体刚度矩阵分析 结构的总体刚度矩阵 K 是由单元刚度矩阵 k 组装而成，公式为[15]
 | (8) |
式中: C e为单元选择矩阵。
因此，求解总体刚度矩阵的核心便是先获得精确的单元刚度矩阵。
单元刚度矩阵的计算实际上是从单元的尺寸、材料的性质到单元刚度矩阵的映射问题，而深度学习的一大优势便是建立输入与输出之间的复杂映射关系。对单元刚度矩阵进行分析，可以减少卷积神经网络不必要的输入和输出，减小网络的规模，缩短网络训练时间，同时也更容易获得一个比较好的学习效果。
设单元的形函数为Ni(x, y, z)(i=1, 2, …, m)，m为单元节点数，单元应变矩阵 B 为
 | (9) |
根据平面弹性力学等有限单元法，可知[13]：
 | (10) |
 | (11) |
 | (12) |
平面单元刚度矩阵为
 | (13) |
式中: D 为弹性矩阵；t为单元厚度。
由式(10)~式(12)可知，当平面单元在厚度保持不变时，节点坐标((x1, y1), (x2, y2), …, (xm, ym))等比例放大n倍，单元的形函数 B 会变为原来的1/n，而|J|会变为原来的n2倍。由式(13)可知，平面单元的刚度矩阵不会随着几何参数的缩放而改变。
由上述分析可知：
1) 对于平面单元，单元刚度矩阵不会随着几何参数整体缩放而改变。
2) 单元刚度矩阵为对称方阵[2]。
3) 单元刚度矩阵的每行元素之和为0[2]。
4) 从单元刚度矩阵的物理意义可以看出，其只与单元各节点的相对坐标、材料特性、形函数等有关，与节点的绝对坐标无关[3]。
3 卷积神经网络的实现 对于线弹性问题，单元刚度矩阵的计算实际上就是从单元的节点坐标到单元刚度矩阵元素的映射。
3.1 数据前处理 对于任意平面四边形单元(见图 2)，单元刚度矩阵只与单元各节点的相对坐标、材料特性、形函数等有关，与节点的绝对坐标无关, 且单元刚度矩阵并不会随着单元节点坐标的整体缩放而改变，特作如下前处理操作：
 |
| 图 2 平面四边形单元 Fig. 2 Planar quadrilateral unit |
| 图选项 |
1) 平移。将平面四边形单元的中心平移到原点。
2) 缩放。将平移后的单元节点坐标等比缩放，使得最大节点坐标值的绝对值为1。
经过平移和缩放之后的平面四边形单元如图 3所示。
 |
| 图 3 数据处理之后的平面四边形单元 Fig. 3 Planar quadrilateral unit after data processing |
| 图选项 |
前处理的目的为：①可以有效地减小输入单元节点坐标分布的离散程度，从而更有利于卷积神经网络学习由单元节点坐标到刚度矩阵元素之间的映射规律；②可以增加卷积神经网络在预测任意平面四边形单元刚度矩阵时的通用性。
3.2 确定卷积神经网络输入 由于卷积神经网络的输入是以图像为代表的RGB三通道矩阵数据，确定平面四边形单元节点坐标的输入数据格式便成为首要任务。
众所周知，彩色图像上的任何一个像素点都由RGB 3个通道的3个数据组成，这3个数据共同确定了该像素点的具体颜色。图 4为图像的RGB三通道数据存储方式。同理，对于组成结构的节点而言，每个节点在三维空间内的位置都由x、y、z三个方向的坐标数据组成，这3个数据共同确定了该节点的具体位置。对于平面四边形单元，单元的节点坐标由x、y两个方向上的坐标数据组成，这2个数据共同确定了该节点在平面中的具体位置。根据这一思想，可将平面四边形单元的节点坐标存储格式确定为图 5。
 |
| 图 4 图像的RGB三通道数据存储 Fig. 4 Data storage of a three-channel RGB picture |
| 图选项 |
 |
| 图 5 平面四边形单元的数据存储 Fig. 5 Data storage of planar quadrilateral unit |
| 图选项 |
该存储格式的优点为：有效保留了节点之间的相对位置信息，同时具有很强的扩展性。但是，由于2×2×2的矩阵数据输入对卷积神经网络而言太小，不利于加深卷积神经网络的深度，同时为了更好地利用单个节点的坐标信息，特通过补零操作将平面四边形单元的数据存储确定为图 6。
 |
| 图 6 补零操作后平面四边形单元的数据存储 Fig. 6 Data storage of planar quadrilateral unit after padding zero |
| 图选项 |
3.3 确定卷积神经网络输出 根据刚度矩阵的特性，即单元刚度矩阵为对称矩阵, 单元刚度矩阵的每行元素之和为零, 可知，只需求解出单元刚度矩阵的左上部分元素数值，便可推导出整个单元刚度矩阵的所有元素数值。对于平面四边形单元所需求解的刚度矩阵元素，见图 7中虚线围绕的元素。
 |
| 图 7 平面四边形单元所需求解的刚度矩阵元素 Fig. 7 Stiffness matrix elements to be solved for a planar quadrilateral unit |
| 图选项 |
求解单元刚度矩阵卷积神经网络的示意图如图 8所示。图中：输入为2@4×4，表示输入的数据为2×4×4的矩阵，补零操作是为了保证在卷积操作前后，数据的尺寸保持不变，这样可以加深卷积神经网络的深度，有效提升网络的学习能力，特征图为卷积核在数据上进行卷积操作之后所得到的结果，一个卷积核对应一个特征图，2×2卷积核表示卷积核尺寸的大小。
 |
| 图 8 卷积神经网络示意图 Fig. 8 Schematic diagram of convolutional neural network |
| 图选项 |
4 算例 为了方便讨论卷积神经网络的学习性能与卷积核数目、训练样本数目之间的关系，依据卷积神经网络的核心思想，即卷积神经网络首先提取少量的简单特征，然后通过不断地组合成更多的复杂特征，完成对原始数据中复杂含义的提取。每一个卷积核都对应一种特征组合，假定卷积神经网络的卷积核数目不断增加，卷积神经网络满足后一层卷积核的数目是前一层卷积核数目的2倍条件。
卷积神经网络的学习能力由测试数据误差表示。如果测试数据上的预测误差越小，则说明卷积神经网络的学习能力越强，预测的结果也就越准确；反之，测试数据上的预测误差越大，则说明卷积神经网络的学习能力越弱，预测的结果也就越不准确。
平面四边形单元的无量纲化基本信息为：单元厚度t=1, 弹性模量E=1, 泊松比μ=0.25。
算例在台式机上运行，对硬件的要求并不高，具体的硬件配置为：处理器Intel(R) Xeon(R) E2176M, 6核12线程CPU, 4GB内存。如果在专业的显卡上运行程序，则可以显著提升程序的运行速度。
4.1 卷积神经网络的学习效果与卷积核数目之间的关系 训练样本数目和卷积神经网络结构相同，仅改变卷积神经网络的卷积核数目时，对比卷积神经网络的学习能力与卷积核数目之间的关系。
卷积神经网络的参数设定为：学习率lr=0.001, 正则化系数λ=0，训练卷积神经网络的样本总数都为10000(其中80%为训练样本，20%为验证样本)。每次训练的样本数为100，迭代10000次以后，卷积神经网络的第一层卷积核数目从5, 10, 15，…，50依次递增，卷积神经网络学习的效果对比如图 9所示。
 |
| 图 9 预测误差与卷积核数目的关系 Fig. 9 Relationship between prediction error and number of convolution kernels |
| 图选项 |
从图 9可以看出，卷积神经网络在测试数据机上的预测误差，随着卷积核数目的增加而不断减小。卷积核数目增加到一定程度时，卷积神经网络学习效果提升的幅度不断减小，这是因为训练样本的数目保持不变，卷积神经网络的学习能力不断增强，训练样本就会显得不足，使得学习能力更强的卷积神经网络无法展现出应有的学习效果。
综上，在相同的卷积神经网络结构和训练样本数目情况下，增加卷积神经网络的卷积核数目可以有效提升网络的学习能力。
4.2 卷积神经网络的学习效果与训练样本数目之间的关系 卷积神经网络结构相同，仅改变训练卷积神经网络的样本数目，对比卷积神经网络的学习能力与训练样本数目之间的关系。
卷积神经网络第一层卷积核数目为40，学习率lr=0.0003, 正则化系数λ=0。小批量的样本数为100，迭代10000次以后，样本数目从5000, 10000, 15000, …, 50000依次递增(其中80%为训练样本，20%为验证样本)，卷积神经网络学习能力的对比如图 10所示。
 |
| 图 10 预测误差与训练样本数目的关系 Fig. 10 Relationship between prediction error and number of training samples |
| 图选项 |
从图 10可以看出，卷积神经网络在测试数据机上的预测误差，随着训练样本数目的增加而不断减小。训练样本数目增加到一定程度时，卷积神经网络学习能力提升的幅度几乎保持不变。说明不断增加训练样本数目，始终是提升预测准确度的有效途径。
综上，在相同的卷积神经网络结构和训练参数下，增加训练样本数目，可以有效提升卷积神经网络的学习效果。
上述分析为有效提升卷积神经网络的学习效果提供了依据：适量增大卷积神经网络的卷积核数目，同时增加训练样本数目。本文采用的卷积神经网络求解单元刚度矩阵的参数为：网络第一层卷积核数目为40，学习率lr=0.001, 正则化系数λ=0.0007，每次训练卷积神经网络的样本数为100，总样本数为100000。经过30000次迭代后，最终的训练结果为：单个样本上的载荷误差为0.016 5。每个样本共有28个要预测的元素数值，在单个元素数值上的误差为0.0006。
4.3 单元刚度矩阵计算结果对比 从表 1可以看出，用全连接神经网络求解单元刚度矩阵，误差在0.001以下，且其预测精度可以根据具体要求，通过调整卷积神经网络的卷积核数目和训练样本数目而不断优化，满足实际工程要求。
表 1 平面四边形单元所需求解单元刚度矩阵元素的真实值与预测值对比 Table 1 Comparison between real and predicted values of stiffness matrix elements to be solved for a planar quadrilateral unit
| 单元刚度矩阵元素 | 真实值 | 预测值 |
| k12 | 0.167 | 0.167 |
| k13 | -0.893 | -0.893 |
| k14 | -0.033 | -0.033 |
| k15 | -0.485 | -0.485 |
| k16 | -0.167 | -0.167 |
| k17 | 0.408 | 0.407 |
| k18 | 0.033 | 0.034 |
| k23 | 0.033 | 0.033 |
| k24 | -0.276 | -0.276 |
| k25 | -0.167 | -0.167 |
| k26 | -0.241 | -0.241 |
| k27 | -0.033 | -0.033 |
| k28 | 0.035 | 0.034 |
| k34 | -0.167 | -0.167 |
| k35 | 0.408 | 0.409 |
| k36 | -0.033 | -0.033 |
| k37 | -0.485 | -0.486 |
| k38 | 0.167 | 0.167 |
| k45 | 0.033 | 0.033 |
| k46 | 0.035 | 0.036 |
| k47 | 0.167 | 0.167 |
| k48 | -0.241 | -0.242 |
| k56 | 0.167 | 0.167 |
| k57 | -0.893 | -0.894 |
| k58 | -0.033 | -0.033 |
| k67 | 0.033 | 0.033 |
| k68 | -0.276 | -0.276 |
| k78 | -0.167 | -0.167 |
| 注：真实值的输入为经过平移和归一化操作之后的平面四边形单元的4个节点坐标值，分别对应：x1=-0.387, y1=-1；x2=0.387, y2=-1；x3=0.387, y3=1；x4=-0.387, y4=1。 | ||
表选项
5 结论 1) 采用卷积神经网络实现单元刚度矩阵的计算是可行的，而且卷积神经网络的求解精度可以根据要求设定，但随着问题求解精度的增加，训练网络的时间将会显著增加，在实际应用时应根据具体需要而确定求解的精度。要达到更高的预测精度，网络训练阶段所需时间较长，需要的训练数据也就越多。在卷积神经网络训练完成后，计算新的单元刚度矩阵则是实时的。
2) 在相同的卷积神经网络结构和训练样本数目下，增加卷积神经网络的卷积核数目可以有效提升网络的学习能力。
3) 在相同的卷积神经网络结构下，增加训练样本数目可以有效提升卷积神经网络的学习效果。
4) 对于平面线弹性单元，单元刚度矩阵不会随着单元节点坐标的等比例缩放而改变。
参考文献
| [1] | GOLOVKO V A. Deep learning:An overview and main paradigms[J]. Optical Memory and Neural Networks, 2017, 26(1): 1-17. |
| [2] | 孙道恒, 罗萍平, 胡俏, 等. 有限元单刚矩阵计算的神经网络法[J]. 东北大学学报(自然科学版), 1997, 18(4): 431-434. SUN D H, LUO P P, HU Q, et al. Neural network method for finite element single stiffness matrix calculation[J]. Journal of Northeastern University(Natural Science), 1997, 18(4): 431-434. DOI:10.3321/j.issn:1005-3026.1997.04.005 (in Chinese) |
| [3] | 孙道恒, 胡俏. 弹性力学的实时神经计算原理与数值仿真[J]. 力学学报, 1998, 30(3): 348-353. SUN D H, HU Q. The principle and numerical simulation of real-time neural computing for elastic mechanics[J]. Chinese Journal of Theoretical and Applied Mechanics, 1998, 30(3): 348-353. DOI:10.3321/j.issn:0459-1879.1998.03.012 (in Chinese) |
| [4] | 丁耀武. 线弹性结构静动态有限元法[M]. 沈阳: 东北工学院出版社, 1990: 31-54. DING Y W. Static and dynamic FEM for linear elastic structural analysis[M]. Shenyang: Northeastern University Press, 1990: 31-54. (in Chinese) |
| [5] | 蒋维城. 固体力学有限元分析[M]. 北京: 北京理工大学出版社, 1989: 308-312. JIANG W C. Finite element analysis of solid mechanics[M]. Beijing: Beijing Institute of Technology Press, 1989: 308-312. (in Chinese) |
| [6] | 姜晋庆, 张铎. 结构弹塑性有限元分析法[M]. 北京: 宇航出版社, 1990: 1-10. JIANG J Q, ZHANG D. Structural elastoplastic finite element analysis method[M]. Beijing: Aerospace Publishing House, 1990: 1-10. (in Chinese) |
| [7] | 曾攀. 材料的概率疲劳损伤特性及现代结构分析原理[M]. 北京: 科学技术文献出版社, 1933: 140-141. ZEN P. Probabilistic fatigue damage properties of materials and modern structural analysis principle[M]. Beijing: Science and Technology Literature Press, 1933: 140-141. (in Chinese) |
| [8] | YOSHIOKA T, KARITA S, NAKATANI T.Far-field speech recognition using CNN-DNN-HMM with convolution in time[C]//IEEE International Conference on Acoustics, Speech and Signal Processing.Piscataway, NJ: IEEE Press, 2015: 15361744. |
| [9] | YOSHIOKA T, OHNISHI K, FANG F, et al.Noise robust speech recognition using recent developments in neural networks for computer vision[C]//IEEE International Conference on Acoustics, Speech and Signal Processing.Piscataway, NJ: IEEE Press, 2016: 16021176. |
| [10] | GLOROT X, BENGIO Y.Understanding the difficulty of training deep feedforward neural networks[C]//Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010: 249-256. |
| [11] | ZHANG X, ZOU Y, SHI W.Dilated convolution neural network with LeakyReLU for environmental sound classification[C]//International Conference on Digital Signal Processing.Piscataway, NJ: IEEE Press, 2017: 1-5. |
| [12] | ZHAO P, ZHANG T.Accelerating minibatch stochastic gradient descent using stratified sampling[EB/OL].(2014-05-13)[2019-03-28].https://arxiv.org/abs/1405.3080. |
| [13] | NG A Y.Feature selection, L1 vs.L2 regularization, and rotational invariance[C]//Proceedings of the Twenty-first International Conference on Machine Learning.New York: ACM, 2004: 78. |
| [14] | KINGMA D, BA J.Adam: A method for stochastic optimization[EB/OL].(2014-12-22)[2019-03-28].https://arxiv.org/abs/1412.6980. |
| [15] | 陈国荣. 有限单元法原理及应用[M]. 北京: 科学出版社, 2016: 91-116. CHEN G R. Principle and application of finite element method[M]. Beijing: Science Press, 2016: 91-116. (in Chinese) |
