CN106096102A - 基于fem‑kl的非平稳随机动态响应分析方法 - Google Patents

Abstract

本发明公开了一种非平稳随机动态响应分析方法，包括以下步骤：1、确定非平稳随机载荷的均值和自协方差矩阵；2、计算自协方差矩阵的特征值和特征向量，并获得特征值的截断数；3、建立结构的有限元模型，并采用商业有限元中的瞬态分析方法，计算以载荷均值和载荷自协方差矩阵的特征向量分别作为载荷下的均值响应函数和特征向量响应函数；4、基于KL展开获得结构的非平稳随机动态响应，包括响应的方差和自协方差函数。本发明克服传统非平稳随机动态响应分析方法只能用于简单结构的局限性，提供了一种适合于复杂结构非平稳随机动态响应分析方法，同时可以显著提高非平稳随机动态响应分析的计算效率。

Description

基于FEM-KL的非平稳随机动态响应分析方法

技术领域

本发明涉及非平稳随机动态响应分析技术领域，具体涉及一种基于FEM-KL的非平稳随机动态载荷分析方法

背景技术

非平稳随机载荷是工程结构中普遍存在的一种载荷，例如：地震载荷、潮汐载荷、飞行器在服役期间由发动机点火和各种飞行姿态调整引起的振动和冲击动力学载荷等。在工程应用中，由于非平稳随机动态响应分析方法的局限性，常将非平稳随机动态载荷简化为平稳随机动态载荷，然而这样的简化方式会对后续的随机动态响应分析带来不可忽视的误差。

非平稳随机动态响应分析是当前随机动态响应分析的一个技术难点。目前对于非平稳随机动态响应分析，通常采用Karhunen-Loeve(KL)展开和Polynomial Chaos(PC)展开等谱随机有限元技术将非平稳随机动态载荷分解为一系列确定性随机变量后，利用蒙特卡罗法进行动态响应分析。然而该方法仅适用于简单结构的非平稳随机动态响应分析，对于大型复杂结构，由于具有较大的自由度，目前缺乏有效率的分析方法。

发明内容

发明目的：针对现有技术中存在的问题，本发明公开了一种非平稳随机动态响应分析方法，该方法适合于复杂结构的非平稳随机动态响应分析，具有重要的工程应用价值。

技术方案：本发明公开了一种非平稳随机动态响应分析方法，包括如下步骤:

(1)确定非平稳随机动态载荷的均值和自协方差矩阵；

(2)计算自协方差矩阵的特征值和特征向量，并获得特征值的截断数；

(3)建立结构的有限元模型，并采用商业有限元软件中的瞬态分析方法，计算以载荷均值和特征向量分别作为载荷下的均值响应函数和特征向量响应函数；

(4)根据均值响应函数和特征向量响应函数以及KL展开获得结构的非平稳随机动态响应，包括响应的方差和自协方差函数。

进一步地，所述步骤(1)中非平稳随机载荷X(t)的均值μ(t)和自协方差矩阵C(t₁,t₂)计算公式为：

μ(t)＝E[X(t)] (1)

C(t₁,t₂)＝E[(X(t₁)-μ(t₁))(X(t₂)-μ(t₂))] (2)

其中t，t₁,t₂均为时间变量，E[·]表示求期望。

进一步地，步骤(2)包括如下步骤：

(201)将时间t分成m个时间段[t_k-1,t_k]；

其中k＝1，2，3，…，m-1，m；m的取值应大于或等于非平稳随机动态载荷X(t)的时间步数；

(202)选择分段常函数h_k(t)作为正交基，其表达式为：

(203)求解第二类Fredholm积分方程，获得自协方差矩阵的特征值和特征向量；其中第二类Fredholm积分方程为：

MΦ＝ΛNΦ (4)

式中，矩阵M的元素为矩阵N中的元素为矩阵Λ的元素为Λ_ij＝δ_ijλ_i，矩阵Φ＝[φ₁(t),φ₂(t),...,φ_i(t),...,φ_m(t)]^T，φ_i(t)为自协方差矩阵C(t₁,t₂)的第i阶特征向量，λ_i是φ_i(t)对应的特征值，t_min和t_max分别为分析时间的上下界，δ_ij为克罗内克函数，定义如式(5)；i,j＝1,2,……,m；

${\mathrm{\δ}}_{ij}=\left\{\begin{array}{ccc}0,& if& i\≠j\\ 1,& if& i=j\end{array}\right.---\left(5\right)$

(204)获得特征值的截断数n，即自大到小的前n个特征值之和大于所有特征值之和的95％时，在第n阶处截断。

进一步地，所述步骤(4)包括如下步骤：

(401)基于KL展开获得结构在p点的响应为：

$x_p\left(t\right)=A\left(t\right)+\underset{i=1}{\overset{n}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{\mathrm{\ξ}}_i{m}_i\left(t\right)---\left(6\right)$

式中A(t)为以载荷均值μ(t)作为载荷下的均值响应函数；m_i(t)为以第i阶特征向量φ_i(t)作为载荷下的特征向量响应函数；所述ξ_i表示一组标准正态的随机变量，具有均值为0、方差为1的性质；

(402)计算获得响应的方差和自协方差函数R_p(t_k-1,t_k)分别如下：

${\mathrm{\σ}}_p^2\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{m}_i^2\left(t\right)---\left(7\right)$

$R_p(t_{k-1},t_k)=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i{m}_i\left(t_{k-1}\right)m_i\left(t_k\right)---\left(8\right)$

有益效果：本发明公开的非平稳随机动态响应分析方法，将FEM(Finite ElementMethod，有限元方法)和KL方法结合，克服了传统非平稳随机动态响应分析方法只能用于简单结构的局限性，适合于复杂结构非平稳随机动态响应的分析，同时可以显著提高非平稳随机动态响应分析的计算效率。

附图说明

图1为本发明方法的逻辑流程框图；

图2为加筋板的有限元模型；

图3为加筋板A处的位移响应标准差；

图4为加筋板A处的位移响应自协方差函数。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

表1加筋板的几何和材料属性

参数	数值
		加筋板长度	1.1684m
加筋板宽度	1.1684m
		加筋板厚度	0.005m
筋条高度	0.0577m
		筋条厚度	0.003m
筋条与加筋板边的距离	0.400m
		弹性模量	73.08×109N/m2
泊松比	0.33
		密度	2.70×103kg/m3
阻尼比	0.02

以一四边简支加筋板为例，施加均值为零，自协方差为贝塞尔形式的随机载荷进行非平稳随机动态响应分析。载荷步数为1000，载荷时长为1s。加筋板的几何尺寸和材料属性如表1所示。

步骤1：确定随机动态载荷X(t)的均值μ(t)和自协方差矩阵C(t₁,t₂)，分别如式(9)和式(10)：

μ(t)＝0 (9)

$C(t_1,t_2)=10J_0\left(\frac{|t_1-t_2|}{0.1}\right)---\left(10\right)$

其中J₀表示0阶第一类贝塞尔函数；

步骤2：计算自协方差矩阵C(t₁,t₂)的特征值和特征向量，并获得特征值和特征向量的截断数n，具体步骤如下：

201、将时间1s分成1000个时间段[t_k-1,t_k]，其中k＝1，2，3，…，999，1000；

202、选择分段常函数h_k(t)作为正交基，其表达式为：

203、求解第二类Fredholm积分方程，获得自协方差矩阵的特征值λ_i和特征向量φ_i(t)；其中第二类Fredholm积分方程为：

MΦ＝ΛNΦ

其中，矩阵M中的元素矩阵N中的元素矩阵Λ中的元素Λ_ij＝δ_ijλ_i，矩阵Φ＝[φ₁(t),φ₂(t),...,φ_i(t),...,φ₁₀₀₀(t)]^T，φ_i(t)为自协方差矩阵C(t₁,t₂)的第i阶特征向量，λ_i是φ_i(t)对应的特征值，δ_ij为克罗内克函数，i,j＝1,2,……,1000。

204、获得特征值的截断数n。在此实施例中前7阶特征值之和为0.98，所以在第7阶处截断。

步骤3：建立加筋板的有限元模型如图2所示，并采用商业有限元软件中的瞬态分析方法，计算以载荷均值μ(t)和特征向量φ_i(t)分别作为载荷下加筋板A处的位移响应函数A(t)和m_i(t)；

步骤4：基于KL展开获得结构的非平稳随机动态响应，包括响应的方差和自协方差函数。

401、基于KL展开获得加筋板在A处的位移响应为：

$x_A\left(t\right)=0+\underset{i=1}{\overset{7}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{\mathrm{\ξ}}_i{m}_i\left(t\right)---\left(11\right)$

其中ξ_i表示一组标准正态的随机变量，具有均值为0、方差为1的性质；

402、计算获得位移响应的方差和自协方差函数R_p(t_k-1,t_k)分别如以下两式所示，将二者表示为图形的形式，分别如图3和图4所示。

${\mathrm{\σ}}_p^2\left(t\right)=\underset{i=1}{\overset{7}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{m}_i^2\left(t\right)$

$R_p(t_{k-1},t_k)=\underset{i=1}{\overset{7}{\Σ}}{\mathrm{\λ}}_i{m}_i\left(t_{k-1}\right)m_i\left(t_k\right)$

Claims (4)

1.一种非平稳随机动态响应分析方法，其特征在于包括以下步骤：

(1)确定非平稳随机动态载荷的均值和自协方差矩阵；

(2)计算自协方差矩阵的特征值和特征向量，并获得特征值的截断数；

(3)建立结构的有限元模型，并采用商业有限元软件中的瞬态分析方法，计算以载荷均值和特征向量分别作为载荷下的均值响应函数和特征向量响应函数；

(4)根据均值响应函数和特征向量响应函数以及KL展开获得结构的非平稳随机动态响应，包括响应的方差和自协方差函数。

2.根据权利要求1所述的非平稳随机动态响应分析方法，其特征在于，所述步骤(1)中非平稳随机载荷X(t)的均值μ(t)和自协方差矩阵C(t₁,t₂)计算公式为：

μ(t)＝E[X(t)]

C(t₁,t₂)＝E[(X(t₁)-μ(t₁))(X(t₂)-μ(t₂))]

其中t，t₁,t₂均为时间变量，E[·]表示求期望。

3.根据权利要求1所述的非平稳随机动态响应分析方法，其特征在于，所述步骤(2)包括以下步骤：

(201)将时间t分成m个时间段[t_k-1,t_k]；

其中k＝1，2，3，…，m-1，m；m的取值应大于或等于非平稳随机动态载荷X(t)的时间步数；

(202)选择分段常函数h_k(t)作为正交基，其表达式为：

(203)求解第二类Fredholm积分方程，获得自协方差矩阵的特征值和特征向量；其中第二类Fredholm积分方程为：

MΦ＝ΛNΦ；式中，矩阵M的元素为矩阵N中的元素为矩阵Λ的元素为Λ_ij＝δ_ijλ_i，矩阵Φ＝[φ₁(t),φ₂(t),...,φ_i(t),...,φ_m(t)]^T，φ_i(t)为自协方差矩阵C(t₁,t₂)的第i阶特征向量，λ_i是φ_i(t)对应的特征值，t_min和t_max分别为分析时间的上下界，δ_ij为克罗内克函数，i,j＝1,2,……,m；

(204)获得特征值的截断数n，即自大到小的前n个特征值之和大于所有特征值之和的95％时，在第n阶处截断。

4.根据权利要求1所述的非平稳随机动态响应分析方法，其特征在于，所述步骤(4)包括以下步骤：

(401)基于KL展开获得结构在p点的响应为：

$x_p\left(t\right)=A\left(t\right)+\underset{i=1}{\overset{n}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{\mathrm{\ξ}}_i{m}_i\left(t\right);$

式中A(t)为以载荷均值μ(t)作为载荷下的均值响应函数；m_i(t)为以第i阶特征向量φ_i(t)作为载荷下的特征向量响应函数；所述ξ_i表示一组标准正态的随机变量，具有均值为0、方差为1的性质；

(402)计算获得响应的方差和自协方差函数R_p(t_k-1,t_k)分别如下：

${\mathrm{\σ}}_p^2\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}\sqrt{{\mathrm{\λ}}_i}{m}_i^2\left(t\right),$

$R_p(t_{k-1},t_k)=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i{m}_i\left(t_{k-1}\right)m_i\left(t_k\right).$