1.引言

随机有限单元法（SFEM）是基于概率理论对有限单元法（FEM）的扩展补充，将不确定性分析引入有限元计算中，采用随机事件去描述模型的边界条件、载荷以及材料参数。计算结果不再是一个确定的单一数据，而是符合随机场理论的一组数据的集合，通过概率密度关系能够得到一定置信区间内的计算结果。

 

Karhunen-Loeve展开法（后文简称K-L法）是广泛使用的随机场离散方法之一，最早在结构的随机有限元分析中首次使用K-L法，并推导了协方差函数的特征函数和特征值需满足的条件。K-L法最初应用于信号学，本质是用有限阶数去截断无穷级数，采用低阶特征值与特征函数去近似描述随机场[3]。截断误差会导致计算结果与理论值存在偏差，计算时需依据相关理论合理地确定截断阶数。采用逐点误差的平均值对K-L展开过程中的截断误差进行了分析，发现任意离散方案下，截断误差均随展开阶数的增大而迅速降低，低阶截断阶数下计算精度仍符合分析要求。分析了小波法、镜像阵法等方法的计算效率及精度，发现算法结构对K-L展开的计算效率存在较大影响。地基土体是复杂地质活动的产物，土性参数的分布是典型的随机场问题，谭晓慧[8]等基于K-L法建立了地基土体的随机场模型，结果表明基于随机场的边坡的可靠度分析结果更为合理。

 

K-L法的广泛应用主要有两个难点，第一是求解第二类Fredholm积分方程获得特征函数与特征值[9~11]，第二是内嵌到现有商用有限元软件。ABAQUS是达索公司开发的大型通用有限元软件，提供了多种计算机语言的二次开发接口。本文基于数值积分算法提出了第二类Fredholm积分方程的一种数值解法，并采用Fortran计算机语言编写了USDFLD子程序与ABAQUS的计算内核联立，实现了K-L展开在ABAQUS中的实现。

 

2.Karhunen-Loeve法的数值解法

2.1Karhunen-Loeve法基础理论

相关函数是描述材料参数相关性的函数，是分析区域内部物理参数与空间坐标的随机关系，相关函数与协方差函数的关系式可表达为：

 

 

 

2.2第二类Fredholm积分方程的数值解法

 

 

 

 

图1计算示意图

 

 

 

 

 

 

 

 

3.合理性验证

 

 

 

 

图2特征值n与展开阶数关系

 

3.2特征函数

K-L法展开的本质是依据式（2），采用有限个特征去描述随机场的整体特征，低阶特征值与特征向量的精度保证了总体展开的合理性。提取前4阶的特征函数()nfk的数值解与理论解（见图3）和计算误差（见表1）。由图3可知，()nfk的数值解与理论解的函数型式一致，分布规律相同。由表1可知，()nfk的数值解计算误差较小，最小误差为0.01?，最大误差为5.46?，满足分析精度的要求。

 

 

（a）特征函数数值解

 

 

（b）特征函数解析解

图3特征函数曲线

 

表1特征函数误差表（?）

 

 

 

3.3协方差函数分析

Ghanem在文献[2]中证明了，当式（2）成立时，协方差函数的谱表达式为：

 

 

 

式中，参数含义同前。

依据式（22）计算各阶展开下的协方差函数，并绘制协方差函数4阶展开的理论解和数值解（见图4）；绘制4、8、12和16阶展开时数值解与理论解的误差图（见图5）。由图4可知，协方差函数的数值解与解析解函数型式一致，当tk=时，协方差函数恒大于零，并以此为分界线，分割成两个对称的区域，矩阵型式的协方差函数的正定性得到验证。

由图5可知，数值解与理论解的误差较小，最小误差仅为-0.00071%。误差图像以tk=为分界呈对称分布，且最大误差区域位于tk=的直线上。由式（22）易知，数值解的累计误差随着展开阶数的增大而增加，4阶展开时最大误差为0.015%，当展开阶数增大到16时，最大误差增加到0.076%，但仍在误差要求范围内。

 

 

图4四阶展开协方差函数

 

 

图5协方差函数计算误差（%）

 

3.4相关距离的影响

计算不同相关距离时的K-L法展开结果，对比不同展开阶数下特征值的数值解与理论解（见图6）。由图6可知，一阶展开对应的特征值n及其收敛速度随着相关距离的增大而增加。展开阶数相同时，大的相关距离下的计算精度更高，结果是一致的。

 

 

（a）数值解

 

 

（b）解析解

图6特征值曲线

 

4.ABAQUS二次开发

4.1算法架构

ABAQUS是法国达索公司开发的大型通用有限元软件，提供了Python和Fortran计算机语言的二次开发接口。Python的应用场景是前后处理阶段，达索公司开发了大量的通用函数，极大地提高了用户的编码效率。Fortran的重心则是求解器的二次开发，作为底层语言能在相同的算法架构下保证最高的计算效率。本文采用Python与Fortran混合编码的方式将K-L法内嵌至ABAQUS（具体架构见图7），其中K-L法关键的特征函数与特征值计算采用Python编码实现，并将展开的随机场以空间坐标的方式输出为csv文件。然后将Python输出的csv文件数据输入ABAQUS前处理?？槟冢谰莘治銮蚴且晃?、二维还是三维，添加场变量对应csv文件中的x坐标、y坐标和z坐标数据。再编制USDFLD子程序将各材料积分点的场变量与坐标对应即可实现K-L法在ABAQUS中的内嵌。

 

 

图7算法架构图

 

4.2USDFLD子程序

USDFLD子程序是ABAQUS内置的场变量控制子程序[15]，允许用户编制Fortran代码将材料积分点的材料参数与场变量联系在一起，场变量可以是一个或多个，数学表达为：

 

 

 

4.3特征函数与特征值计算

Python的Numpy库提供了np.linalg.eig（）函数计算矩阵的特征值和特征向量，在生成K矩阵后直接调用函数即可。协方差函数的特征函数包含于K矩阵无穷多特征向量的集合，软件计算得到的特征向量通常以单位向量表示，尚不满足式（4），需依据式（15~19）计算特征函数的长度。

 

 

 

 

4.4材料参数输出

有限元软件通常不支持将材料参数输出到后处理数据库内，用户使用ABAQUS时只能在后处理?？橄允境”淞縁V的云图（见图8）。

 

 

图8场变量

 

 

 

载初期材料积分点的应力和应变，通过弹性矩阵反算出弹性模量E，将其存储到状态变量中即可得到K-L展开后弹性模量E的分布规律，对于K-L随机场展开，同一分析区域的不同材料参数在空间上的分布规律是一致的，弹性模量E分布的合理性即可对其他参数的分布进行验证（图9为依据式（25）计算得到的地基岩土体的弹性模量随即展开?。?。

 

 

图9弹性模量

 

5.结论

本文基于数值分析理论提出了Karhunen-Loeve展开的一种数值算法，编制了相应程序将其内嵌至ABAQUS软件中，并与理论解答进行了对比，得到以下重要结论：

 

（1）采用Python和Fortran77混合编码的方式将Karhunen-Loeve展开内嵌至ABAQUS软件效果较好，扩展了ABAQUS软件的随机有限元分析功能。

 

（2）有限元软件不能直接在后处理?？槭涑霾牧喜问?，能通过ABAQUS内置的GETVRM函数获取各积分点的应力张量与应变张量，进而通过弹性力学理论生成材料参数的随机场云图。

 

（3）通过数值求积算法，将求解第二类Fredholm积分方程转换为求解K矩阵的特征值与特征向量是可行的。

 

（4）当展开阶数增大，特征值迅速衰减，低阶特征值远大于高阶特征值，保证了有限截断下随机场的计算精度。

 

（5）特征值于特征函数的数值解于理论解基本一致，误差满足分析精度要求。

 

（6）协方差函数的数值解与理论解一致，图像呈对称分布，误差随展开阶数增大而增加，各阶展开下的最大误差均位于中线上。