A PSEUDO-SPECTRAL METHOD FOR THE GENERALIZED DENSITY EVOLUTION EQUATION
-
摘要: 概率密度演化方法(probability density evolution equation, PDEM)为非线性随机结构的动力响应分析提供了新的途径. 通过PDEM获得结构响应概率密度函数(probability density function, PDF)的关键步骤是求解广义概率密度演化方程(generalized probability density evolution equation, GDEE). 对于GDEE的求解通常采用有限差分法, 然而, 由于GDEE是初始条件间断的变系数一阶双曲偏微分方程, 通过有限差分法求解GDEE可能会面临网格敏感性问题、数值色散和数值耗散现象. 文章从全局逼近的角度出发, 基于Chebyshev拟谱法为GDEE构造了全局插值格式, 解决了数值色散、数值耗散以及网格敏感性问题. 考虑GDEE的系数在每个时间步长均为常数, 推导了GDEE在每一个时间步长内时域上的序列矩阵指数解. 由于序列矩阵指数解形式上是解析的, 从而很好地克服了数值稳定性问题. 两个数值算例表明, 通过Chebyshev拟谱法结合时域的序列矩阵指数解求解GDEE得到的结果与精确解以及Monte Carlo模拟的结果非常吻合, 且数值耗散和数值色散现象几乎可以忽略. 此外, 拟谱法具有高效的收敛性且序列矩阵指数解不受CFL (Courant-Friedrichs-Lewy)条件的限制, 因此该方法具有良好的数值稳定性和计算效率.
-
关键词:
- 概率密度演化方法 /
- 广义概率密度演化方程 /
- 拟谱方法 /
- 蒙特卡洛模拟
Abstract: Probability density evolution method (PDEM) provides a new approach for the dynamic response analysis of nonlinear stochastic structures. One of the key steps in obtaining the probability density function (PDF) of the dynamic response of nonlinear stochastic structures through PDEM is to solve the generalized probability density evolution equation (GDEE). The finite difference method is usually used to solve GDEE. Unfortunately, since GDEE is a first order hyperbolic partial differential equation with variable coefficients and discontinuous initial conditions, solving GDEE by finite difference method may encounter grid sensitivity problems, numerical dispersion and numerical dissipation. In this work, from the perspective of global approximation, a global interpolation scheme was constructed for GDEE based on the Chebyshev pseudo-spectral method, which overcame the problems of numerical dispersion, numerical dissipation, and grid sensitivity. Considering that the coefficients of GDEE remain constant at each time step, a sequential matrix exponential solution for GDEE in the time domain was derived. Since the sequent matrix exponential solution of GDEE is analytically formulated, this method can overcome numerical stability issues. The numerical results from two examples demonstrate that GDEE can be effectively solved using a combination of the sequential matrix exponential method for time integration and the Chebyshev pseudo-spectral method for spatial discretization. The numerical solution exhibited excellent agreement with both the exact solution and the results from Monte Carlo simulations, with negligible numerical dissipation and numerical dispersion. In addition, due to the high convergence of pseudo-spectral method and the fact that the sequential matrix exponential method is not restricted by CFL (Courant-Friedrichs-Lewy) condition, the proposed method has excellent numerical stability and computational efficiency. -
引 言
随机振动分析的目的是为了评估随机结构或系统在随机激励作用下响应的统计量. 线性随机系统的振动分析, 已经在概率论的框架内开发了多种方法, 包括功率谱法[1]、矩方程法[2]和虚拟激励法[3]等高效精确的方法, 这一领域的理论研究已趋于成熟. 由于叠加原理不能用于非线性系统, 因此用于线性随机振动分析的方法对于非线性随机结构不再有效. 此外, 非线性系统在高斯激励下响应的统计量通常不再是高斯分布, 其响应的前两阶矩很难反映结构响应的全部概率信息. 因此, 非线性随机结构的振动分析需要获得结构响应的高阶矩, 甚至是结构响应的概率密度函数. 当前, 非线性随机结构的振动分析方法主要包括随机摄动法[4-5]、FPK方程法[6-7]、随机平均法[8-9]、随机响应面法[10-11]和等效线性化法[12-13]等. 然而, 这些方法因其适用范围有限或求解困难而无法应用于工程实际. 此外, 尽管蒙特卡罗模拟[14](Monte Carlo simulations, MCS)被视为非线性随机振动分析的通用解决方案, 但其高昂的计算成本往往限制其实际应用.
近年来, 李杰等[15-18]提出了随机动力系统分析的概率密度演化方法(probability density evolution method, PDEM), 并基于概率守恒原理推导了广义密度演化方程(generalized density evolution equation, GDEE). 该方法为复杂的非线性随机振动分析提供了一种新的解决方案并且已经成功应用于多个领域. 李杰等[19]基于广义密度演化方程和Pontryagin极大值原理提出了一种随机最优控制方法. 贾东林等[20]为评估随机腐蚀作用下在役钢桥的疲劳抗力及其演化特性, 在广义概率密度演化方程的基础上发展了腐蚀疲劳抗力概率密度演化方程. Xu[21]从概率守恒原理角度出发, 导出了疲劳损伤与随机参数联合概率密度函数的演化方程. 并通过恒幅加载疲劳试验和变幅加载疲劳实验结果验证了其有效性和准确性. Xiao等[22]在考虑轨道结构随机影响的情况下, 提出了一种基于概率密度演化法的时变车辆-轨道-桥梁系统随机动力学分析方法. 并以实际建设的车辆-轨道-桥梁系统为例, 验证了所提出的随机模型和评价方法的可行性和可靠性. Peng等[23]将概率密度演化方法与基于方差的灵敏度分析方法相结合, 开发了基础隔震系统结构优化设计方法. Song等[24]采用概率密度演化方法对海上浮式风机在风浪联合作用下的短期动力稳定性进行了评估并对其可靠性进行了评价.
求解GDEE是通过概率密度演化方法获得结构响应响应概率密度函数(probability density function, PDF)的关键之一. GDEE本质上是一个具有时变系数且初始条件间断的一阶双曲偏微分方程. 虽然某些特殊类型的随机振动系统存在GDEE的解析解[25], 但由于实际工程的随机性和复杂性, 通常采用数值方法来求解GDEE. 为了求解GDEE, Li等[18]在有限差分的框架下为GDEE构造了单边差分、L-W双边差分和总变差递减(TVD) 3种数值格式, 计算表明TVD格式能有效地抑制由非光滑初始条件引起的数值耗散和数值色散. Papadopoulos等[26]通过有限元法对GDEE进行了求解, 并通过静力随机有限元问题进行了验证, 结果表明该方法不受CFL条件(Courant-Friedrichs-Lewy 条件, 是双曲型方程显式差分格式的收敛性和稳定性条件)的限制, 提高了计算效率. 然而, 该方法仅对平缓的演变速度进行了验证, 在实际工程中结构在地震、爆炸和风等载荷作用下响应的演变速度通常是剧烈的. 最近, Tao等[27]改进了原来的有限差分方法, 并利用小波函数的多分辨率特性改善了其网格敏感性问题. 上述方法全是从局部插值的角度对GDEE进行求解. 事实上, 采用差分法求解GDEE网格的敏感性问题是不可避免的, 网格过疏, 会导致较大的数值耗散. 而过密的网格, 又会引起剧烈的数值振荡. 此外, 有限差分法很难构造高阶格式, 也会导致结果精度不高.
鉴于此, 本文从全局插值的角度出发, 基于拟谱法对GDEE进行求解. 首先, 通过半离散格式将GDEE转化为一系列的常微分方程, 在合理的假定下推导了GDEE时域的序列矩阵指数解. 然后, 通过Chebyshev拟谱法对GDEE进行空间离散. 最后, 结合时域的序列矩阵指数解和拟谱空间离散格式对GDEE进行求解, 并通过两个数值算例对所提出方法的计算精度和有效性进行了验证.
1. 概率密度演化理论
一般随机动力系统的控制方程和初始条件的形式为[15]
$$ {\dot {\boldsymbol{X}} = {\boldsymbol{G}}\left( {{\boldsymbol{X}},{\boldsymbol{\varTheta}} ,{{t}}} \right),}\qquad{{\boldsymbol{X}}\left( {{{{t}}_0}} \right) = {{\boldsymbol{X}}_0}} $$ (1) 式中, $ {\boldsymbol{X}} = {\left( {{X_1},{X_2}, \cdots ,{X_{{n}}}} \right)^{\mathrm{T}}} $为系统响应的n维向量, ${\boldsymbol{\varTheta}} = {\left( {{\theta _1},{\theta _2}, \cdots ,{\theta _{{n}}}} \right)^{\mathrm{T}}}$为随机参数的n维向量, t为时间.
假定方程(1)是适定的, 其解的形式能被记为
$$ {\boldsymbol{X}} = {\boldsymbol{H}}\left( {{\boldsymbol{\varTheta}} ,{{t}}} \right) $$ (2) 根据概率守恒原理, 可得到其概率密度演化方程
$$ \frac{{\partial {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,t} \right)}}{{\partial t}} + \sum\limits_{l = 1}^{{n}} {{{\dot H}_l}\left( {\theta ,t} \right)\frac{{\partial {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,t} \right)}}{{\partial {x_l}}}} = 0 $$ (3) 式中, $ {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,t} \right) $为t时刻$ \left( {{\boldsymbol{X}},{\boldsymbol{\varTheta}} } \right) $的联合概率密度函数.
对于高维动力系统, 获得所有变量的联合概率密度函数是困难的. 因此, 一维的概率密度演化方程应用更为广泛, 其形式为
$$ \frac{{\partial {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,t} \right)}}{{\partial t}} + \dot H\left( {\theta ,t} \right)\frac{{\partial {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,t} \right)}}{{\partial x}} = 0 $$ (4) 初始条件为
$$ {p_{{\boldsymbol{X}}{\boldsymbol{\varTheta}} }}\left( {x,\theta ,{t_0}} \right) = \delta \left( {x - {x_0}} \right){p_{\boldsymbol{\varTheta}} }\left( \theta \right) $$ (5) 式中,$ \delta (\cdot) $为狄拉克$\delta $函数.
2. GDEE的数值方法
2.1 GDEE的半离散化
考虑到GDEE是一个具有时变系数的一阶双曲型方程, 通过半离散技术容易将概率密度演化方程转换为
$$ \dot {\boldsymbol{p}}{\mathbf{ = }}{\mathcal{L}_k}{\boldsymbol{p}} $$ (6) 式中, $ {\mathcal{L}_k} $为包括p相对于空间坐标的导数的微分时变算子, $ \dot {\boldsymbol{p}} $为相对于时间的导数.
2.2 GDEE的序列解法
对于一维常系数一阶双曲方程, 控制方程的半离散形式为
$$ \dot{\boldsymbol{p}} = \mathcal{L} \boldsymbol{p},\left.\qquad \boldsymbol{p}\right|_{t = 0} = \boldsymbol{p}_0 $$ (7) 式中, $ \mathcal{L} $为微分时不变算子, 其包括${\boldsymbol{p}}$相对于空间坐标的导数, $ \dot {\boldsymbol{p}} $为相对于时间的导数.
方程(7)的精确解是容易得到的, 其形式为
$$ {\boldsymbol{p}}\left( t \right){\mathbf{ = }}{\boldsymbol{p}}\left( {{t_0}} \right){{\mathrm{e}}^{t\mathcal{L}}} $$ (8) 由于GDEE的系数是时变的, 因此方程(7)的结论不能直接应用于概率密度演化方程. 事实上, 在合理的假设下, 离散时间间隔内的概率密度演化方程可以应用方程(7)的结论来求解.
对GDEE的时间变量进行等间距离散, 使得 ${t_k} \in \left[ {0,T} \right],k = 0,1,2,\cdots,n$, 且$\Delta t = {t_{k + 1}} - {t_k}$, 当$\Delta t$足够小时, 可认为在${t_k}$到${t_{k + 1}}$微段内GDEE的演化速度为常数, 此时, 在每一个微段内可直接应用方程(7)的结论. 并将${t_k}$时间步的解作为${t_{k + 1}}$步的初始条件, 可得GDEE的序列精确解
$$ {\boldsymbol{p}}\left( {{{{t}}_{k + 1}}} \right){\mathbf{ = }}{\boldsymbol{p}}\left( {{{{t}}_k}} \right){{\mathrm{e}}^{\Delta t{\mathcal{L}_k}}} $$ (9) 需要注意的是, 上述推导中假定在每一个时间步内GDEE的系数是时不变的, 因此, 采用序列精确解进行GDEE的时间推进时, 其计算精度和时间步长直接相关, 尤其是对于非平稳随机激励下结构的随机响应分析, 这类问题对应的结构的响应演化速度变化更为剧烈, 与之相应的GDEE的系数项具有更强的非线性, 为保证计算精度和算法稳定性, 应选择较小的时间步长.
2.3 Chebyshev 多项式
在谱方法中Chebyshev多项式常被用于处理非周期性问题. Chebyshev多项式${T_n}\left( x \right)$本质上是奇异Sturm-Liouville问题[28]的特征函数
$$ \frac{{\mathrm{d}}}{{{\mathrm{d}}x}}\left( {\sqrt {1 - {x^2}} \frac{{{\mathrm{d}}{T_n}\left( x \right)}}{{{\mathrm{d}}x}}} \right) + \frac{{{n^2}}}{{\sqrt {1 - {x^2}} }}{T_n}\left( x \right) = 0 $$ (10) 事实上, ${T_n}\left( x \right)$恰好是n次多项式, 可通过幂级数展开为
$$ {T_n}\left( x \right) = \frac{n}{2}\sum\limits_{i = 0}^{\left\lfloor {{n \mathord{\left/ {\vphantom {n 2}} \right. } 2}} \right\rfloor } {{{\left( { - 1} \right)}^n}\frac{{\left( {n - i - 1} \right)!}}{{i!\left( {n - 2i} \right)!}}{{\left( {2x} \right)}^{n - 2i}}} $$ (11) 式中, $\left\lfloor {{n \mathord{\left/ {\vphantom {n 2}} \right. } 2}} \right\rfloor $表示${n \mathord{\left/ {\vphantom {n 2}} \right. } 2}$的整数部分.
通过引入参变量并借助余弦三角函数, Chebyshev多项式有更简单的形式
$$ {T_n}\left( x \right) = \cos \left( {n\arccos x } \right) $$ (12) 根据三角函数恒等式关系
$$ \begin{array}{*{20}{c}} {\cos \left[\left( {n + 1} \right)\theta \right]+ \cos \left[\left( {n - 1} \right)\theta\right] = 2\cos \theta \cos (n\theta) ,}&{n \geqslant 1} \end{array} $$ (13) 前几个Chebyshev 多项式以及更高次的3项递推关系式能容易获得
$$ \left.\begin{split} &{T}_{0}\left(x\right) = 1\text{, }\quad {T}_{1}\left(x\right) = x\\ &{T}_{n + 1}\left(x\right) = 2x{T}_{0}\left(x\right)-{T}_{n-1}\left(x\right)\text{, }\quad n\geqslant 1\end{split}\right\} $$ (14) 并且, 通过归纳法Chebyshev 多项式导数的递推关于也能容易获得
$$ \left.\begin{split} &{T}_{0}\left(x\right) = {{T}^{\prime }_{1}}\left(x\right)\text{, }\quad 2{T}_{1}\left(x\right) = 0.5{{T}^{\prime }_{2}}\left(x\right)\\ &{T}_{n}\left(x\right) = \frac{1}{2\left(n + 1\right)}{{T}}'_{n + 1}\left(x\right)-\frac{1}{2\left(n-1\right)}{{T}}'_{n-1}\left(x\right)\text{, }\quad n\geqslant 1\end{split}\right\} $$ (15) 谱方法中之所以使用 Chebyshev 多项式作为试探函数, 原因之一是它在标准区间[−1,1]上满足以下正交性条件
$$ \int_{ - 1}^1 {{T_n}\left( x \right){T_m}\left( x \right)} \omega \left( x \right){\mathrm{d}}x = \frac{{{d_{{n}}}\text{π} }}{2}{\delta _{nm}} $$ (16) 式中, ${\delta _{nm}}$是Kronecker $\delta $函数, 权函数$ \omega \left( x \right) $和常数${d_n}$分别取为
$$\qquad\qquad\qquad \omega \left( x \right) = \frac{1}{{\sqrt {1 - {x^2}} }} $$ (17) $$\qquad\qquad\qquad {d_n} = \left\{\begin{aligned} & 2,\quad {n = 0} \\ & 1,\quad {n \geqslant 1} \end{aligned} \right. $$ (18) 由于这种正交性, Chebyshev 多项式可以构成多项式空间的一组基函数.
2.4 GDEE的空间离散
拟谱法[29]本质上是一种插值技术, 在插值点处的数值解完全满足方程, 为了求解GDEE选择Chebyshev多项式为基函数来构造一个全局插值的格式
$$ p\left( x \right) = \sum\limits_{i = 0}^N {{{\hat p}_i}} {T_i}\left( x \right) $$ (19) 式中,$ {\hat p_i} $是Chebyshev系数, ${T_n}\left( x \right)$是Chebyshev多项式.
对式(19)两边同时关于${T_n}\left( x \right)$做内积, 并考虑 Chebyshev 多项式的正交性条件式(16), 即可得到
$$ {\hat p_i} = \frac{{{{\left[ {p\left( x \right),{T_i}\left( x \right)} \right]}_\omega }}}{{{{\left[ {{T_i}\left( x \right),{T_i}\left( x \right)} \right]}_\omega }}} $$ (20) 为了获得$ p\left( x \right) $对空间坐标x的导数, 对式(19)两边同时取导数, 可得到
$$ p'\left( x \right) = \sum\limits_{i = 0}^N {{{\hat p}_i}} {T'_i}\left( x \right) $$ (21) 将Chebyshev 系数$ {\hat p_i} $用 Chebyshev-Gauss-Lobatto点$ {x_m} $节点处的函数值表示, 并将其代入式(21), 即可得到
$$ {{x_m} = - \cos \left( {\frac{{m\text{π} }}{N}} \right),}\quad {m = 1,2,\cdots,N} $$ (22) $$\begin{split} & p'\left( {{x_m}} \right) = \sum\limits_{i = 0}^N {\left( {\sum\limits_{j = 0}^N {\frac{2}{{{d_j}\text{π} }}{{T'_j}}\left( {{x_m}} \right){T_j}\left( {{x_i}} \right){\omega _i}} } \right)p\left( {{x_i}} \right)} = \\ &\qquad \sum\limits_{i = 0}^N {{D_{m,i}}p\left( {{x_i}} \right)} \end{split} $$ (23) 式中, $ {D_{m,i}} $为Chebyshev微分矩阵, 其显式表达如下
$$ {D_{m,i}} = \left\{ \begin{split} & {\frac{{{c_i}{{\left( { - 1} \right)}^{i + m}}}}{{{c_m}\left( {{x_i} - {x_m}} \right)}},}\quad {i \ne m} \\ &{ - \frac{{{x_m}}}{{2\left( {1 - x_m^2} \right)}},}\quad {1 \leqslant i = k \leqslant N - 1} \\ & {\frac{{2{N^2} + 1}}{6},}\quad {i = m = 0} \\ & {\frac{{2{N^2} + 1}}{6},}\quad {i = m = N} \end{split} \right. $$ (24) 考虑概率密度演化方程(4), 并结合上述的微分矩阵$ {\boldsymbol{D}} $, 能重写概率密度演化方程为
$$ \dot {\boldsymbol{p}} = - \dot H\left( {{t_k}} \right){\boldsymbol{Dp}} $$ (25) 式中, $ {\boldsymbol{p}} = {\left( {{p_{k,0}}\left( t \right),{p_{k,1}}\left( t \right),{p_{k,2}}\left( t \right),\cdots,{p_{k,N}}\left( t \right)} \right)^{\mathrm{T}}} $,向量中元素的第1个下标$ k $表示时间微段的序列, 第2个下标对应Chebyshev-Gauss-Lobatto插值点.
基于式(23), 易得GDEE对应的微分算子$ {\mathcal{L}_k} $
$$ {\mathcal{L}_k} = - \dot H\left( {{t_k}} \right){\boldsymbol{D}} $$ (26) 2.5 初始条件的离散
需要注意的是, GDEE的初始条件是狄拉克函数, 在数值求解时通常采用高斯函数来近似狄拉克函数, 然而对于高斯函数参数标准差的选取通常是根据经验得到的, 太小会引起数值振荡, 太大会导致精度降低, 对于其取值, 本文建议以相邻节点间最大的间距为标准差, 如此不仅能抑制数值振荡同时能保证较高的数值精度.
2.6 数值算法
步骤1: 采用数论方法对概率空间进行离散, 获得向量${\theta _q},q = 1,2,\cdots,{N_{{\mathrm{sel}}}}$, 并对其赋予相应的概率${\omega _q},q = 1,2,\cdots,{N_{{\mathrm{sel}}}}$;
步骤2: 进行确定性求解获得演化速度$ \dot {\boldsymbol{X}} = {\boldsymbol{H}}\left( {{\boldsymbol{\varTheta}} ,t} \right) $;
步骤3: 采用Chebyshev拟谱法对GDEE进行空间离散, 采用序列精确解进行时间推进, 矩阵指数的求解采用Padé方法[30];
步骤4: 根据下式获得累积概率密度函数
$$ p\left( {{x_j},{t_{k + 1}}} \right) = \sum\limits_{q = 1}^{{N_{{\mathrm{sel}}}}} {{w_q}} {p^{\left( q \right)}}\left( {{x_j},{t_{k + 1}}} \right) $$ (27) 3. 数值算例
3.1 无阻尼的单自由度振动系统
考察一个无阻尼的SDOF线性系统, 其自由振动的运动方程为
$$ \ddot X(t) + {\omega ^2}X = 0 $$ (28) 系统的初始条件为
$$ {X(t)\left| {_{t = 0}} \right. = 0.1,}\quad {\dot X(t)\left| {_{t = 0}} \right. = 0} $$ (29) 系统的圆频率是一个随机变量, 在区间$\left[ {5\text{π} /4,7\text{π} /4} \right]$内服从均匀分布. 位移响应的PDF的解析解在文献[18]中有详细的推导.
采用均匀网格选取100个代表性点, 对运动方程进行求解获得确定性响应. 选取300阶的Chebyshev多项式, 计算时间步长取为0.002 s. 根据前面推导的数值方法对GDEE进行求解. 图1是由Chebyshev拟谱法求解GDEE得到的位移响应概率密度曲面. 图2和图3是该数值方法计算得到的概率密度函数和解析解在0.9和1.1 s处的对比, 由图2和图3可以看出该数值方法计算的结果和解析解几乎完全重合, 可见, 采用Chebyshev拟谱法求解GDEE可获得高精度的解. 值得被强调的是本文提出的数值方法计算得到的结果几乎没有数值耗散、数值色散现象, 以及差分法中的网格敏感性问题. 且序列矩阵指数解形式上是解析的, 采用该方法进行时间积分时不受CFL条件的限制, 因此相比其他数值方法本文提出的方法具有良好的数值稳定性. 此外, 拟谱法具有 “指数”收敛性, 拟谱法达到与有限差分法相同精度时所需的节点数更少, 可有效降低矩阵维数, 增加计算效率.
3.2 多自由度的非线性振动系统
考察一个两跨8层的剪切框架结构在地震激励下位移响应的概率密度演化过程. 该结构如图4所示, 一楼柱截面大小为500 mm × 500 mm, 层高为4 m, 其余楼层柱截面大小为400 mm × 400 mm, 层高为3 m. 外部激励为El Centro E-W地震加速度. 结构的质量参数、刚度参数以及加速度峰值具有随机性, 其均值和变异性见表1. 采用Rayleigh阻尼模型, 假定前两阶模态阻尼比为5%.
表 1 随机参数的概率信息Table 1. Probabilistic information of random variablesRandom variable Distribution Mean value COV M normal 2.6 × 105 kg 0.2 E normal 3.0 × 1010 Pa 0.2 PGA normal 1.2 m/s2 0.2 COV: coefficient of variation 采用连续可微Bouc-Wen-Baber-Noori滞回模型[31]模拟层间恢复力, 即
$$ F\left( {u,z} \right) = \alpha ku + \left( {1 - \alpha } \right)kz $$ (30) $$ \dot z = \frac{{h\left( z \right)}}{\eta }\dot u\left\{ {A - \nu [\beta {\mathrm{sign}}(\dot u){{\left| z \right|}^{n - 1}}z + \gamma {{\left| z \right|}^n}]} \right\} $$ (31) 式中, $ F $为恢复力. $\nu $, $\eta $和$h\left( z \right)$分别为与强度退化、刚度退化以及捏缩相关的参数, 参数$\nu $和$\eta $与滞回能$\varepsilon $相关, 其表达式为
$$\qquad\qquad\qquad \nu \left( \varepsilon \right) = 1 + {\delta _\nu }\varepsilon $$ (32) $$\qquad\qquad\qquad \eta \left( \varepsilon \right) = 1 + {\delta _\eta }\varepsilon $$ (33) $$ \qquad\qquad\qquad \varepsilon \left( t \right) = \int_0^t {z\dot u} {\mathrm{d}}t $$ (34) 捏缩函数$h\left( z \right)$的定义为
$$ h\left( z \right) = 1 - {\zeta _1}{{\mathrm{e}}^{ - \tfrac{{{{\left[ {z{\mathrm{sign}}(\dot u) - q{z_u}} \right]}^2}}}{{\zeta _2^2}}}} $$ (35) 式中, ${z_u}$, ${\zeta _1}$ 和 ${\zeta _2}$表达式分别为
$$\qquad\qquad\quad {z_u} = {\left[ {\frac{A}{{\nu (\beta + \gamma )}}} \right]^{1/n}} $$ (36) $$\qquad\qquad\quad {\zeta _1}(\varepsilon ) = {\zeta _s}\left( {1 - {{\mathrm{e}}^{ - p\varepsilon }}} \right) $$ (37) $$\qquad\qquad\quad {\zeta _2}\left( \varepsilon \right) = \left( {\psi + {\delta _\psi }\varepsilon } \right)\left( {\lambda + {\zeta _1}} \right) $$ (38) Bouc-Wen-Baber-Noori模型的参数见表2, 典型的滞回曲线如图5所示. 采用数论方法选择900个代表点进行确定性求解, 选取一层层间位移作为目标响应, 并运用本文提出的方法对GDEE进行数值求解.
表 2 BWBN模型参数Table 2. Parameters of the BWBN modelParameters $ \alpha $ $ \beta $ $ \gamma $ $ n $ $ \delta_{v} $ $ A $ $ \delta_{i} $ Value 0.02 14 6 1 9 0.5 9 Parameters p $ \psi $ $ \delta_{\psi} $ $ \zeta_{I} $ $ \lambda $ q Value 1 1.5 0.6 0.5 7 0.1 为了进一步验证数值方法的准确性, 将5万次Monte Carlo模拟的结果与本文数值方法计算的PDF对比, 从典型时刻的概率密度函数对比图(图6和图7)可以看出二者非常吻合. 另外从结构的滞回曲线图5可以看出结构已进入非线性阶段. 将线性结构的位移概率密度曲面(图1)和非线结构的位移概率密度曲面(图8)进行对比可以看出, 非线性结构的位移概率分布不再是规则形状且随时间显著变化.
![]()
图 8 一层层间位移的概率密度曲面Figure 8. Probability density surface of the first story interlayer displacement4. 结论
在谱方法的框架下, 对GDEE发展了一种空间离散格式, 并在合理的假定条件下推导了GDEE的时间推进的解析解. 给出了相关数值算法的详细介绍. 以线性结构和非线性结构为例, 进行了随机响应分析并与解析解和Monte Carlo模拟的结果进行了对比. 研究表明本文提出的数值算法能较为准确地求解GDEE同时具有较高的数值稳定性. 此外, 拟谱法能有效抑制数值耗散和数值色散, 且没有网格敏感性问题. 当然, 由于时间推进采用了矩阵指数函数, 因此计算效率依然有必要进一步做出研究.
-
![]()
图 8 一层层间位移的概率密度曲面
Figure 8. Probability density surface of the first story interlayer displacement
表 1 随机参数的概率信息
Table 1 Probabilistic information of random variables
Random variable Distribution Mean value COV M normal 2.6 × 105 kg 0.2 E normal 3.0 × 1010 Pa 0.2 PGA normal 1.2 m/s2 0.2 COV: coefficient of variation 
下载: 导出CSV
表 2 BWBN模型参数
Table 2 Parameters of the BWBN model
Parameters $ \alpha $ $ \beta $ $ \gamma $ $ n $ $ \delta_{v} $ $ A $ $ \delta_{i} $ Value 0.02 14 6 1 9 0.5 9 Parameters p $ \psi $ $ \delta_{\psi} $ $ \zeta_{I} $ $ \lambda $ q Value 1 1.5 0.6 0.5 7 0.1
下载: 导出CSV
-
[1] Priestley MB. Power spectral analysis of non-stationary random processes. Journal of Sound and Vibration, 1967, 6: 86-97
[2] Frangopol DM. Review of random vibration of mechanical and structural systems by T.T. Soong and Mircea Grigoriu. Journal of Engineering Mechanics, 1996, 122: 184-184
[3] Lin JH, Zhang YH, Li QS, et al. Seismic spatial effects for long-span bridges, using the pseudo excitation method. Engineering Structures, 2004, 26: 1207-1216
[4] Crandall SH. Perturbation techniques for random vibration of nonlinear systems. The Journal of the Acoustical Society of America, 2005, 35: 1700-1705
[5] Kamiński M. Uncertainty analysis in solid mechanics with uniform and triangular distributions using stochastic perturbation-based finite element method. Finite Elements in Analysis and Design, 2022, 200: 103648
[6] Li J, Jiang ZM. A data-based CR-FPK method for nonlinear structural dynamic systems. Theoretical and Applied Mechanics Letters, 2018, 8: 231-244
[7] Cui W, Caracoglia L, Zhao L, et al. Examination of occurrence probability of vortex-induced vibration of long-span bridge decks by Fokker-Planck-Kolmogorov equation. Structural Safety, 2023, 105: 102369
[8] Liu YL, Liu LQ, Dostal L, et al. The applicability of stochastic averaging method to solve the ship rolling response excited by narrow-band waves. Ocean Engineering, 2022, 251: 111109
[9] Roberts JB, Spanos PD. Stochastic averaging: An approximate method of solving random vibration problems. International Journal of Non-Linear Mechanics, 1986, 21: 111-134
[10] 左彦飞, 江志农, 冯坤等. 转子支承系统临界转速概率分析的随机响应面法. 振动与冲击, 2019, 38(24): 85-90, 114 (Zuo Yanfei, Jiang Zhinong, Feng Kun, et al. A stochastic response surface method for probability analysis of critical speeds of a rotor systems. Journal of Vibration and Shock, 2019, 38(24): 85-90, 114 (in Chinese) Zuo Yanfei, Jiang Zhinong, Feng Kun, et al. A stochastic response surface method for probability analysis of critical speeds of a rotor systems. Journal of Vibration and Shock, 2019, 38(24): 85-90, 114 (in Chinese)
[11] He JJ, Huang M, Wang W, et al. An asymptotic stochastic response surface approach to reliability assessment under multi-source heterogeneous uncertainties. Reliability Engineering & System Safety, 2021, 215: 107804
[12] Su C, Xian JH, Huang H. An iterative equivalent linearization approach for stochastic sensitivity analysis of hysteretic systems under seismic excitations based on explicit time-domain method. Computers & Structures, 2021, 242: 106396
[13] Rastehkenari SF, Ghadiri M. Nonlinear random vibrations of functionally graded porous nanobeams using equivalent linearization method. Applied Mathematical Modelling, 2021, 89: 1847-1859
[14] Shinozuka M. Monte Carlo solution of structural dynamics. Computers & Structures, 1972, 2: 855-874
[15] 李杰, 陈建兵. 随机结构动力反应分析的概率密度演化方法. 力学学报, 2003, 35(4): 437-442 (Li Jie, Chen Jianbing. Probability density evolution method for analysis of stochastic structural dynamic response. Chinese Journal of Theoretical and Applied Mechanics, 2003, 35(4): 437-442 (in Chinese) Li Jie, Chen Jianbing. Probability density evolution method for analysis of stochastic structural dynamic response. Chinese Journal of Theoretical and Applied Mechanics, 2003, 35(4): 437-442 (in Chinese)
[16] Li J. Probability density evolution method: Background, significance and recent developments. Probabilistic Engineering Mechanics, 2016, 44: 111-117
[17] 陈建兵, 李杰. 非线性随机结构动力可靠度的密度演化方法. 力学学报, 2004, 36(2): 196-201 (Chen Jianbing, Li Jie. The probability density evolution method for dynamic reliability assessment nonlinear stochastic structures. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(2): 196-201 (in Chinese) Chen Jianbing, Li Jie. The probability density evolution method for dynamic reliability assessment nonlinear stochastic structures. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(2): 196-201 (in Chinese)
[18] Li J, Chen JB. Probability density evolution method for dynamic response analysis of structures with uncertain parameters. Computational Mechanics, 2004, 34: 400-409
[19] 李杰, 彭勇波. 基于广义密度演化方程的结构随机最优控制. 计算力学学报, 2010, 27(6): 976-982 (Li Jie, Peng Yongbo. Generalized density evolution equation based structural stochastic optimal control. Chinese Journal of Computational Mechanics, 2010, 27(6): 976-982 (in Chinese) Li Jie, Peng Yongbo. Generalized density evolution equation based structural stochastic optimal control. Chinese Journal of Computational Mechanics, 2010, 27(6): 976-982 (in Chinese)
[20] 贾东林, 张清华, 陈李桥等. 随机腐蚀作用下在役钢桥疲劳抗力概率密度演化方法. 中国公路学报, 2023, 36(5): 163-174 (Jia Donglin, Zhang Qinghua, Chen Liqiao, et al. Probability density evolution method of Fatigue resistance of existing steel Bridge under stochastic corrosion effect. China Journal of Highway and Transport, 2023, 36(5): 163-174 (in Chinese) Jia Donglin, Zhang Qinghua, Chen Liqiao, et al. Probability density evolution method of Fatigue resistance of existing steel Bridge under stochastic corrosion effect. China Journal of Highway and Transport, 2023, 36(5): 163-174 (in Chinese)
[21] Xu YZ. Fatigue reliability evaluation using probability density evolution method. Probabilistic Engineering Mechanics, 2015, 42: 1-6
[22] Xiao X, Yan Y, Chen B. Stochastic dynamic analysis for vehicle-track-bridge system based on probability density evolution method. Engineering Structures 2019, 188: 745-761
[23] Peng YB, Ma YY, Huang TC, et al. Reliability-based design optimization of adaptive sliding base isolation system for improving seismic performance of structures. Reliability Engineering & System Safety, 2021, 205: 107167
[24] Song YP, Basu B, Zhang ZL, et al. Dynamic reliability analysis of a floating offshore wind turbine under wind-wave joint excitations via probability density evolution method. Renewable Energy, 2020, 168: 991-1014
[25] 蒋仲铭, 李杰. 三类随机系统广义概率密度演化方程的解析解. 力学学报, 2016, 48(2): 413-421 (Jiang Zhongming, Li Jie. Analytical solutions of generalized probability density evolution equation of three classes of stochastic systems. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 413-421 (in Chinese) Jiang Zhongming, Li Jie. Analytical solutions of generalized probability density evolution equation of three classes of stochastic systems. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 413-421 (in Chinese)
[26] Papadopoulos V, Kalogeris I. A Galerkin-based formulation of the probability density evolution method for general stochastic finite element systems. Computational Mechanics, 2016, 57: 701-716
[27] Tao WF, Li J. A difference-wavelet method for solving generalized density evolution equation in stochastic structural analysis. International Journal of Structural Stability and Dynamics, 2016, 17: 1750055
[28] Gana S. Numerical computation of spectral solutions for sturm-liouville eigenvalue problems. International Journal of Analysis and Applications, 2023, 21: 86
[29] Huang WZ, Sloan DM. The pseudospectral method for solving differential eigenvalue problems. Journal of Computational Physics, 1994, 111: 399-409
[30] Arioli M, Codenotti B, Fassino C. The Padé method for computing the matrix exponential. Linear Algebra and Its Applications, 1996, 240: 111-130
[31] Yu B, Ning CL, Li B. Hysteretic model for shear-critical reinforced concrete columns. Journal of Structural Engineering, 2016, 142: 04016056
计量
- 文章访问数: 128
- HTML全文浏览量: 59
- PDF下载量: 30
