CN113436690A - 基于有限数据的横观各向同性材料参数随机样本生成方法 - Google Patents

Abstract

本发明公开了一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法，包括：将四阶弹性刚度张量采用Mandel标记法；将横观各向同性弹性张量用随机张量模拟；将随机张量分解为一组基张量和与材料参数相关的随机系数；基于最大熵原理，创建由随机系数构成的随机向量概率模型，计算其概率密度函数；明确各随机系数样本的生成方法；采用基于M‑H的马尔科夫链蒙特卡洛方法及伽马分布理论，计算概率密度函数的超参数，并生成随机系数样本；将随机系数样本反推成材料参数样本。本发明考虑了材料参数表现出的差异性、离散性和不确定性，可以较为全面地评估材料的客观变异性对有限元模型计算结果的影响，并做出相应的统计特性研究。

Description

基于有限数据的横观各向同性材料参数随机样本生成方法

技术领域

本发明属于随机有限元技术领域，更具体地，涉及一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法。

背景技术

经典有限元法对确定性结构分析的有效性已经得到普遍认可，而在此基础之上发展起来的随机有限元法，使用随机变量或随机过程来描述现实世界中的不确定性因素，因而更加贴近于实际工程问题。随着计算机技术的普及和发展，以Monte-Carlo法为基础的各类数值仿真方法一直起着极其重要的作用。通过对随机变量样本使用有限元程序反复计算，再对结果进行统计分析，可以较为全面地考虑影响有限元模型计算结果的客观变异性，并做出相应的统计特性分析。

常规的有限元法作为一种确定性方法，未考虑问题所涉及的系统中各力学参数的随机性对有限元模型计算结果的影响。事实上，工程领域中的各类构件及结构，由于内、外界结构参数等各种因素，往往存在一定程度的不确定性。

例如，某些复合材料层合板，由于受制造加工工艺的影响，其材料性质通常是异质且各向异性的，材料的力学性能在板平面内基本相同，垂直于板平面方向(厚度方向)不同，可采用横观各向同性本构关系对其力学行为建模。而在宏观尺度上，试验样品之间测量的物理参数表现出较强的差异性和离散性。为了考虑有限元模型中各力学参数的离散性和不确定性，可将材料参数采用随机变量通过参数化方法构造概率模型。因此，如何基于有限试验数据的横观各向同性材料力学参数生成随机样本显得尤为重要。

发明内容

针对现有技术的以上缺陷或改进需求，考虑材料参数表现出的差异性和离散性对有限元模型计算结果的影响，基于随机有限元理论，提供一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法。

为实现上述目的，本发明提供了一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法，包括：

(1)在弹性阶段，横观各向同性材料的应力-应变关系由广义胡克定律给出，其四阶弹性刚度张量采用Mandel标记法，以二阶6×6矩阵形式表示；

(2)将横观各向同性弹性张量用随机张量模拟；

(3)将随机张量分解为一组基张量和与材料参数相关的随机系数；

(4)基于最大熵原理，对随机系数构成的随机向量创建概率模型，计算其概率密度函数；

(5)通过对随机系数间的统计相关性研究，明确各随机系数样本的生成方法；

(6)采用基于Metropolis-Hasting算法的马尔科夫链蒙特卡洛方法及伽马分布理论，计算概率密度函数的超参数，并生成随机系数样本；

(7)将生成的随机系数样本反推成材料参数样本。

在一些可选的实施方案中，在步骤(1)中，采用Mandel标记法以二阶6×6矩阵形式表示四阶弹性刚度张量具体如下：

横观各向同性材料在弹性阶段应力-应变关系为：

其中，

为横向剪切模量，

为横向压缩模量，故横观各向同性弹性张量取决于5个独立的力学参数：纵向杨氏模量E^L、横向杨氏模量E^T，纵向泊松比ν^L、横向泊松比ν^T及纵向剪切模量G^L，σ_ij表示应力分量，ε_kl表示应变分量，i,j,k,l＝1,2,3。

在一些可选的实施方案中，步骤(3)包括：将横观各向同性弹性张量用随机张量模拟表示为：

[C]＝c₁[E⁽¹⁾]+c₂[E⁽²⁾]+c₃([E⁽³⁾]+[E⁽⁴⁾])+c₄[E⁽⁵⁾]+c₅[E⁽⁶⁾]，其中，c₁,...,c₅是一组由材料参数表示的系数，c₁＝E^L+4(ν^L)²k^T,c₂＝2k^T,

c₄＝2G^T,c₅＝2G^L，[E⁽¹⁾],...,[E⁽⁶⁾]是一组四阶基张量，采用Mandel标记法，以二阶6×6的矩阵形式表示为：

在一些可选的实施方案中，步骤(4)包括：

确定随机向量概率密度函数为：

其中，(C₁,C₂,C₃)是统计相关的随机变量，其联合概率密度函数为：

C₄,C₅是统计独立的随机变量，且服从伽马分布，其边缘概率密度函数为：

表示(c₁,c₂,c₃)的取值范围，k₁₂₃为归一化系数，拉格朗日乘子λ＝(λ₁,...,λ₅,λ)，R⁺表示正实数集，R表示实数集，1_S123(c₁,c₂,c₃)表示指示函数，当(c₁,c₂,c₃)∈S₁₂₃时取1，否者取0，k₄、k₅表示为使概率密度函数全域积分为1的归一化系数。

在一些可选的实施方案中，步骤(6)包括：

对目标向量f^target＝(m_C,v_C)和评价向量

间的误差函数J(λ)取极小值：

其中，m_C＝E{C}＝c＝(c ₁,...,c ₅)为随机向量C的均值，随机向量C＝(C₁,...,C₅)，v_C＝E{ln(det([C]))}，由已知的试验数据计算获得，[C]＝C₁[E⁽¹⁾]+C₂[E⁽²⁾]+C₃([E⁽³⁾]+[E⁽⁴⁾])+C₄[E⁽⁵⁾]+C₅[E⁽⁶⁾]，E{·}表示数学期望；

为基于Metropolis-Hasting算法马尔科夫链蒙特卡洛方法构造样本生成器，生成随机样本并对其计算m_C及v_C，α∈[0,1]为自由常数，A_λ表示拉格朗日乘子λ＝(λ₁,...,λ₅,λ)的取值范围。

在一些可选的实施方案中，步骤(7)包括：

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1、考虑了试验数据间的离散性和不确定性，通过对材料参数采用随机变量进行概率建模，并对随机变量样本使用有限元程序反复计算，再对结果进行统计分析，可以较为全面地考虑影响有限元模型计算结果的客观变异性，并做出相应的统计特性研究；

2、仅基于有限的试验数据，可生成分布在边界范围内与试验数据吻合良好的随机材料参数样本。

附图说明

图1是本发明实施例提供的一种方法流程示意图；

图2是本发明实施例提供的一种部分随机系数(C₁,C₂,C₃)边界示意图；

图3是本发明实施例提供的一种以生成(C₁,C₂,C₃)样本计算评价向量分量

及

的收敛性示意图；

图4是本发明实施例提供的一种生成统计相关的(C₁,C₂,C₃)样本示意图；

图5是本发明实施例提供的一种生成统计独立的C₄样本示意图；

图6是本发明实施例提供的一种生成统计独立的C₅样本示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

如图1所示是本发明实施例提供的一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法的流程示意图，包括以下步骤：

步骤S1：采用Mandel标记法，将四阶横观各向同性弹性张量以二阶6×6矩阵形式表示；

步骤S2：将横观各向同性弹性张量用随机张量模拟，并分解为一组基张量和与材料参数相关的随机系数；

步骤S3：基于最大熵原理，对随机系数构成的随机向量创建概率模型，计算其概率密度函数；

步骤S4：研究随机系数间的统计相关性，明确各随机系数样本的生成方法；

步骤S5：基于Metropolis-Hasting算法的马尔科夫链蒙特卡洛方法及伽马分布理论，计算概率密度函数的超参数，并生成随机系数样本；

步骤S6：将生成的随机系数样本反推成材料参数样本。

在本发明实施例中，具体过程如下：

横观各向同性材料在弹性阶段的应力-应变关系，由广义胡克定律给出：

σ_ij＝C_ijklε_kl,其中,i,j,k,l＝1,2,3

上式采用了Mandel标记法，将四阶弹性刚度张量C_ijkl(用双中括号标记，如

)以二阶6×6的矩阵(用单中括号标记，如[C])形式表示。其中，

为横向剪切模量，

为横向压缩模量。故横观各向同性弹性张量取决于5个独立的力学参数：纵向杨氏模量E^L、横向杨氏模量E^T，纵向泊松比ν^L、横向泊松比ν^T及纵向剪切模量G^L，σ_ij表示应力分量，ε_kl表示应变分量。

为了研究横观各向同性弹性张量各分量之间的统计相关性，将四阶弹性张量分解为：

其中{c₁,...,c₅}是一组与材料参数{E^L,E^T,G^L,ν^L,ν^T}相关的系数，

为一组四阶基张量定义如下：

其中

[q]＝[I]-[p]为两个二阶对称张量，n表示垂直于正交各向同性面的单位向量，[I]为二阶单位张量，[[I]]＝[I]⊙[I]为四阶单位张量。符号

与⊙分别表示张量积及对称张量积，定义如下：

对矢量α与β间运算，定义有

对二阶张量[A]与[B]间运算，定义有

及

将四阶张量采用Mandel标记法，以二阶6×6的矩阵形式表示为：

[C]＝c₁[E⁽¹⁾]+c₂[E⁽²⁾]+c₃([E⁽³⁾]+[E⁽⁴⁾])+c₄[E⁽⁵⁾]+c₅[E⁽⁶⁾]

式中基矩阵：

及材料系数：

c₁＝E^L+4(ν^L)²k^T,c₂＝2k^T,

c₄＝2G^T,c₅＝2G^L

将与[C]相应的随机刚度矩阵记为[C]，将系数{c₁,...,c₅}采用随机变量{C₁,...,C₅}建模；定义随机向量C＝(C₁,...,C₅)，对应于可以将随机矩阵[C]以(确定性)基矩阵{[E⁽¹⁾],...,[E⁽⁶⁾]}及随机变量{C₁,...,C₅}为基本单位分解：

[C]＝C₁[E⁽¹⁾]+C₂[E⁽²⁾]+C₃([E⁽³⁾]+[E⁽⁴⁾])+C₄[E⁽⁵⁾]+C₅[E⁽⁶⁾]

因此，构建随机刚度矩阵[C]的概率模型将严格等价于对随机向量C＝(C₁,...,C₅)的概率建模。将随机向量C的概率密度函数及概率分布函数(均未知待求)分别记为p_C(c)及P_C(c)＝p_C(c)dc，其取值区间记为S。基于最大熵原理，在一组定义可用信息的约束条件下，建立随机向量C的概率模型。相关约束条件以数学形式定义如下：

其数学期望平均值记为c，为已知量(可由试验数据评估)

E{C}＝c,其中c＝(c ₁,...,c ₅)

随机弹性矩阵[C]及其逆矩阵[C]^-1(柔度矩阵)，具有有限的二阶矩

E{ln(det([C]))}＝v_C,其中|v_C|＜+∞

其概率密度函数p_C(c)满足归一化条件

∫_Sp_C(c)dc＝1

随机向量C的概率密度函数p_C(c)的求解，需引入如下拉格朗日乘子λ＝(λ₁,...,λ₅,λ)，具体过程不在此赘述，给出数学表达式如下：

(C₁,C₂,C₃)是统计相关的三元随机变量，其联合概率密度函数

由下式给出：

其中

表示(c₁,c₂,c₃)的取值范围，其边界表面为圆锥形，如图2所示。k₁₂₃为归一化系数，R⁺表示正实数集，R表示实数集，

表示指示函数，当(c₁,c₂,c₃)∈S₁₂₃时取1，否者取0。

C₄,C₅是统计独立的随机变量，且服从伽马分布，其边缘概率密度函数为：

p_C4(c₄)＝1_R+(c₄)k₄c₄ ^-2λexp(-λ₄c₄)

p_C5(c₅)＝1_R+(c₅)k₅c₅ ^-2λexp(-λ₅c₅)

对参数(λ₁,λ₂,λ₃,λ)的计算及三元随机变量(C₁,C₂,C₃)样本值的生成，可采用基于M-H采样算法的马尔科夫链蒙特卡洛(MCMC)方法；拉格朗日乘子λ₄及λ₅的值可由试验数据c₄,c₅的平均值及λ确定，k₄、k₅表示为使概率密度函数全域积分为1的归一化系数。

引入目标向量f^target＝(m_C,v_C)，其中m_C＝E{C}＝c＝(c ₁,...,c ₅)为随机向量C的均值(由试验数据计算)，v_C＝E{ln(det([C]))}；定义评估向量

为基于M-H算法的MCMC方法生成随机样本，并对其m_C及v_C进行估算(其值随λ＝(λ₁,...,λ₅,λ)变化，故加上标∧以示区分)。拉格朗日乘子λ的数值最优解可通过对矢量f^target及f^estimate间的误差函数J(λ)取极小值实现：

α∈[0,1]为自由常数，建议取值α＝1/2，A_λ表示拉格朗日乘子λ＝(λ₁,...,λ₅,λ)的取值范围。

联合概率密度函数

可改写为：

其中

假设随机变量(C₁,C₂,C₃)样本值的众数(Mode或modal value)与已获得的试验数据均值(c ₁,c ₂,c ₃)一致：

λ₄及λ₅由伽马分布性质给出：

M-H算法需要预先定义提议分布，一种简单合理的选择就是取正态分布。选取均值向量μ∈R³及协方差矩阵[Σ]为三阶实对称正定方阵的三维正态分布作为提议分布。随机向量X＝(X₁,X₂,X₃)的多元正态分布的概率密度函数p_X(x,μ,[Σ])为：

均值向量μ及协方差矩阵[Σ]取值为：

横观各向同性材料参数E^L、E^T、ν^L、ν^T、G^L、G^T、k^T，可通过以下关系式由生成的随机变量{C₁,...,C₅}样本反推：

本发明某一具体应用实例举例如下。应当说明的是，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

根据某批次27个刨花板样品的材性力学试验结果，将直接测量的材料参数{E^L,E^T,G^L,ν^L,ν^T}转换成系数{c₁,...,c₅}，并计算其均值(c ₁,...,c ₅)及v_C：

c ₁＝0.1269,c ₂＝2.2899,c ₃＝0.0187,

c ₄＝1.5104,c ₅＝0.1065,v_C＝-4.9273

对参数λ赋予初始值λ＝λ⁰，对λ⁰的合理猜测可以通过参数研究确定。试算表明，设置λ⁰＝-100产生的误差函数J(λ)的初始值很小。通过MCMC方法M-H算法生成(C₁,C₂,C₃)样本，拉格朗日乘子的最优解计算为λ＝-100.3622，对应于误差函数J(λ)＝6.286×10^-5，并计算其样本均值

及参数

如下：

f^estimate＝(0.1289,2.2748,0.0158,1.5178,0.1069,-4.8893)

图3展示了在计算概率模型超参数时，生成的(C₁,C₂,C₃)样本均值及参数

的收敛性随样本数量的变化，其中，图3中左图表示生成的(C₁,C₂,C₃)样本均值

随样本数量的变化曲线，右图表示参数

随样本数量的变化曲线。

将确定的拉格朗日乘子λ＝(λ₁,...,λ₅,λ)代入概率模型，再次通过MCMC方法M-H算法生成10000个(C₁,C₂,C₃)样本，如图4，其中十字为试验数据，点为生成的样本，在边界范围内二者分布吻合良好。C₄,C₅为服从伽马分布的独立随机变量，确定其概率模型参数后，其样本的生成相对容易得多。图5、图6展示了其概率密度函数

及生成的样本，其中，图5中左图表示随机变量C₄的概率密度函数，右图表示生成的10000个C₄样本，图6中左图表示随机变量C₅的概率密度函数，右图表示生成的10000个C₅样本。

需要指出，根据实施的需要，可将本申请中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于有限试验数据的横观各向同性材料力学参数随机样本生成方法，其特征在于，包括：

(1)在弹性阶段，横观各向同性材料的应力-应变关系由广义胡克定律给出，其四阶弹性刚度张量采用Mandel标记法，以二阶6×6矩阵形式表示；

(2)将横观各向同性弹性张量用随机张量模拟；

(3)将随机张量分解为一组基张量和与材料参数相关的随机系数；

(4)基于最大熵原理，对随机系数构成的随机向量创建概率模型，计算其概率密度函数；

(5)通过对随机系数间的统计相关性研究，明确各随机系数样本的生成方法；

(6)采用基于Metropolis-Hasting算法的马尔科夫链蒙特卡洛方法及伽马分布理论，计算概率密度函数的超参数，并生成随机系数样本；

(7)将生成的随机系数样本反推成材料参数样本。

2.根据权利要求1所述的方法，其特征在于，在步骤(1)中，采用Mandel标记法以二阶6×6矩阵形式表示四阶弹性刚度张量具体如下：

横观各向同性材料在弹性阶段应力-应变关系为：

其中，

为横向剪切模量，

为横向压缩模量，故横观各向同性弹性张量取决于5个独立的力学参数：纵向杨氏模量E^L、横向杨氏模量E^T，纵向泊松比ν^L、横向泊松比ν^T及纵向剪切模量G^L，σ_ij表示应力分量，ε_kl表示应变分量，i，j，k，l＝1，2，3。

3.根据权利要求2所述的方法，其特征在于，步骤(3)包括：将横观各向同性弹性张量用随机张量模拟表示为：

[C]＝c₁[E⁽¹⁾]+c₂[E⁽²⁾]+c₃([E⁽³⁾]+[E⁽⁴⁾])+c₄[E⁽⁵⁾]+c₅[E⁽⁶⁾]，其中，c₁，...，c₅是一组由材料参数表示的系数，c₁＝E^L+4(v^L)²k^T，c₂＝2k^T，

c₄＝2G^T，c₅＝2G^L，[E⁽¹⁾]，...，[E⁽⁶⁾]是一组四阶基张量，采用Mandel标记法，以二阶6×6的矩阵形式表示为：

4.根据权利要求3所述的方法，其特征在于，步骤(4)包括：

确定随机向量概率密度函数为：

其中，(C₁，C₂，C₃)是统计相关的随机变量，其联合概率密度函数为：

C₄，C₅是统计独立的随机变量，且服从伽马分布，其边缘概率密度函数为：

表示(c₁，c₂，c₃)的取值范围，k₁₂₃为归一化系数，拉格朗日乘子λ＝(λ₁，...，λ₅，λ)，R⁺表示正实数集，R表示实数集，

表示指示函数，当(c₁，c₂，c₃)∈S₁₂₃时取1，否者取0，k₄、k₅表示为使概率密度函数全域积分为1的归一化系数。

5.根据权利要求4所述的方法，其特征在于，步骤(6)包括：

对目标向量f^target＝(m_C，v_C)和评价向量

间的误差函数J(λ)取极小值：

其中，m_C＝E{C}＝c＝(c ₁，...，c ₅)为随机向量C的均值，随机向量C＝(C₁，...，C₅)，v_C＝E{ln(det([C]))}，由已知的试验数据计算获得，[C]＝C₁[E⁽¹⁾]+C₂[E⁽²⁾]+C₃([E⁽³⁾]+[E⁽⁴⁾])+C₄[E⁽⁵⁾]+C₅[E⁽⁶⁾]，E{·}表示数学期望；

为基于Metropolis-Hasting算法马尔科夫链蒙特卡洛方法构造样本生成器，生成随机样本并对其计算m_C及v_C，α∈[0，1]为自由常数，A_λ表示拉格朗日乘子λ＝(λ₁，...，λ₅，λ)的取值范围。

6.根据权利要求5所述的方法，其特征在于，步骤(7)包括：

Images