一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法

摘要

本发明公开了一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法，包括：建立非线性随机结构系统状态响应的矩函数方程；采用高斯截断技术，建立状态响应的封闭矩函数方程；建立时域内的多尺度离散区间；采用Runge‑Kutta方法，在多尺度区间内迭代求解状态响应的封闭矩函数方程，最终获得非线性随机结构系统状态响应的统计特征。本发明的多尺度迭代方法可以避免同时求解三个随机微分方程，不仅节省了大量的计算机存储资源，而且极大地提高了计算效率。此外，通过误差控制和初值设置可以自动调节迭代次数，实现求解过程的加速收敛。本发明为大规模非线性随机结构系统状态分析提供了高效的求解方法，具有实际工程应用价值。

说明书

技术领域

本发明主要是适用于大规模非线性随机结构系统状态响应统计特征的高效求解计算，具体涉及一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法。

背景技术

地震，风载以及海上风浪是常见的环境载荷。在长期受到环境载荷作用下，人造结构的动态响应会影响结构的功能和安全性。因此，结构在环境载荷作用下的动态响应分析尤为重要。由于环境载荷通常处理为随机过程，在动态响应分析时需要建立并求解随机微分方程。

Fokker-Planck-Kolmogorov(FPK)方法广泛用于建立非线性随机系统微分方程，该方法建立的方程本质上等价于随机微分方程。然而，在平稳相关激励情况下，只有当外部随机激励可表示为白噪声滤波，FPK方法才能够使用。尽管借助随机平均方法和滤波方法促进了FPK方法的应用，然而该方法仍然存在两方面缺点：一方面，外部随机激励必须为delta相关随机过程或者可表示为白噪声滤波；另一方面，该方法增加了系统的状态数量。在这种情况下，一些学者更加重视统计矩函数方法的研究。例如，Lutes等研究了在delta相关随机激励下非线性系统的动态响应，并提出了一种建立响应矩函数方程和累积量函数方程的相对直接表达式。Di Paola等采用多项式形式的高斯滤波作为外部随机激励，研究了非线性系统的响应，获得了非高斯响应的统计特征。正如上面提到的，许多研究工作集中在非线性动力响应分析。然而，绝大多数理论研究忽略了矩函数方程数值求解效率的改进研究。显然，低效的数值求解算法会导致非线性动力响应分析方法在实际应用时遇到困难，特别是大规模非线性随机结构系统。此外，delta相关随机激励常常被用于建立非线性系统的矩函数方程，但是这种随机过程并不能准确描述环境载荷的特征。实际上，邱志平等已经研究了受环境载荷作用的动力系统并指出采用直接概率法建立状态响应的矩函数方程十分有效，然而，该研究工作仅适用于线性随机结构系统问题。

本发明针对受到环境载荷作用的非线性随机结构系统，采用直接概率法建立状态响应的矩函数方程。建立时域内的多尺度离散区间，在多尺度区间内迭代求解状态响应的封闭矩函数方程，最终获得状态响应的统计特征。本发明为大规模非线性随机结构系统状态分析提供了高效的求解方法，具有实际工程应用价值。

发明内容

本发明的目的在于克服现有技术的不足，提供一种求解大规模非线性随机结构系统状态的多尺度迭代方法，避免同时求解三个随机微分方程，实现求解过程的快速高效计算，为大规模非线性随机结构系统状态分析提供了有效的计算方法，具有现实的工程应用价值。

为实现上述目的，本发明提供的一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法，该方法实现步骤如下：

步骤一、基于非线性随机结构系统状态方程，通过直接概率法建立非线性系统状态响应的矩函数方程；

步骤二、采用高斯截断技术对状态响应矩函数方程中的非封闭项进行截断，建立状态响应的封闭矩函数方程；

步骤三、将整个分析时间段离散为多个粗尺度时间区间，再在每一个粗尺度时间区间内进行细尺度时间离散，最终建立时域内的多尺度离散区间；

步骤四、在粗尺度时间区间[t_i,t_i+1]上，采用Runge-Kutta方法同时求解封闭的状态响应的均值矩函数方程和封闭的状态响应的均方值矩函数方程，得到状态响应均值和状态响应均方值

步骤五、判断状态响应均值和状态响应均方值的收敛性，如果结果收敛，跳转到步骤七；如果结果不收敛，继续执行步骤六；

步骤六、将步骤四中得到的状态响应均值和状态响应均方值代入封闭的状态响应—激励的相关矩函数方程，计算状态响应—激励相关矩函数令k＝k+1，跳转到步骤四；

步骤七、判断时间条件是否满足t≤T，如果满足条件，令i＝i+1且k＝1，跳转到步骤四；如果不满足条件，输出状态响应均值状态响应均方值和状态响应—激励相关矩函数迭代终止。

其中，步骤一中采用直接概率法建立非线性系统状态响应的矩函数方程，包括：状态响应的均值矩函数方程、状态响应的均方值矩函数方程和状态响应—激励相关矩函数方程，具体表示为：

状态响应的均值矩函数方程：

式中，A₁，A₃为系统矩阵，B为输入矩阵，m_x(t)为状态响应的均值函数，m_u(t)为外部随机激励的均值函数，R_xοxοx(t)为状态响应的联合矩函数，符号“ο”为两个向量的Hadamard乘积，表示两个向量各对应元素相乘，乘积结果仍然为维数相同的向量。

状态响应的均方值矩函数方程：

式中，R_xx(t,t)为状态响应的均方值函数，为状态响应—激励的相关矩函数，其余参数含义同上所述。

状态响应—激励相关矩函数方程：

式中，R_ux(t,t)为状态响应—激励的相关矩函数，R_uu(t,t)为外部随机激励的均方值函数，其余参数含义同上所述。

其中，步骤二中采用高斯截断技术建立状态响应的封闭矩函数方程，包括：封闭的状态响应的均值矩函数方程、封闭的状态响应的均方值矩函数方程和封闭的状态响应—激励的相关矩函数方程，具体表示为：

封闭的状态响应的均值矩函数方程：

式中，各个参数含义同上所述。

封闭的状态响应的均方值矩函数方程：

式中，各个参数含义同上所述。

封闭的状态响应—激励的相关矩函数方程：

式中，各个参数含义同上所述。

其中，步骤三中采用等间隔离散方式建立时域内的多尺度离散区间。首先，将整个分析时间段[t₀,T]离散为多个粗尺度时间区间[t_i,t_i+1]，具体表示为：

D_coarse＝{t_i:t₀＝t₀<t₁<…<t_I＝T}

式中，t₀为初始时刻，t_I＝T为终止时刻，t_i,i＝0,1,…,I为粗尺度时间区间的离散时间点，I为粗尺度时间区间个数。

其次，在每一个粗尺度时间区间[t_i,t_i+1]内进行细尺度时间离散，得到细尺度时间区间[t_i,j,t_i,j+1]，具体表示为：

$D_{fine}^{\left(i\right)}=\{t_{i,j}:t_i=t_{i,0}<t_{i,1}<...<t_{i,J_i}=t_{i+1}\}$

式中，t_i＝t_i,0为第i+1个粗尺度时间区间内的初始时刻，为第i+1个粗尺度时间区间内的终止时刻，t_i,j,j＝0,1,…,J_i为细尺度时间区间的离散时间点，J_i为第i+1个粗尺度时间区间内的细尺度时间区间个数。

其中，步骤四中Runge-Kutta方法的时间步长是细尺度时间区间[t_i,j,t_i,j+1]，状态响应的均方值矩函数方程中的响应—激励相关矩函数取值为具体表示为：

式中，各个参数含义同上所述。

其中，步骤五中判断结果收敛性的条件包括状态响应均值收敛条件和状态响应均方值收敛条件，具体表示为：

$\left(\begin{array}{c}\left|\right|m_x^{\left(k\right)}\left(t\right)-m_x^{(k-1)}\left(t\right)\left|\right|<{ϵ}_1\\ \left|\right|R_{xx}^{\left(k\right)}(t,t)-R_{xx}^{(k-1)}(t,t)\left|\right|<{ϵ}_2\end{array}\right)$

式中，ε₁为状态响应均值允许误差，ε₂为状态响应均方值允许误差。

其中，步骤六中以细尺度时间区间[t_i,j,t_i,j+1]为时间步长，采用Runge-Kutta方法单独求解封闭的状态响应—激励的相关矩函数方程，相关矩函数方程中的状态响应均值取值为相关矩函数方程中的状态响应均方值取值为具体表示为：

式中，各个参数含义同上所述。

本发明与现有技术相比的优点在于：本发明提供了一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法，主要优点如下：

(1)本发明采用直接概率法建立非线性随机结构系统状态响应的矩函数方程，不需要考虑FPK方法涉及的复杂计算。系统的状态响应统计特征可以根据初始条件和外部随机激励的统计特征直接获得，对于环境载荷激励同样适用；

(2)本发明提出了求解非线性随机结构系统状态响应的矩函数方程的新算法，将分析时间段离散为粗尺度时间区间和细尺度时间区间，在多尺度区间内迭代求解矩函数方程，避免了同时求解三个随机微分方程，节省了大量的计算机存储资源；

(3)本发明通过误差控制和初值设置可以根据收敛条件自动调节迭代次数，实现求解过程的加速收敛，极大地提高了计算效率，特别适用于求解大规模非线性随机结构系统状态分析问题，具有实际工程应用价值。

附图说明

图1是本发明多尺度迭代方法实现流程图；

图2是本发明实例中输电塔架结构；

图3是本发明实例中输电塔架结构3端垂直速度的均值；其中图3(a)为3端垂直速度均值的时间历程；图3(b)为图3(a)的局部放大图；

图4是本发明实例中输电塔架结构3端垂直位移的均值；其中图4(a)为3端垂直位移均值的时间历程；图4(b)为图4(a)的局部放大图；

图5是本发明实例中多尺度迭代方法和传统方法的求解时间对比。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，本发明提出了一种高效求解大规模非线性随机结构系统状态的多尺度迭代方法，其具体实现步骤是：

(1)基于非线性随机结构系统状态方程，采用直接概率法建立非线性系统状态响应的矩函数方程，包括：状态响应的均值矩函数方程、状态响应的均方值矩函数方程和状态响应—激励相关矩函数方程。

一个n自由度立方非线性随机结构系统的控制方程表示为：

式中，M为系统质量矩阵，C₁，C₃为系统阻尼矩阵，K₁，K₃为系统刚度矩阵，Q(t)为外部随机激励，y(t)，和分别为系统位移向量，速度向量和加速度向量，符号“ο”为两个向量的Hadamard乘积，表示两个向量各对应元素相乘，乘积结果仍然为维数相同的向量。

引入2n维状态向量和2n维外部随机激励向量表示为：

$x\left(t\right)=\left(\begin{array}{c}\stackrel{·}{y}\left(t\right)\\ y\left(t\right)\end{array}\right),u\left(t\right)=\left(\begin{array}{c}0\\ Q\left(t\right)\end{array}\right)---\left(2\right)$

非线性随机结构系统的控制方程(1)可以改写成非线性随机系统状态方程，具体表示为：

式中，A₁，A₃为系统矩阵，B为输入矩阵，x(t)为系统状态向量，u(t)为外部随机激励向量，具体表示为：

$A_1=\left(\begin{array}{cc}-M^{-1}{C}_1& -M^{-1}{K}_1\\ I& 0\end{array}\right),A_3=\left(\begin{array}{cc}-M^{-1}{C}_3& -M^{-1}{K}_3\\ 0& 0\end{array}\right),B=\left(\begin{array}{cc}0& M^{-1}\\ 0& 0\end{array}\right)---\left(4\right)$

对式(3)取数学期望，可以得到：

定义状态响应的数学期望表达式：

$E\left\{\frac{dx\left(t\right)}{dt}\right\}=\frac{dE\left\{x\left(t\right)\right\}}{dt}=\frac{{\mathrm{dm}}_x\left(t\right)}{dt}---\left(6\right)$

式(5)可以改写为状态响应的均值矩函数方程，具体表示为：

式中，

根据随机向量函数均方值的定义，可以得到：

$\frac{d\left[u\left(t\right)x^T\left(t\right)\right]}{dt}=\frac{d\left(u\right(t\left)\right)}{dt}\left(x^T\right(t\left)\right)+\left(u\right(t\left)\right)\frac{d\left(x^T\right(t\left)\right)}{dt}---\left(9\right)$

将式(3)及其转置表达式代入式(9)，可以得到：

对式(10)取数学期望，可以得到状态响应的均方值矩函数方程，具体表示为：

式中，R_xx(t,t)为状态响应的均方值函数，为状态响应—激励的相关矩函数，具体表示为：

类似的，根据状态响应—激励的相关矩的定义，可以得到：

$\frac{d\left[u\left(t\right)x^T\left(t\right)\right]}{dt}=\frac{d\left(u\right(t\left)\right)}{dt}\left(x^T\right(t\left)\right)+\left(u\right(t\left)\right)\frac{d\left(x^T\right(t\left)\right)}{dt}---\left(13\right)$

将式(3)的转置代入式(13)，可以得到：

对式(14)取数学期望，可以得到状态响应—激励的相关矩函数方程，具体表示为：

式中，R_ux(t,t)为状态响应—激励的相关矩函数，R_uu(t,t)为外部随机激励的均方值函数，具体表示为：

$\begin{array}{c}{R}_{ux}(t,t)=E\left\{u\left(t\right)x^T\left(t\right)\right\}\\ R_{uu}(t,t)=E\left\{u\left(t\right)u^T\left(t\right)\right\}\end{array}---\left(16\right)$

(2)采用高斯截断技术建立状态响应的封闭矩函数方程，包括：封闭的状态响应的均值矩函数方程、封闭的状态响应的均方值矩函数方程和封闭的状态响应—激励的相关矩函数方程。

式(7)，式(11)和式(15)中的高阶非封闭项分别为R_xοxοx(t)，R_(xοxοx)x(t,t)和R_u(xοxοx)(t,t)，对高阶非封闭项应用高斯截断技术可以得到：

将式(17)分别代入式(7)，式(11)和式(15)中，可以得到：

封闭的状态响应的均值矩函数方程：

封闭的状态响应的均方值矩函数方程：

封闭的状态响应—激励的相关矩函数方程：

(3)建立时域内的多尺度离散区间。首先，将整个分析时间段[t₀,T]平均离散为等间隔的时间区间[t_i,t_i+1]，称为粗尺度时间区间，具体表示为：

D_coarse＝{t_i:t₀＝t₀<t₁<…<t_I＝T}(21)

式中，t₀为初始时刻，t_I＝T为终止时刻，t_i,i＝0,1,…,I为粗尺度时间区间的离散时间点，I为粗尺度时间区间个数。

其次，再将每一个粗尺度时间区间[t_i,t_i+1]平均离散为等间隔的子区间[t_i,j,t_i,j+1]，称为细尺度时间区间，具体表示为：

$D_{fine}^{\left(i\right)}=\{t_{i,j}:t_i=t_{i,0}<t_{i,1}<...<t_{i,J_i}=t_{i+1}\}---\left(22\right)$

式中，t_i＝t_i,0为第i+1个粗尺度时间区间内的初始时刻，为第i+1个粗尺度时间区间内的终止时刻，t_i,j,j＝0,1,…,J_i为细尺度时间区间的离散时间点，J_i为第i+1个粗尺度时间区间内的细尺度时间区间个数。

(4)在粗尺度时间区间[t_i,t_i+1]上，采用Runge-Kutta方法同时求解封闭的状态响应的均值矩函数方程和封闭的状态响应的均方值矩函数方程。

Runge-Kutta方法的时间步长取为细尺度时间区间长度Δt＝t_i，j+1-t_i，j，在进行第k次迭代计算时，状态响应的均方值矩函数方程中的响应—激励相关矩函数取值为具体表示为：

在粗尺度时间区间[t_i,t_i+1]上进行第1次迭代计算时(即k＝1)，响应—激励相关矩函数取值为该时间区间初始时刻的响应—激励相关矩R_ux(t_i,t_i)，即

(5)判断状态响应均值和状态响应均方值的收敛性，判断收敛性的条件具体表示为：

$\left(\begin{array}{c}\left|\right|m_x^{\left(k\right)}\left(t\right)-m_x^{(k-1)}\left(t\right)\left|\right|<{ϵ}_1\\ \left|\right|R_{xx}^{\left(k\right)}(t,t)-R_{xx}^{(k-1)}(t,t)\left|\right|<{ϵ}_2\end{array}\right)---\left(24\right)$

式中，ε₁为状态响应均值允许误差，ε₂为状态响应均方值允许误差。

如果相邻两次迭代结果满足式(24)的不等式条件，则结果收敛，跳转到步骤(7)；否则结果不收敛，继续执行步骤(6)。

(6)根据步骤(4)中得到的状态响应均值和状态响应均方值在粗尺度时间区间[t_i,t_i+1]上采用Runge-Kutta方法单独求解封闭的状态响应—激励的相关矩函数方程。

Runge-Kutta方法的时间步长取为细尺度时间区间长度△t＝t_i,j+1-t_i,j，在进行第k次迭代计算时，状态响应—激励的相关矩函数方程中的状态响应均值和状态响应均方值分别取值为和具体表示为：

当式(25)求解完成后，令k＝k+1，跳转到步骤(4)，继续下一次迭代求解。

(7)判断时间条件是否满足t≤T，如果满足条件，令i＝i+1且k＝1，跳转到步骤(4)，继续下一个粗尺度时间区间的第一次迭代求解；如果不满足条件，输出状态响应均值状态响应均方值和状态响应—激励相关矩函数迭代终止。

实施例：

1.结构参数及模型介绍

为了更充分的了解该发明的特点及其对工程实际的适用性，本发明以图2所示的输电塔架结构为例说明多尺度迭代方法高效求解大规模非线性随机结构系统状态响应的有效性。图4中的塔架结构主要由管材组成，截面分为两种，四条边上为主要立柱，矩形截面为□120×120×4mm，其余的管截面为Φ50×5mm。材料的弹性模量为2.06×10⁵MPa，泊松比为0.3。塔架在底部四处基座处固定，6处悬臂端部承受向下的集中载荷作用F(t)＝-psin(4πt)，其中p为随机参数，其均值和标准差分别为μ_p＝1960N和σ_p＝100。假设结构的刚度是非线性的，结构系统的非线性控制微分方程表示为：

式中，M和K分别为线性质量矩阵和线性刚度矩阵，F(t)为外部随机激励。

2.多尺度迭代程序及运行环境

根据多尺度迭代方法，结合现有的确定性有限元分析程序，采用Matlab设计了随机有限元分析程序。在该程序中，结构的几何模型和有限元网格模型由确定性有限元分析程序产生。在本实施例中，输电塔架结构的线性质量矩阵和线性刚度矩阵采用MSC/Nastran有限元软件建立，并采用Nastran中的DMAP提取。系统状态响应的统计特征分析在Matlab R2012上执行，计算机的处理器为3.40GHz Intel(R)Core(TM)i7四核处理器。

3.系统状态响应的统计特征

基于非线性随机结构系统状态方程，采用直接概率法可以获得在外部随机激励作用下的结构系统状态响应的统计特征。图3和图4分别给出了悬臂端3处垂直速度均值的时间历程和垂直位移均值的时间历程。为了比较结果的有效性，图中同时给出了Monte Carlo模拟获得的结果。在随机结构分析中，Monte Carlo模拟结果被认为是精确结果。从图3(a)和图4(a)中可以看到直接概率方法获得的状态响应统计特征与Monte Carlo模拟结果具有很好的一致性。从局部放大图可以清晰地看到两种方法结果的微小差异，这种差异来源于高斯截断法引起的截断误差。对于工程问题，这种微小的差异完全可以忽略不计，因此，直接概率法得到的状态响应统计特征具有一定的工程应用价值。

4.多尺度迭代方法与传统方法对比分析

采用多尺度迭代方法可以避免同时求解三个随机微分方程，提高非线性系统状态响应统计特征的求解效率，从而为大规模非线性随机结构系统状态分析提供了高效的计算方法，促进了随机结构分析方法的工程应用。

为了说明多尺度迭代方法的计算效率，本实施例考虑了五种不同规模的塔架结构有限元模型。五种模型的有限元单元数量从488增加到920，相应的自由度数量从1500增加到2844。在本实施例的多尺度迭代方法中，粗尺度时间取值为△t_coarse＝0.05s，细尺度时间取值为△t_fine＝0.0125s。表1列出了不同规模问题的多尺度迭代方法和传统方法的求解时间以及两种方法耗时之比。从表中可以看到，随着求解问题规模的增加，传统方法的计算时间从510.544s到633.114s，而多尺度迭代方法的计算时间从409.021s到487.245s。相比传统方法，多尺度迭代方法花费了较少的计算时间。此外，表1中的多尺度迭代方法与传统方法的耗时之比大约在76％到80％之间。由此可见，多尺度迭代方法节省了大量的计算时间，特别是大规模工程问题，节省的时间是非常显著的。根据表1中的数据，图5给出了多尺度迭代方法和传统方法的求解时间对比的曲线图。

表1

针对输电塔架结构，采用本发明多尺度迭代方法求解了状态响应的封闭矩函数方程，获得了状态响应的统计特征。该方法避免了同时求解三个随机微分方程，实现了求解过程的快速高效计算。以上实例验证了本发明方法高效求解大规模非线性随机结构系统状态响应的可行性和优越性。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制。

本发明未详细阐述部分属于本领域技术人员的公知技术。