基于Kriging响应面的有限元模型确认方法与流程

本发明属于模型确认技术领域，尤其设计一种基于kriging响应面的有限元模型确认方法。

背景技术：

在现代工程结构的研发和设计过程中以及重大装备或工程设施在服役阶段的安全性和可靠性评估方面，计算仿真正发挥着越来越重要的作用，特别是产品设计阶段以及复杂系统很难进行全系统试验的情形下，就必须依靠计算仿真进行评估，因此仿真模型的精度和置信度至关重要。模型确认在工程界的最初概念由美国能源部提出，主要应用于战略武器存储管理的可靠性评估和决策中，并于1998年由aiaa起草了计算流体动力学仿真模型验证和确认的指南，oberkampf对此进行了系统总结，并综述了机械工程领域仿真模型确认的发展。jung将模型确认的分层体系应用到一般工程产品的开发中。在计算仿真过程中，由于各种不确定性的存在使得计算仿真的预测结果与试验结果之间往往存在很大差异，因此模型修正技术也得到了广泛的应用，重点对不确定性对模型预测的影响和修正方法进行讨论。模型确认在国内也逐渐受到重视，模型确认的概念最早由张令弥引入到国内，郭勤涛等将模型确认应用到具体研究中，王瑞利和邓小刚分别对模型确认在计算机程序和流体动力学中的应用进行讨论。刘信恩对模型确认中的贝叶斯框架进行讨论和简化。然而，总体来看，模型确认的基本流程尚不明确，还处在研究阶段。

对于复杂结构的输入参数和输出相应特征之间的转换关系，可以通过构建代理模型来代替有限元模型，可以在短时间进行大量随机抽样，大大减少计算量。近年来，响应面模型、kriging模型和神经网络等代理模型快速发展。然而各种代理模型有其各自的优缺点。响应面模型易于实现，但是逼近非线性问题能力较差；kriging模型对非线性问题有较高的准确度，但模型的获取和使用难度较大，神经网络模型的精度和泛化能力依赖其网络结构和大量的学习样本，有时存在“过学习”现象。同时，由于复杂模型的不确定性参数众多，如何选取对目标响应相关性最大的参数也是研究的热点。

技术实现要素：

本发明的发明目的是：为了解决以上问题，本发明提出了一种基于kriging响应面的有限元模型确认方法，为复杂系统计算仿真中的模型确认工作提供便利，同时为产品的结构设计和健康评估等提供依据。

本发明的技术方案是：一种基于kriging响应面的有限元模型确认方法，包括以下步骤：

a、对压力容器模型进行子系统划分，建立有限元模型；

b、对步骤a中有限元模型的模型参数进行参数相关性和参数灵敏度分析；

c、选取与最大等效应力相关性最大的参数，采用拉丁超立方实验设计方法构建压力容器模型最大等效应力的kriging响应面；

d、对步骤c中kriging响应面进行随机抽样，采用核密度估计方法对抽样样本概率分布和置信区间进行分析；

e、将单次有限元计算结果与有限元仿真结果进行对比，实现有限元模型确认。

进一步地，所述步骤b对有限元模型的模型参数进行参数相关性和参数灵敏度分析具体为：

将有限元模型的模型参数分为控制参数和几何材料参数，采用实验设计模块计算几何材料参数对最大等效应力和最大变形量的灵敏度，得到模型参数灵敏度柱状图；采用实验设计模块计算几何材料参数对最大等效应力的相关性，得到模型参数相关性矩阵。

进一步地，所述步骤c中压力容器模型最大等效应力的kriging响应面的数学模型具体为：

其中，为kriging响应面在样本点x处的响应预测值，为估计后的基函数系数，f(x)为拟合函数的基函数，r为相关系数矩阵，yd为样本的观测值，f为f(x)的观测值。

进一步地，所述步骤d对kriging响应面进行随机抽样，采用核密度估计方法对抽样样本概率分布和置信区间进行分析具体为：

采用蒙特卡洛及随机抽样方法，选取5000个半径r、厚度th的参数样本组合，对kriging响应面进行抽样得到5000组最大等效应力幅值，利用直方图法对最大等效应力进行分析；采用非参数核密度估计方法计算等效应力最大值的概率密度函数，通过设置不同的置信度水平，得到压力容器模型的最大等效应力的置信区间。

进一步地，所述步骤e将单次有限元计算结果与有限元仿真结果进行对比，实现有限元模型确认具体为：

根据压力容器模型设计低水平和高水平大气压力试验，构建两组kriging响应面，对每个kriging响应面进行抽样得到认证样本，将两组得到的认证样本利用核密度估计方法进行拟合，得到kde曲线；将kde曲线与两种实验条件下的有限元仿真得到的最大等效应力值进行比较，实现有限元模型确认。

本发明的有益效果是：本发明通过对压力容器模型系统进行分层设计，建立参数化有限元模型，并对模型参数进行参数相关性和参数灵敏度分析，选取与目标响应相关性最大的参数构建kriging响应面，对响应面随机抽样实现对目标响应的预测，利用核密度估计方法对抽样样本概率分布和置信区间进行分析，并与单次有限元计算结果进行对比实现模型确认，从而为模型确认提供便利，同时为产品的结构设计和健康评估等提供依据。

附图说明

图1是本发明的基于kriging响应面的有限元模型确认方法的流程示意图。

图2是本发明实施例中压力容器模型确认分层示意图。

图3是本发明实施例中压力容器模型试验设计载荷施加示意图。

图4是本发明实施例中参数灵敏度柱状图。

图5是本发明实施例中参数相关性矩阵图。

图6是本发明实施例中压力容器模型最大等效应力kriging响应面示意图。

图7是本发明实施例中压力容器模型最大等效应力核密度估计示意图。

图8是本发明实施例中认证实验一的kriging响应面示意图。

图9是本发明实施例中认证实验一的kde曲线示意图。

图10是本发明实施例中认证实验一的cdf曲线示意图。

图11是本发明实施例中认证实验二的kriging响应面示意图。

图12是本发明实施例中认证实验二的kde曲线示意图。

图13是本发明实施例中认证实验二的cdf曲线示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

如图1所示，为本发明的基于kriging响应面的有限元模型确认方法的流程示意图。一种基于kriging响应面的有限元模型确认方法，包括以下步骤：

a、对压力容器模型进行子系统划分，建立有限元模型；

b、对步骤a中有限元模型的模型参数进行参数相关性和参数灵敏度分析；

c、选取与最大等效应力相关性最大的参数，采用拉丁超立方实验设计方法构建压力容器模型最大等效应力的kriging响应面；

d、对步骤c中kriging响应面进行随机抽样，采用核密度估计方法对抽样样本概率分布和置信区间进行分析；

e、将单次有限元计算结果与有限元仿真结果进行对比，实现有限元模型确认。

在步骤a中，由于压力容器模型存在固液气耦合问题，不确定参数较多，因此本发明对压力容器模型进行子系统划分，并建立相关的有限元模型；划分后的子系统拥有较少的不确定参数，在分析过程中可以大大减少不确定参数的误差耦合产生的影响。本发明通过对压力容器的实际载荷和约束情况进行分析，从环境复杂度和系统复杂度两个方面进行考虑，对仿真模型进行子系统划分。如图2所示，为本发明实施例中压力容器模型确认分层示意图。

本发明中压力容器为一个密闭容器，由一个圆柱形壳体和两个半球形壳体构成，在圆柱面和球面连接端施加固定约束，压力容器内有一定的大气压强和一定高度的特殊液体。在pro/e软件中建立压力容器的参数化三维实体模型，并通过与ansysworkbench14.0的无缝链接建立有限元分析模型。

本发明根据有限元模型确认流程和分层思想，对压力容器模型进行子系统划分和试验设计。试验设计分为校准试验、确认试验和认证试验，分别施加不同层次的载荷进行试验。如图3所示，为本发明实施例中压力容器模型试验设计载荷施加示意图，其中左为校准试验，中为确认试验，右为认证试验。

(1)校准试验：常压；

(2)确认试验：施加一定压强p载荷；

(3)认证试验：施加一定压强p和一定高度h液体混合载荷。

校准试验是通过试验对子结构模型参数(如材料弹性模量、泊松比、密度、比重等)进行直接或间接测量，对参数的误差进行估计；确认试验通过构建kriging响应面模型，并对目标响应的概率分布和置信区间进行分析；认证试验在确认过的子系统模型基础上，构建认证试验有限元模型，并将有限元模型仿真结果与响应面预测结果进行比较。

在步骤b中，本发明将有限元模型的模型参数分为控制参数和几何材料参数。控制参数包括压强、液体高度、液体比重，几何材料参数包括长、半径、厚度、弹性模量等。

利用ansysworkbench中的实验设计模块(doe)对确认试验中的不确定参数l、r、th、e进行灵敏度分析，将参数l、r、th、e作为输入参数，将最大等效应力和最大变形量作为输出参数。假设设结构设计参数为x～n(μ,σ)的正态分布，μ为各设计参数初值，σ＝μα为变量标准差，α为设计参数的变异系数，五个设计参数的变异系数均取5％。如表1所示为压力容器模型进行doe分析的不确定性参数。

表1压力容器模型不确定性参数取值

利用doe模块产生25组设计点，并且对每个设计点进行求解，得到25组压力容器的最大等效应力和最大总变形量的样本点。如图4所示，为本发明实施例中参数灵敏度柱状图，可以看出半径r和厚度th对最大应力的灵敏度性是最高的，模型的长度l和弹性模量e的灵敏度则相当小，基本可以忽略不计；而材料参数弹性模量e对最大应变量的灵敏度是最大的。

本发明利用workbench的doe模块进行实验设计可以得到详细的输入参数与输出参数的具体相关性数值。在参数相关性分析当中，若输入参数对输出参数的相关度大于80％，则表示相关性较高，反之，若相关度小于10％，则表示相关性较低。如图5所示为本发明实施例中参数相关性矩阵图，如表2所示为参数相关性矩阵表，可以发现输入参数中的压力容器的厚度th与最大等效应力的相关度达到了81.83％，半径r与最大等效应力的相关度为54.93％，其余输入参数对最大等效应力的相关度均小于10％，可以认为压力容器的半径和厚度对最大等效应力的相关性相对较高，压力容器的长度和材料弹性模量对最大等效应力的相关性很低。利用相关性分析可以选取与关心量(最大等效应力)相关性最大的参数(半径r和厚度th)来构建响应面。

表2参数相关性矩阵表

在步骤c中，本发明采用拉丁超立方试验设计方法(latinhypercubedesign,ld)来构建kriging响应面。选取了25组设计样本点，将输入参数的x轴为压力容器半径r，输入参数的y轴为压力容器壳体厚度th，z轴为最大等效应力(maxeqv.stress)。如图6所示为本发明实施例中压力容器模型最大等效应力kriging响应面示意图

本发明的kriging响应面模型以多项式逼近的方式对计算仿真模型进行近似描述，其表达式为：

其中，fj(x)是拟合函数的基函数，βj是基函数的系数，z(x)是用于拟合的偏差函数。kriging插值方法一般认为不同插值点处的拟合偏差量不是相互独立的，并假设偏差函数是一种随机过程z(x)，随机过程的均值为0，方差为σ²，并且协方差不为0。随机两点t和u处的协方差函数定义为：

cov[z(t),z(u)]＝σ²ρ(t,u；θ)

其中，ρ(t,u；θ)为两点间的相关函数；θ为相关函数的参数，该参数用来衡量两个样本点t与u之间的相关性随着两个样本点间距增加的衰减度，此参数越小，则所构建的响应面就越光滑。

在确定相关函数ρ和样本的观测值yd＝[y(x1),y(x2),…y(xn)]^t后，还要依据样本观察值计算随机过程方差σ²、基函数系数β以及相关函数参数θ，将估计后的三个参数记作即为kriging模型的超参数(hyperparameters)。利用最大似然估计法构建观测样本值的联合概率分布函数p(yd)，通过将联合概率分布函数最大化来得到超参数具体数值：

l(β,σ²,θ)＝p(yd)

其中，

f＝[f(x1),f(x2),…f(xn)]^t

r为相关系数矩阵，rij＝ρ[xi,xi；θ],(1≤i,j≤n)。

参数β、θ表达式为

得到：

相关系数θ利用数值方法进行求解，对极大似然函数求极值得到其数值的大小。记r(x)＝[ρ(x,x1),ρ(x,x2),…,ρ(x,xn),]^t，则kriging响应面模型在任意设计样本点x处的响应预测表达式为：

为了对比不同的响应面实验设计方法对响应面的预测精度的影响，在相同样本点数(25组)的条件下，采用另外两种实验设计方法(水平全因子设计和复合中心设计)来构建kriging响应面模型。同时，利用响应面精度检验指标，对所构建的三种响应面精度进行检验。本发明利用5个随机样本点作为检验点，如表3所示为kriging响应面精度检验结果。

表3kriging响应面精度检验结果

对比三种试验设计方法构建的kriging响应面精度检验指标r²、rmse可以发现，拉丁超立方设计方法可以在相同数量的样本点条件下，构建更加精确的kriging响应面模型。

综上分析，本发明利用拉丁超立方实验设计方法来构建压力容器最大等效应力的kriging响应面是合理有效的。

在步骤d中，本发明的kriging响应面以试验设计参数样本点为基础，来对非试验设计参数样本组合下压力容器壳体表面最大等效应力幅值进行预测。结合蒙特卡洛及随机抽样方法对响应面进行抽样，选5000个半径r、厚度th的参数样本组合，对响应面进行抽样得到5000组最大等效应力幅值。首先利用直方图法对最大等效应力进行分析，观察其分布情况。

由于压力容器最大等效应力的分布类型未知，采用非参数核密度估计方法(kerneldensityestimation,kde)计算得到等效应力最大值的概率密度函数，如图7所示为本发明实施例中压力容器模型最大等效应力核密度估计示意图；并通过设置不同的置信度水平，得到压力容器的最大等效应力的置信区间，如表4所示为等效应力幅值在不同置信度下的置信区间，可以发现在不同置信度条件下，等效应力最大值的上下限并未随置信度的下降发生显著变化。

表4等效应力幅值在不同置信度下的置信区间

在步骤e中，本发明在建立的有限元模型中加入实验特例对建立的有限元模型进行认证评估。本发明设计低水平和高水平大气压力的试验来对确认试验进行认证评估，选取压力容器测试数据中的两组数据，分别设计认证实验一：p＝36.725psi，χ＝0.5，h＝13in；认证实验二：p＝52.884psi，χ＝0.7，h＝45in进行试验，依照确认试验的流程得到在两组控制参数下的kriging响应面，对响应面进行蒙特卡洛随机抽样后得到频率累积直方图、核密度估计曲线和累积分布函数曲线。

从确认过的子结构进行实验设计，在两组试验条件下，构建两组kriging响应面，对每个响应面代理模型进行抽样5000次，就得到5000个认证样本，将两组得到的认证样本利用核密度估计方法进行拟合，得到kde曲线。同时，将两种实验条件下的有限元仿真得到的最大等效应力值与kde曲线进行对比。如图8所示为本发明实施例中认证实验一的kriging响应面示意图。如图9所示为本发明实施例中认证实验一的kde曲线示意图。如图10所示为本发明实施例中认证实验一的cdf曲线示意图。如图11所示为本发明实施例中认证实验二的kriging响应面示意图。如图12所示为本发明实施例中认证实验二的kde曲线示意图。如图13所示为本发明实施例中认证实验二的cdf曲线示意图。通过kde曲线和有限元仿真得到的最大等效应力值进行比较可知，响应面预测数据与仿真数据有一定的偏差，但相差不太大。大部分的响应面预测结果都要比有限元仿真结果要小，因此利用kriging响应面代理模型进行预测的结果偏保守，从而实现有限元模型确认。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。