有限体积法的实体结构静力学随机分析

陈卫东，于艳春，张丰超，路胜卓，贾平，王巍，严涵

(哈尔滨工程大学航天与建筑工程学院，黑龙江哈尔滨150001)

自然界的各种事物是非常复杂的，随机性是它们的固有特性，因而概率设计的思想也被引入到结构的设计中，能够对结构的安全提供更加符合客观实际的评估［1］。近年来，关于结构的随机分析研究已成为结构分析领域里的热点问题之一［2-4］。目前研究结构随机分析的方法有很多，常见的有Monte-Carlo法、随机有限元法等。Monte-Carlo法适用范围广，且结果十分接近精确值，但这种方法因计算工作量很大而限制了它的实际应用。随机有限元法在结构静力学随机分析中得到广泛应用，但由于该方法的基本方程需经过变分原理和加权余量法推导，理论分析相对于有限体积法比较复杂。而有限体积法对方程的离散更简便，不需经过变分原理和加权余量法，直接在控制体上对微分方程进行积分，且在求解三维实体结构模型时，有限体积法在计算节点内力时将体积分运算转换为面积分运算，因而提高了计算效率，且推导的离散方程能更好的保持原微分方程的守恒性、各项物理意义明确。因此，提出了将有限体积法与摄动法相结合，解决实体结构的随机静力学分析问题。另外，有限体积法在计算流体力学领域应用广泛，主要用于计算流体流动和传热等问题［5-7］，但在固体力学中的应用尚不多见，只解决了固体力学里的一些线性和非线性问题［8-11］，从而也拓宽了有限体积法在固体力学领域的应用范围。

1 有限体积法的基本方程

有限体积法的网格分为结构化网格和非结构化网格。结构化网格的特点是网格规则，便于构造，如矩形，长方体。但对于复杂的求解域，无法准确模拟边界。非结构化网格的特点是网格不规则，不便于构造，如三角形、多边形和四面体等。但对于复杂的求解域适应性强。

本文的三维实体结构问题采用图1所示的三维实体结构的非结构实体单元网格:连接四面体的重心、各面的重心及边的中点可以组成6个平面，这6个面将四面体网格分为4个部分，而每个部分又属于与其相连的节点所代表的控制体的一部分。

图1 非结构实体单元网格Fig.1 Unstructured entity unit grid

将有限体积法的基本方程在四面体网格节点所代表的控制体上进行积分得

式中:ρ为结构材料密度。若忽略体力，并利用高斯公式将体积分转为面积分，则上式可写为

式中:nj(j=1，2，3)为控制体表面外法线方向单位矢量的分量。

假设四面体网格内的应力是常量，则控制体积边界上的应力也为常量;又因为控制体积的表面是由平面组成的，则其外法线方向单位矢量的分量也为常数，以图1中的节点1来说，与节点1连接的有m个四面体网格，需将连接的m个网格内的积分累加，式(2)右端的最终形式可以表示为

式中:a1、a2、a3是以节点1为控制体积中心，与其连接的第n个四面体网格的系数，其值与四面体网格节点坐标有关。对于其他节点有同样形式的表达式:

2 摄动有限体积法的基本方程

对于有限体积法中的实体模型，控制体内所包含的某一个四面体n来说，可以将上述推导的有限体积法的基本控制方程(3)写成矩阵形式:

其中:

式中:D为弹性矩阵，B为应变矩阵，δ为位移列阵:

式中:Vn为四面体网格的体积;常系数b1、c1、d1如下

对于某一个控制体积(包含nk4个四面体)来说，可将有限体积法的控制方程写成如下的矩阵形式:

对于每一个控制体积都可建立如上式所示的矩阵形式，边界条件可以通过边界上的条件引进，最终可得如下形式的总体有限体积法求解公式:

根据摄动法的基本原理，假设ai为随机变量Xi在均值点的微小摄动量，显然:

因此，

根据二阶摄动法可得［13］

由式(16)～(18)可以得出位移的均值和方差:

由于响应量的高阶矩较难获得，且误差较大，故这里位移的方差仅取第一项。这样，位移的均值和方差能够得到，进而可求出应力的均值和方差。

3 实体结构算例分析

图2为实体结构的计算模型图，材料的密度为7 800 kg/m3，弹性模量为 210 GPa，泊松比为 0.3。

图2 实体结构计算模型图Fig.2 Calculation model of entity structure

载荷的形式为

式中:a和b作为相互独立的基本随机变量。a的分布为正态分布，其均值为 200.0，变异系数为 0.2，b的分布也为正态分布，其均值为50.0，变异系数为0.2。

由于载荷的随机性对结构响应量的统计特征的影响较大，因此本文仅考虑载荷是随机的，将有限体积法与摄动法相结合求解了实体结构静力学随机分析问题，得出了实体结构响应量对基本随机变量的均值和方差等统计特征，并采用Monte-Carlo模拟来验证提出的方法的正确性。显然，Monte-Carlo法的计算精度与抽样次数N有关，即抽样次数越多，误差越小。为了使Monte-Carlo的结果达到一定的精度，根据文献［12］中建议抽样次数应满足一定的条件，即

式中:Pf为结构的失效概率。通过摄动有限体积法计算得到的响应量统计特征预估了实体结构的失效概率约为45%，因此本文选用5 000次抽样进行了Monte-Carlo模拟，并与摄动有限体积法(P-FVM)的计算结果进行了对比。图3为2种方法计算的响应量的均值和方差对比图。由于节点和单元的数目很大，本文选用部分节点和单元的计算结果进行对比说明。

图3 节点和单元的统计特征Fig.3 Statistical characterization of nodes and elements

通过计算结果对比图可以看出，摄动有限体积法的计算结果和Monte-Carlo模拟的结果吻合较好。为了更好地证明本文提出方法的正确性和计算精度，下面将给出具有典型代表的节点和单元的计算结果与Monte-Carlo模拟的结果的相对误差对比，如表1～4。

表1～4中的计算结果显示，响应量的统计特征的最大相对误差仅为0.59%，能够满足允许的误差范围值。因此，可以证明该方法具有较好的计算精度。从表中的结果也可以看出，响应量的方差的相对误差比均值的大，这是由于响应量的高阶矩较难获得，且误差较大，这里只近似的取方差表达式的第一项，从而影响了它的计算精度。

除此之外，在 CPU 为 i5-2300、2.99 GB内存的计算机上计算了实体结构的静力学随机分析问题。并将2种计算方法的计算时间进行了统计，摄动有限体积法计算耗时为 20 s，而 Monte-Carlo法为110 000 s，可见摄动有限体积法有较高的计算效率。

表1 节点y向位移均值的相对误差Table 1 Relative error of mean value of node displacement in y direction

表2 节点y向位移方差的相对误差Table 2 Relative error of variance of node displacement in y direction

表3 单元等效应力均值的相对误差Table 3 Relative error of mean value of element equivalent stress

表4 单元等效应力方差的相对误差Table 4 Relative error of variance of element equivalent stress

4 结论

采用摄动有限体积法研究了实体结构的静力学随机分析问题:

1)2种方法的计算结果吻合较好。且给出了一些典型节点和单元的相对误差值，满足允许的误差范围值，能够证明摄动有限体积法具有较好的计算精度。

2)统计了2种计算方法的计算时间，可以证明摄动有限体积法具有较高的计算效率。

3)证明了摄动有限体积法能够较为准确、有效地计算实体结构静力学的随机分析问题。为实体结构静力学的随机分析及可靠性分析提出了一种新方法，同时也拓宽了有限体积法在固体力学领域的应用范围。

4)摄动有限体积法也存在一些需要改善的问题，如在计算过程中假设网格内的应力是常量，从而限制了一定的计算精度，发展高精度网格是未来急需解决的问题，且响应量的方差的计算精度较低，在以后的研究中有待提高。

［1］陈虬，刘先斌.随机有限元法及其工程应用［M］.成都:西南交通大学出版社，1993:10-34.

CHEN Qiu，LIU Xianbin.Stochastic finite element method and its application in engineering［M］.Chengdu:Southwest Jiaotong University Press，1993:10-34.

［2］安伟光，朱卫兵，严心池.随机有限元法在不确定性分析中的应用［J］.哈尔滨工程大学学报，2002，23(1):132-135.

AN Weiguang，ZHU Weibing，YAN Xinchi.Application of stochastic finite element method to analysis uncertainties［J］.Journal of Harbin Engineering University，2002，23(1):132-135.

［3］朱位秋，任永坚.随机场的局部平均与随机有限元［J］.航空学报，1986，7(6):604-611.

ZHU Weiqiu，REN Yongjian.Local average of random fields and stochastic finite element method［J］.Acta Aeronautica et Astronautica Sinica，1986，7(6):604-611.

［4］朱位秋，任永坚.基于随机场局部平均的随机有限元法［J］.固体力学学报，1988(4):261-271.

ZHU Weiqiu，REN Yongjian.The stochastic finite element method based on Local average of random field［J］.Acta Mechanica Solida Sinica，1988(4):261-271.

［5］TANIGUCHI N，KOBAYASHI T.Finite volume method on the unstructured grid system ［J］.Computers ＆ Fluids，1991，19(3/4):287-295.

［6］DARWISH M S，WHITEMAN J R，BEVIS M J.Numerical modelling of viscoelastic liquids using a finite-volume method［J］.Journal of Non-Newtonian Fluid Mechanics，1992，45(3):311-337.

［7］JENNY P，LEE S H，TCHELEPI H A.Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media［J］.Journal of Computational Physics，2006，217(2):627-641

［8］FRYER Y D，BAILEY C，CROSS M.A control volume procedure for solving the elastic stress-strain equations on an unstructured mesh ［J］.Appl Math Modelling，1991，15:639-645.

［9］WHEEL M A.A finite volume method for analyzing the bending deformation of thick and thin plates［J］.Comput Meth Appl Mech Eng，1997，147:199-208.

［10］杨海天.有限体积法在固体力学中的应用［J］.大连理工大学学报，1995，35(3):296-298.

YANG Haitian.Application of finite volume method in solid mechanics［J］.Journal of Dalian University of Technology，1995，35(3):296-298.

［11］杨海天.有限体积法在二维固体力学问题中的应用［J］.计算结构力学及其应用，1996，13(3):84-92.

YANG Haitian.Application of finite volume method in 2-dimensional problems of solid mechanics［J］.Computational Structural Mechanics and Applications，1996，13(3):84-92.

［12］安伟光，蔡荫林，陈卫东.随机结构系统可靠性分析与优化设计［M］.哈尔滨:哈尔滨工程大学出版社，2005:36-70.

AN Weiguang，CAI Yinlin，CHEN Weidong.Reliability analysis and optimal design of stochastic structural system［M］.Harbin:Harbin Engineering University Press，2005:36-70.