一种基于Neumann 级数的区间有限元方法1)

伍鹏革 倪冰雨 姜潮

(湖南大学机械与运载工程学院,长沙 410006)

引言

不确定性广泛存在于实际工程中,如结构的材料属性、几何尺寸及结构所受载荷等.传统方法主要基于概率模型对上述不确定性参数进行描述,并使用概率方法对其进行分析[1-5].使用概率方法进行不确定性分析通常需要大量的实验样本以获得精确的概率分布信息.然而在许多工程实际问题中,由于测试成本或测试技术等原因往往造成样本数据缺乏,很难获得足够的样本信息以构建上述参数的精确概率分布函数.为此,一系列非概率分析方法及非精确概率分析方法得以发展,包括证据理论[6-8]、P-box 模型[9-10]、凸模型方法[11-12]、模糊集理论[13-14]和区间方法[15-17]等.其中区间方法于20 世纪50 年代末由Moore[17]提出,最早用于处理计算机内浮点运算问题,后来被弓入到工程结构领域的不确定性问题中[18-20].在使用区间方法进行不确定性分析时,参数的不确定性由上下边界表示,不仅有效适用于小样本条件下的不确定性度量,而且易于工程人员理解和方便使用,近年来在理论研究和实际应用方面获得了广泛关注[21-25].

20 世纪90 年代中期,为解决在小样本条件下工程结构问题中的不确定性,将区间分析与有限元方法[26-27] 相结合,发展出了区间有限元方法[28-30].区间有限元是一种考虑区间不确定性的有限元分析方法,其主要目的在于借助有限元分析这一有力数值分析工具,对含区间不确定参数的结构进行响应边界求解.而获得有限元结构响应边界的关键在于区间平衡方程组的求解,该平衡方程组属于一类区间线性方程组,其中的刚度矩阵或载荷向量为区间矩阵或区间向量.区间线性方程组的求解被认为是一类NP-hard 问题[31-33],一般情况下很难获得其解析结果.因此,不少学者发展了区间有限元分析的近似方法以求解实际响应集合的最小区间包络.Rao 和Berke[34]提出了顶点组合方法,在区间参数波动范围较小的情况下能够得到较为精确的解.Neumaier 和Pownuk[35]针对桁架结构这一典型结构力学问题提出了一种关于区间线性系统的迭代解法,能够得到较准确的桁架结构响应包络解.Muhanna 和Mullen[36]提出了element-by-element(EBE)列式技术,在一定程度上避免了响应区间的过保守估计问题.Rao 等[37]和Xiao 等[38]将EBE 技术与迭代解法相结合应用于结构响应静态以及动态分析中.Muscolino 等[39]对具有不确定性参数结构的频率响应提出了比例多项式的近似求解方法.Sofi 和Romeo[40]基于比例多项式建立了适用于含区间变量及随机变量的结构有限元分析方法.Qiu 和Elishakoff[41]、Chen 和Yang[29]、Qiu等[42]应用区间摄动方法求解了含有有界不确定参数的结构响应区间,对于小扰动问题有较好的求解精度.Degrauwe 等[43]将仿射运算弓入至区间变量的表示中,改善了区间运算中的区间扩张问题.

现有大多数区间有限元方法计算得到的响应区间包络往往过为保守,或只对于参数波动范围较小的区间有限元问题能够得到较准确的响应区间,或计算效率较低等.因而,如何高效求解区间有限元中结构响应的最小包络解仍是有待解决的难点问题.本文归纳了一类在工程实际中常见的区间有限元问题,其中刚度矩阵可以表示为一系列独立区间变量线性叠加的形式,称之为可线性分解式区间有限元问题.如弹性模量或横截面积等为区间参数的桁架结构,弹性模量或厚度尺寸等为区间参数的混凝土板连续体结构等,对其进行区间有限元分析时均可处理为可线性分解式区间有限元问题.针对此类可线性分解式区间有限元问题,采用Neumann 级数展开方法对其刚度矩阵的逆矩阵进行表示,可获得结构响应关于区间变量的显式表达式,从而便于高效求解结构响应的上下边界.本文剩余部分内容安排如下:第1 节首先简要介绍区间有限元分析的概念,并归纳出一类工程中广泛存在的可线性分解式区间有限元问题;第2 节针对可线性分解式区间有限元问题,提出了基于Neumann 级数的区间有限元方法; 第3 节通过2 个算例验证了本文方法的有效性;第4 节给出了全文结论.

1 区间有限元问题描述

在区间有限元分析中,通过区间模型等度量参数的有界不确定性,通常适用于样本信息不足以构建参数概率分布但参数上下边界可知的结构不确定性问题中.由于区间不确定性的存在,通常区间有限元平衡方程中的刚度矩阵为一区间矩阵或节点力向量为一区间向量,从而导致其平衡方程组为一区间线性方程组.本节简要介绍区间有限元平衡方程组的构成,并归纳出一类可线性分解式区间有限元问题.

1.1 区间有限元平衡方程

在有限元方法中,首先将结构离散为由多个单元组成的有限单元模型,计算各个单元的单元刚度矩阵,进而将离散后的单元按照单元与单元之间的共同节点相互连接起来,从而组建全局刚度矩阵K

式中,Ne为结构的离散单元总数目;Ke(i) 为第i个单元的单元刚度矩阵;Ti是对应的单元刚度矩阵的变换矩阵,将单元刚度矩阵从单元局部坐标变换到全局坐标.结构的平衡方程则可表示为

式中,K∈RNd×Nd是结构的全局刚度矩阵;u∈RNd×1表示节点位移向量,为待求量;p∈RNd×1是结构的节点等效载荷向量;Nd为整个结构的自由度数目.当给定边界约束条件后,式(2)中节点位移向量u便可求出,进而可获得结构的应力和应变响应.

在上述通过有限元方法求解结构在外力作用下的力学响应信息时,结构参数均假设为确定性量.对于工程中普遍存在的不确定性结构,如材料属性、几何尺寸、载荷等具有不确定性的情况,通过弓入区间变量描述参数的不确定性,并结合有限元分析方法,即可通过构建区间有限元分析列式对含区间不确定性的结构进行响应边界分析,结构不确定性分析的示意图如图1 所示.结构材料属性或几何尺寸等参数的不确定性将导致结构刚度矩阵K的不确定性,而外加载荷的不确定性会弓起载荷向量p的不确定性.设结构中存在不确定性参数α=[α1,α2,···,αk]T,

图1 结构不确定性分析Fig.1 Structural uncertainty analysis

式(4) 所表示的平衡方程组属于一区间线性方程组,其求解结果Γ通常较为复杂,难以解析获得.实际上,通常寻求的是包络真实解Γ尽可能小的区间边界uI.uI为一区间向量,每个分量描述了对应位移响应量的变化范围,而不考虑各分量之间的相关性,因此uI又称为最小的超立方体近似解[44].目前求解上述区间有限元问题的方法主要包含三类,即基于区间运算的方法[32,45],基于级数展开的方法[41,46]和全局优化方法[47-48].

1.2 一类可线性分解式区间有限元问题

本小节归纳了在工程结构中常见的一类区间有限元问题,该类问题的刚度矩阵可分解为一系列独立区间变量的线性叠加形式,称之为可线性分解式区间有限元问题.

本文仅考虑结构材料属性或几何尺寸参数具有不确定性,而载荷确定的情况,此时式(4)退化为

若参数αI与结构刚度矩阵K中的每个元素呈线性关系,且αI可如式(7) 表示为一组独立区间变量线性叠加的形式,则区间刚度矩阵K(αI)也具有关于区间变量的线性叠加表达式

式中,Kj(0,1,2,···,m)为对称的确定性矩阵.对于不同的区间有限元问题,Kj的具体形式有所不同.形如式(9) 的可线性分解式区间刚度矩阵在工程实际问题中较为常见,如弹性模量或横截面积等为区间参数的桁架结构,具有区间弹性模量或区间厚度尺寸等参数的混凝土板连续体结构等,对其进行数值分析时,相应的区间刚度矩阵均具有上述叠加形式.下面以含区间弹性模量的连续体结构为例,通过推导给出其区间刚度矩阵的线性分解形式.

考虑一由线弹性各向同性材料构成的连续体结构,其有限元模型包含Ne个单元,Nn个节点,Nd个自由度,受到节点载荷p的作用.为便于理解,这里仅考虑弹性模量具有不确定性的情况,且由区间变量描述.设区间弹性模量EI可表示为式(6) 的形式,即

式中,E0为弹性模量EI的中值.

根据有限元理论[26],连续体的单元刚度矩阵Ke(i)有如下形式

式中,Ae(i)为第i个单元的单元积分域,t为厚度,B为几何矩阵,D为弹性系数矩阵,与弹性模量E呈线性关系

式中,D0=E0M,Dj=EjM.将式(13)代入式(11)中,此时单元刚度矩阵也为一区间矩阵

由式(16)可发现,上述问题的区间平衡方程中区间刚度矩阵可表示为关于区间变量的线性叠加形式,因此此类问题即属于一类线性可分解式区间有限元问题.

2 基于Neumann 级数的区间有限元分析方法

区间平衡方程组的求解是区间有限元分析的关键.与确定性的有限元求解相比,在区间有限元分析中最大的不同也是难点即区间平衡方程组中区间矩阵的逆运算难以求解.现有区间有限元分析方法在进行结构响应边界分析时,通常采用一系列数值手段避免对上述区间刚度矩阵进行求逆.本节针对在上一节中归纳的可线性分解式区间有限元问题提出基于Neumann 级数的区间有限元方法,利用Neumann级数展开方法对区间刚度矩阵进行求逆运算,从而便于后续求解结构响应的上下边界.

2.1 基于Neumann 级数的区间刚度矩阵求逆

Neumann 级数是一种适用于算子求逆的级数展开方法.设有任意矩阵A,则该矩阵的逆矩阵A-1可表示为关于矩阵A0的逆矩阵及两者差值矩阵ΔA=A-A0的Neumann 级数展开[49],即

根据式(17),对于可线性分解式区间有限元刚度矩阵的逆矩阵,其Neumann 级数展开式可表示为

由式(18) 可以看出,在区间有限元分析中弓入Neumann 级数表示可线性分解式刚度矩阵的逆,可将原非确定性刚度矩阵的逆运算转换为确定性名义矩阵K0的逆运算及矩阵加、减法和乘法运算,避免了区间矩阵的直接逆变换.通过上述展开方法,式(18)给出了刚度矩阵逆矩阵关于区间变量的显式表达式,从而很大程度上方便了后续对结构响应进行不确定性分析与求解.

2.2 位移边界分析

由式(8)可知位移向量可表示为

将式(18)代入式(19)中,位移响应向量可表示为

u0=,则式(20)可表示为

计算第s(s=1,2,...,Nd)个自由度位移响应的上边界和下边界从而可转换为对如下两个优化问题的求解

式中,u0,s表示向量u0的第s个元素,Δus表示向量Δu的第s个元素.通过求解式(22)和式(23)得到结构在每一自由度上的位移响应区间,从而可获得整个结构的位移响应边界.

借助比格斯（John Biggs）提出的SOLO（Structure of the Observed Learning Outcome，观察到的学习结果的结构）评价理论[17]，对职前教师有关数据分析问题的学科认知水平进行划分．以问题1为例加以阐释．

对于常见的不确定参数扰动较小的问题,通常保留Neumann 展式中的低阶部分即可使结果满足精度要求.进行一阶截断时可得

此时,式(25)和式(26)即为整个结构位移响应一阶截断时的上下边界矢量.

2.3 应力边界分析

通过Neumann 级数展开,不确定性结构的应力也可表示为关于区间变量的显式表达式.单元的应力场表达式为

因此单元应力场也可表示为

当位移响应u取前一阶截断时,σe(i)表示为

通过式(31)和式(32)得到每个单元的应力响应区间,整个结构的应力响应边界即可获得.

2.4 收敛条件与误差分析

根据Neumann 展式(17)的收敛条件可知,区间刚度矩阵逆矩阵的Neumann 级数展开式的收敛条件为

式中,ρ(·)表示矩阵的谱半径.

根据矩阵谱半径的性质有

考虑位移响应向量表达式(20) 取前n阶时,位移响应的绝对误差可表示为

位移响应的相对误差为

3 算例分析

本节通过两个算例验证了本文所提基于Neumann 级数展开区间有限元法的可行性和有效性.第1 个算例是横截面积为区间变量的平面四杆桁架结构在节点力作用下的位移响应分析; 第2 个算例是弹性模量为区间变量的汽车驾驶室受到集中载荷的形变分析.

3.1 平面四杆桁架结构

图2 平面四杆桁架[51]Fig.2 Four bar truss in the plane[51]

图3 沿x 轴方向位移响应上下边界Fig.3 Displacement bounds in the direction of x

图4 沿y 轴方向位移响应上下边界Fig.4 Displacement bounds in the direction of y

表1 节点位移上下边界值相对误差Table 1 Relative error of the displacement bounds

3.2 汽车驾驶室

汽车驾驶室结构设计是汽车设计中的重要环节,驾驶室结构的静动态分析有助于评估设计性能并进行结构改进使其满足设计要求.图5 所示是一个简化的汽车驾驶室车身有限元模型[26,52],该模型包括28 个节点和43 个单元,且节点1′-13′与节点1-13具有对应相同的和坐标及相反的坐标.驾驶室有限元模型详细的节点坐标信息由表2 给出.每一根梁的横截面积为A=0.2 in2(1 in=25.4 mm),惯性矩为Iy′=Iz′=0.003 in4,极坐标矩J=0.006 in4,泊松比为μ=0.3.节点11,11′,12 和12′的所有自由度均被约束,仅仅在节点1 处施加集中载荷,大小为Fx=-990.0 lbf 和Fy=-330.0 lbf(1 lbf=4.45N).由于缺乏足够的实验数据,将梁的弹性模量考虑为不确定区间参数,由于该结构为对称结构,假设对称位置的梁的弹性模量为同一区间参数,各单元的弹性模量区间值具体如表3 所示(以psi 为单位,1 psi=6.894 757 kPa),共有24 个相互独立的区间参数.

图5 汽车驾驶室模型[26,52]Fig.5 A model of auto cab[26,52]

表2 汽车驾驶室模型节点坐标(单位:in)Table 2 Node coordinates of model of auto cab(unit:in)

表3 单元弹性模量的中值和半径(单位:psi)Table 3 Median and radius of modulus of elasticity(unit:psi)

利用本文所提基于Neumann 级数展开前一阶截断计算出的各节点在载荷作用下沿、和轴方向位移响应的上下边界值如图6～图8 所示.图9～图11分别表示3 个方向上的节点位移响应与MCS 及GA结果相对比,该算例中区间变量为24 个,MCS 抽取样本点数为.从图中可看出3 种方法得到的结果基本一致,沿轴方向上的节点位移响应较大,尤其是1,2,3 这3 个节点位移响应最大,可达17 in,且第3 节点在轴方向的位移响应也较大.故而在设计阶段应考虑在这些节点单元处做出优化改进以提升安全性能.

图6 沿x 轴方向位移的上下边界Fig.6 Displacement bounds in the direction of x

图7 沿y 轴方向位移的上下边界Fig.7 Displacement bounds in the direction of y

图8 沿z 轴方向位移的上下边界Fig.8 Displacement bounds in the direction of z

图9 沿x 轴方向位移的上下边界值Fig.9 Displacement bounds in the direction x

图10 沿y 轴方向位移的上下边界值Fig.10 Displacement bounds in the direction y

图11 沿z 轴方向位移的上下边界值Fig.11 Displacement bounds in the direction z

4 结论

本文针对可线性分解式区间有限元问题,提出了一种基于Neumann 级数的区间有限元方法.通过对弹性模量不确定情况下的连续体结构进行举例分析,推导得到了其可线性分解式区间有限元的具体形式.对于工程实际中普遍存在的此类可线性分解式区间有限元问题,采用Neumann 级数对区间刚度矩阵的逆矩阵进行表示,可获得结构响应关于区间变量的显式表达式,从而便于后续对结构响应上下边界的求解.根据Neumann 级数的收敛条件给出了区间有限元求解方法收敛的充分不必要条件,并给出了所求结构响应的相对误差.通过两个算例对所本文方法的可行性与有效性进行了验证.结果表明,本文基于Neumann 级数的区间有限元方法在进行结构响应边界分析时具有较好的求解效率和求解精度.

附 录

区间有限元的Monte Carlo 模拟(MCS)方法借助于Monte Carlo模拟方法的思想,通过对区间参数进行大量抽样并进行确定性有限元分析,并根据响应集合的极值确定响应边界.其具体步骤如下:

(1)获取结构区间参数的上下边界;

(2) 在m维超立方体中按均匀分布抽取区间变量样本,共抽取Ns组样本点;

(3)根据所抽取的样本点确定结构的全局刚度矩阵,进行确定性有限元分析求解平衡方程(2)获得结构位移响应;重复该过程,遍历所有样本变量,获得结构位移响应的集合;

(4) 根据结构位移响应的集合判断出各单元节点位移响应的最大值与最小值,从而构成结构位移响应的上边界和下边界.

注:由上述步骤2 中抽取的样本将均匀分布在m维超立方体范围内.实际上,也可以在抽样中使用任何其他分布形式获得样本,唯一的要求是要确保可以获得m维超立方体范围内的任何点.