有限元法是目前计算分析能力最强、应用领域最广的数值计算方法，广泛用于工程结构的不确定性分析，针对有限元法的研究仍是国内外计算力学界的热点.近些年提出的无网格法、广义有限元法、扩展有限元法、多尺度有限元法和光滑有限元法等^[1]都是有限元法的有益补充，其中以光滑有限元取得的进展最为显著.

光滑有限元法是Liu等^[2]将无网格法中的应变光滑措施引入到有限元法中，通过光滑伽辽金弱形式构造的一种计算方法，从而降低了对网格划分的要求，具有高精度和高效率的优点，已广泛应用于工程结构分析中.文献[3]和文献[4]利用edge-based光滑有限元法分别对混合型裂纹问题和压电材料平面问题进行了分析，文献[5]采用node-based光滑有限元法对三维热传递问题进行了分析，文献[6]将光滑有限元法运用到二维声学数值计算中，文献[7]采用Cell-based光滑有限元法对含孔复合材料层合板进行了分析.目前，光滑有限元法所研究的内容主要是确定性问题，其应用范围还有待于进一步拓宽.

对工程结构不确定性——随机性问题的求解，有基于摄动理论的随机有限元法、广义随机有限元法、广义随机光滑有限元法和随机无网格法等^{[8, 9]}，但这些方法不适用于随机参数大变异问题；基于随机场正交展开理论的算法^{[10, 11]}恰恰可以解决该类问题，且具有明显优势.

本文基于随机场正交展开理论，先采用Karhunen-Loève级数将随机场正交分解，再使用混沌多项式将位移随机响应展开，将展开后的随机场及位移响应引入到Cell-based光滑有限元中，导出了基于正交分解理论的Cell-based随机光滑有限元平衡方程；给出了结构位移均值与协方差矩阵的计算公式，求解了材料属性具有随机性的带孔方板的随机响应问题.计算结果与Monte-Carlo模拟5000次的结果进行了对比.

在Cell-based随机光滑有限元法(Cell-based SFEM)中，将定义域Ω离散成N_e个四边形单元，Ω=∪_i=1^N_eΩ_i^e，Ω_i^e∩Ω_j^e=Ø，i≠j，再将Ω_i^e划分成n_e^s∈[1,∞)个四边形光滑域，如图 1所示.

图 1(Figure 1)

Fig. 1

图 1 光滑域的划分及形函数值 Fig. 1 Division of smoothing domain and shape function value

光滑域内任一点x_c处的应变为

式中：ε(x)= ^Δ_sU为FEM中满足位移协调性条件的应变，^Δ_s为微分算子矩阵，U为位移向量；Φ为满足统一性的光滑函数，文中采用下式:

式中：A_c表示第c光滑域面积，.

将式(2)代入式(1)，采用分部积分，得

式中：Γ_c为光滑域Ω_c的边界；n_i,n_j为积分段外法向向量的分量.

将FEM的试函数代入式(3)中，得

式中n_c为光滑域相关节点个数，

式中N_i为形函数.

若光滑单元Ω_c各段边界Γ_ci处采用单高斯点积分，式(5)可改写为

式中：x_b^G和l_b^c分别为光滑边界Γ_cb的高斯点和长度，n_b为每个光滑单元的边界总数.

则光滑单元刚度矩阵为

Cell-based SFEM的控制方程为

式中： ，N 为形函数矩阵，b为体力向量，t为自然边界Γ_t上给定的表面力.

设随机场E(x,θ)为区域Ω上空间坐标点x的函数，则其Karhunen-Loève展开式可描述为

式中：E为随机场E(x,θ)的均值；θ为参数随机性；ξ_i(θ)为互不相关的随机变量；λ_i，ψ_i分别为随机场的协方差函数C_EE(x₁,x₂)的特征值和特征函数；特征函数ψ_i具有正交性，满足下式：

式中δ_ij为Kronecker-delta函数.

λ_i,ψ_i可通过求解第二类Fredholm积分方程得到

解析法求解方程(11)具有一定的局限性，因此利用Galerkin方法求解.

令h_k(x)为Hilbert空间H中完备函数序列，积分核函数C(x₁,x₂)的特征函数可描述为

将式(12)代入式(11)中，截断至n项，将产生误差ε:

由Galerkin法可得，误差ε与Hilbert空间H中基函数h_k(x)正交，即

将式(14)代入式(13)中可得

式(15)可改写为

式中各矩阵元素分别为

将刚度矩阵采用Karhunen-Loève级数展开,得

将式(21)代入式(8)中，得

将位移随机场U采用Karhunen-Loève级数展开,得

式中：χ_j为独立于x的随机变量，g_j,μ_j分别为位移协方差函数C_uu(x₁,x₂)的特征函数和特征值，满足下式：

由于C_uu(x₁,x₂)事先无法知道，不能采用Karhunen-Loève级数分解，因此，可采用混沌多项式展开

式中：a_s为待定系数，s=0,i₁,i₁i₂,i₁i₂i₃,… ；Γ_p为随机变量(ξ_i1,…,ξ_ip)的p阶Hermite多项式，p=0,1,2,….

将混沌多项式(25)截断，取至L项，得

式中：b_l,ψ_l[ξ_r]不同于a_i1i2i3…^k,Γ_p(ξ_i1,…,ξ_ip)，两者之间存在一一对应关系，均匀Hermite混沌多项式ψ_l[ξ_r]形成一组完备正交基，系数b_l可通过Galerkin投影法求得，k维p阶混沌多项式级数的项数L为

表 1给出了混沌多项式级数的项数与阶数和维数的关系.

表 1

Table 1

表 1(Table 1)

表 1 混沌多项式级数的项数与阶数和维数的关系 Table 1 Relationship between the number of terms and the order and dimension of chaotic polynomial series 多项式 维数k 均匀混沌多项式的项数L p=0 p=1 p=2 p=3 p=4 2 1 3 6 10 15 4 1 5 15 35 70 6 1 7 28 83 210

表 1 混沌多项式级数的项数与阶数和维数的关系 Table 1 Relationship between the number of terms and the order and dimension of chaotic polynomial series

将式(26)代入式(23)中，可得

令，可得

将式(29)代入到式(28)中，得

将式(30)两端同乘ψ_m[ξ_r]，m=0,1,2,…,L-1，得

对式(31)取数学期望，得

式中，ψ_l[ξ_r]满足

令c_ilm=ξ_iψ_lξ_rψ_mξ_r，式(32)为

式(34)改写为

由式(35)得

以具有中心圆孔的方板为例，方板在x方向受单向拉伸载荷，σ₀=1.0 N/m，几何模型如图 2所示.方板长度和宽度均为2l=10.0 m，圆孔半径r=1.0 m；设方板处于平面应变状态，泊松比ν=0.3，弹性模量E(x)=E₀(1+α(x))为随机变量，均值E₀=1.0×10³ N/m²，α(x)服从均匀高斯型随机场，均值为0，其自协方差函数为

图 2(Figure 2)

Fig. 2

图 2 有中心圆孔的方板 Fig. 2 Square plate with a circular hole

式中：σ为标准偏差；(x，y)和(x′，y′)为空间x点和x′点的坐标；l_x，l_y分别为随机场沿x，y方向上相关长度.

由于该结构具有对称性，只取板的1/4(右上部分)计算，离散为64个四边形单元，如图 3所示.

图 3(Figure 3)

Fig. 3

图 3 离散单元 Fig. 3 Discrete elements

计算时，每个单元进一步划分为4个四边形光滑域，相关长度l_x = l_y =5.0.

分别对Cell-based随机光滑有限元(CSSFEM)方法的计算结果和Monte-Carlo模拟(MCS)5 000次的结果进行比较.当δ=0.1,0.15,0.2,0.25和0.3时，混沌多项式阶数p=1，随机变量E的Karhunen-Loève展开截断项数L=4，计算得到图 2所示点A，B，C，D，E位移的均值μ_c和标准差σ_c，结果见表 2.由表 2可以看出，CSSFEM仅采用64个单元，就与Monte-Carlo模拟5 000次得到的结果基本一致，具有较高的精度，在δ=0.25,0.3时，CSSFEM的计算结果依旧保持着较高的精度；克服了基于摄动理论的算法因随机量大变异时解的精度低或失效问题；可见，CSSFEM既保持了光滑有限元精度高、网格依赖性小的特点，又具有正交展开理论处理大变异问题的优点.

表 2

Table 2

表 2(Table 2)

表 2 位移的均值μc和标准差σc Table 2 Mean μc and standard deviation σc of the displacement σ 方法 参数 A B C D E ux ux uy uy ux uy 0.1 CSSFEM μc/μm 2 664 5 029 -875 -2 170 4 624 -1 910 σc×106 28.745 51.098 10.964 23.023 46.915 20.020 MCS5000 μc/μm 2 726 5 094 -939.9 -2 189 4 660 -1 940 σc×106 29.720 51.507 11.862 23.482 46.740 20.377 0.15 CSSFEM μc/μm 2 664 5 030 -875 -2 170 4 625 -1 910 σc×106 43.124 76.657 16.449 34.540 70.383 30.035 MCS5000 μc/μm 2 722 5 139 -966 -2 179 4 659 -1 930 σc×106 44.505 77.400E 18.133 34.997 71.202 31.501 0.2 CSSFEM μc/μm 2 665 5 031 -875 -2 170 4 626 -1 910 σc×106 57.512 102 21.938 46.063 93.862 40.053 MCS5000 μc/μm 2 722 5 140 -967 -2 181 4 660 -1 931 σc×106 59.327 109 24.180 46.539 94.127 40.675 0.25 CSSFEM μc/μm 2 666 5 032 -876 -2 170 4 627 -1 910 σc×106 71.909 128 27.431 57.595 117 50.078 MCS5000 μc/μm 2 723 5 142 -967 -2 187 4 661 -1 930 σc×106 74.100 134 30.239 58.289 121 50.921 0.3 CSSFEM μc/μm 2 667 5 034 -876 -2 180 4 628 -1 930 σc×106 86.319 153 32.930 69.136 141 60.111 MCS5000 μc/μm 2 724 5 143 -967 -2 184 4 663 -1 930 σc×106 89.000 159 36.301 69.950 143 61.043

表 2 位移的均值μc和标准差σc Table 2 Mean μc and standard deviation σc of the displacement

图 4和图 5对CSSFEM作了收敛性分析，给出当δ=0.1，取p=1，k=4，单元个数N为64，144，256，400，576，900时，点E在y轴方向位移的均值μ_c和标准差σ_c，并与Monte-Carlo模拟5 000次的结果进行了对比，由图可以看出该方法是正确有效的.

图 4(Figure 4)

Fig. 4

图 4 位置E点处y轴方向位移的均值μc Fig. 4 Mean value μc of y-direction displacement of position E

图 5(Figure 5)

Fig. 5

图 5 位置E点处y轴方向位移的标准差σc Fig. 5

本文提出了基于正交展开理论的Cell-based随机光滑有限元方法，求解了材料属性具有随机性的带孔方板的随机响应问题.通过算例分析，验证了Cell-based随机光滑有限元在处理随机量大变异问题时具有较高的精度，克服了基于摄动理论的算法因随机量大变异时解的精度低或失效问题.该方法对网格依赖性小，继承了光滑有限元的优点.通过与Monte-Carlo模拟结果进行对比，验证了方法的正确性与可行性.