THE CLOSED-FORM SOLUTION OF STEADY STATE RESPONSE OF HYSTERETIC SYSTEM UNDER STOCHASTIC EXCITATION
-
摘要: 滞迟系统属于一类典型的强非线性系统，滞迟力不仅取决于系统的瞬时变形，还与变形历程有关。虽然滞迟系统的随机振动问题已被广泛研究，但至今尚未得到滞迟系统随机响应概率密度函数的精确闭合解。本文运用迭代加权残值法获得了高斯白噪声激励下Bouc-Wen滞迟系统稳态响应概率密度函数的近似闭合解。首先，运用等效线性化法求出系统的稳态高斯概率密度函数；然后以此构造权函数，应用加权残值法求得了系统指数多项式形式的非高斯概率密度函数；最后引入迭代的过程，逐步优化权函数，提高计算所得结果的精度。以随机地震激励下钢纤维陶粒混凝土结构的稳态响应作为算例，其中Bouc-Wen模型的参数是基于拟静力学试验数据，并应用最小二乘法辨识获得。与Monte Carlo模拟结果相比，等效线性化法得到的结果精度较差；由加权残值法得到的结果能够表现出非线性特征，但其精度依然无法令人满意；采用迭代加权残值法得到的近似闭合解与Monte Carlo模拟的结果吻合非常好；对于较强随机激励情形，采用渐进迭代加权残值法具有较高的求解效率，所获得的理论解析解具有较高的精度。结果表明，所获得的近似闭合解不仅对于土木工程领域具有重要的实际应用价值，而且还可作为检验其他非线性系统随机响应预测方法的精度的标准。
-
关键词:
- 滞迟系统 /
- 稳态响应 /
- 闭合解 /
- 迭代加权残值法 /
- 钢纤维陶粒混凝土结构
Abstract: The hysteretic system is one of the typical strongly nonlinear systems. Hysteretic force depends not only on the instantaneous deformation but also on the past history of deformation. In the last few decades, random vibration of hysteretic system has been studied extensively, but no closed-form solution of random response of hysteretic systems is available so far. In this paper, the newly developed nonlinear random vibration scheme called iterative method of weighted residuals is explored to obtain the closed-form solution of steady-state probability density function (PDF) of the Bouc-Wen hysteretic system under Gaussian white noise excitation. First, a Gaussian PDF is obtained with equivalent linearization technique, which is used as a weighting function. Then, the method of weighted residuals is utilized to determine the non-Gaussian PDF of exponential polynomial type. Finally, an iterative procedure is introduced to improve the accuracy of the solutions obtained from the method of weighted residuals. As an illustrative example, the steadystate stochastic response of the steel fiber reinforced ceramsite concrete column under random excitation is studied, in which the hysteretic parameters associated with Bouc-Wen hysteretic model is identified from the pseudo-static test by using the method of least square. Compared to the Monte Carlo results, the accuracy of results obtained from equivalent linearization method is poor. The results obtained from weight residue method can show the nonlinearity of Bouc-Wen systems, but its accuracy is still unsatisfactory. The iterative method of weight residuals can lead to results with higher accuracy. In the case of stronger random excitation, the progressive iterative method of weighted residuals has high efficiency. The obtained solutions agree well with the Monte Carlo simulation data. The proposed closed-form solution of PDF of Bouc-Wen hysteretic system not only is significant to the civil engineering, but also can be a benchmark to examine the accuracy of solutions obtained by other methods. -
引言
滞迟现象如弹塑性[1-3]、铁电性[4]、形状记忆合金材料[5]等常出现在科学研究和工程实际的不同领域，强烈震动载荷下的结构系统通常表现出滞迟现象[6-7].学术界提出了许多模型描述这些滞迟关系，如双线性模型[8-9]，Ramberg-Osgood模型[10]，Bouc-Wen模型[6-7]，Ozdemir模型[11]，Masing模型[12]，Duhem模型[13]，Preisach模型[14]等，其中Bouc-Wen模型是较为通用的一种.
滞迟力不仅与系统当前的状态有关，而且还与系统过去的状态有关.因此，滞迟动力学系统属于一类典型的强非线性系统.近年来，滞迟系统的随机振动问题已被广泛研究，但至今尚未获得精确的闭合解，因此，学术界提出了各种近似解法.如，Iwan[15]采用基于克雷洛夫-包哥留波夫假设的等效线性化法研究了双线性滞迟系统的随机响应. Roberts[16-17]分别采用标准随机平均法及能量包线随机平均法研究了双线性滞迟系统的随机响应. Bouc[7]釆用FPK方程法研究了滞迟系统的随机响应. Wen[7, 18]采用线性化方法研究了滞迟系统的随机响应. Zhu等[19]采用能量包线随机平均法研究分布弹塑性元件为滞迟恢复力模型的滞迟系统的随机响应. Lin等[20]采用能量包线随机平均法研究了Bouc-Wen滞迟系统的随机响应. Ying等[21]采用能量包线随机平均法研究了Duhem滞迟系统的随机响应.对具有非局部记忆特性的Preisach滞迟系统，Mayergoyz和Korman[22]研究了随机输入下Preisach系统的平均输出，Ni等[23]基于方差和非局部记忆滞迟本构模型的切换概率分析，近似得到了Preisach滞迟系统的非线性随机动力响应的二阶统计量. Spanos等[24]和Wang等[25]分别采用标准随机平均和能量包线随机平均法研究了Preisach滞迟系统的随机响应.最近，Jin等[26]运用随机平均法和伽辽金法获得了随机地震激励下滞迟系统的近似瞬态响应.
然而，上述方法都存在一些不足，如随机平均法仅限于能量耗散很小且弱激励情形；等效线性化法能得到较准确的均方速度和均方位移，但是这种方法只局限于高斯统计的情形，在参激情况下常被认为是不充足或不合适的；非高斯闭合法在尾部通常会获得负值概率，特别是参激情形.因此，还需进一步开展滞迟非线性系统的随机响应预测研究.
Er [27-29]提出了一种指数多项式闭合法求解稳态FPK方程，但由于采用了高斯概率密度函数构造权函数，因此指数多项式闭合解法的应用范围受到了很大限制. Di Paola和Sofi [30]改进了指数多项式闭合法，提出采用一种简单有效的迭代过程以提高近似解的精确性，但该方法没有从根本上解决指数闭合解的局限性问题.最近，文献[31]提出了一种求解FPK方程稳态解的新方法------迭代加权残值法，该方法的核心思想是逐步优化权函数，目前已被成功应用于求解多种复杂单自由度强非线性系统[31-32].
本文将迭代加权残值法进一步应用于Bouc-Wen滞迟系统，构造其稳态概率密度函数.首先，将稳态FPK方程的解设为指数多项式，然后借助加权残值法确定假设解中的待定系数；引入迭代方法，逐步优化权函数，提高加权残值法的精度，得到系统稳态响应的近似闭合解.作为应用，本文获得了高斯白噪声激励下钢纤维陶粒混凝土框架结构的近似稳态响应概率密度函数闭合解.研究表明，所获得的闭合解与Monte Carlo模拟结果吻合得较好.
1. 迭代加权残值法
迭代加权残值法主要由两个步骤组成.首先，应用加权残值法获得系统的近似稳态概率密度响应；然后，引入迭代过程，逐步优化权函数，提高加权残值法的精度，最终得到具有较高精度的近似闭合稳态解.
1.1 加权残值法
考虑外激高斯白噪声激励下的Bouc-Wen系统，其运动方程如下
$$\left. \begin{align} &\frac{\text{d}X}{\text{d}t}=Y \\ &\frac{\text{d}Y}{\text{d}t}=-2\xi Y-\alpha X-(1-\alpha )Z+W(t) \\ &\frac{\text{d}Z}{\text{d}t}=\lambda Y-\beta {{\left| Z \right|}^{n-1}}Z\left| Y \right|-\gamma {{\left| Z \right|}^{n}}Y \\ \end{align} \right\}$$ (1) 式中，X, Y与Z分别表示系统位移、速度与滞迟力；ξ为黏滞阻尼率；α∈(0; 1) 是屈服前后的刚度比；λ，β，γ与n为滞回环参数，β与γ控制滞回环的形状，λ控制滞回环的幅值，n控制滞回环曲线的光滑性；W(t)是强度为2Dδ(t)的高斯白噪声.
与系统(1) 相应的稳态FPK方程为
$$\begin{align} &\frac{-\partial \left( py \right)}{\partial x}+\frac{\partial \left\{ \left[ 2\xi y+\alpha x+\left( 1-\alpha \right)z \right]p \right\}}{\partial y}- \\ &\qquad \frac{-\partial \left[ \left( \lambda y-\beta z{{\left| z \right|}^{n-1}}\left| y \right|-\gamma y{{\left| z \right|}^{n}} \right)p \right]}{\partial z}+D\frac{{{\partial }^{2}}p}{\partial {{y}^{2}}}=0 \\ \end{align}$$ (2) 式中，p = p(x, y, z).目前，式(2) 尚未获得精确解析解.
构造如下形式的近似解
$$\bar{p}(x,y,z)={{C}_{0}}\exp [\varphi (x,y,z,{{c}_{ijk}})]$$ (3) 其中，C0为归一化常数，φ(x, y, z, cijk)为关于状态变量的n阶多项式，可表示为
$$\varphi \left( x,y,z,{{c}_{ijk}} \right)=\sum{{{c}_{ijk}}{{x}^{i}}{{y}^{j}}{{z}^{k}}}$$ (4) 其中，cijk为待求系数. p(x, y, z)的存在条件为
$$\left. \begin{align} &\varphi (x,y,z,{{c}_{ijk}})＜0,\ \ \left| x \right|+\left| y \right|+\left| z \right|\to \infty \\ &\varphi (x,y,z,{{c}_{ijk}})\to {{R}^{\mu }}, \\ &\ mu＞-1,\ \ \left| x \right|+\left| y \right|+\left| z \right|\to 0 \\ \end{align} \right\}$$ (5) 其中，x = R sin θ cos φ, y = R sin θ sin φ, z = R cos θ.将方程(3) 代入式(2) 中，得残差
$$\begin{align} &r\left( x,y,z,{{c}_{ijk}} \right)=\frac{-\partial \varphi }{\partial x}+\left[ 2\xi y+\alpha x+\left( 1-\alpha \right)z \right]\frac{\partial \varphi }{\partial y}- \\ &\qquad \left( \lambda y-\beta z{{\left| z \right|}^{n-1}}\left| y \right|-\gamma y{{\left| z \right|}^{n}} \right)\frac{\partial \varphi }{\partial z}+D\frac{{{\partial }^{2}}\varphi }{\partial {{y}^{2}}}+D{{\left( \frac{\partial \varphi }{\partial y} \right)}^{2}}+ \\ &\qquad 2\xi +\beta \left| y \right|{{\left| z \right|}^{n-1}}+n\gamma y{{\left| z \right|}^{n-1}}\text{sgn}(z) \\ \end{align}$$ (6) 由于p(x, y, z)只是p(x, y, z)的近似值.因此，残差r(x, y, z, cijk)通常不为零.根据加权残值法，引入一组权函数Mijk(x, y, z)，使其与残差r(x, y, z, cijk)的乘积在域内的积分为零，即
$$\underset{-\infty }{\overset{+\infty }{\mathop{\int }}}\,\underset{-\infty }{\overset{+\infty }{\mathop{\int }}}\,\underset{-\infty }{\overset{+\infty }{\mathop{\int }}}\,{{M}_{ijk}}(x,y,z)r\left( x,y,z,{{c}_{ijk}} \right)\text{d}x\text{d}y\text{d}z=0$$ (7) $${{M}_{ijk}}\left( x,y,z \right)={{p}_{m}}\left( x,y,z \right){{x}^{i}}{{y}^{j}}{{z}^{k}}$$ (8) 其中, pm(x, y, z)是由等效线性化法或随机平均法等方法获得的系统概率密度函数.值得注意的是，在数值计算时方程(7) 中的积分域通常为有限积分域，但通常难以确定.为此，运用Monte Carlo模拟法粗略估计积分域Ωs.方程(7) 可近似表示为
$$\iiint\limits_{{{\Omega }_{\rm{s}}}}{{{p}_{m}}{{x}^{i}}{{y}^{j}}{{z}^{k}}r\left( x,y,z,{{c}_{ijk}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}\rm{=0}}$$ (9) 其中，i = 0, 1, 2, …；j = 0, 1, 2, …, k = 0, 1, 2, …，即
$$\left. \begin{align} & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{i}}{{y}^{j}}{{z}^{k}}r\left( x,y,z,{{c}_{ijk}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{1}}{{y}^{0}}{{z}^{0}}r\left( x,y,z,{{c}_{100}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{0}}{{y}^{1}}{{z}^{0}}r\left( x,y,z,{{c}_{010}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{0}}{{y}^{0}}{{z}^{1}}r\left( x,y,z,{{c}_{001}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{1}}{{y}^{1}}{{z}^{0}}r\left( x,y,z,{{c}_{110}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \iint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{1}}{{y}^{0}}{{z}^{1}}r\left( x,y,z,{{c}_{101}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ & \quad \quad \quad \quad \quad \quad \quad \cdots \\ & \iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{p}_{m}}{{x}^{i}}{{y}^{j}}{{z}^{k}}r\left( x,y,z,{{c}_{ijk}} \right)\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}=0} \\ \end{align} \right\}$$ (10) 其中，0＜i + j + k≤l, l为解的阶数.将式(10) 进行数值积分，即可得到一组关于cijk的二次非线性代数方程组.采用牛顿法求解该方程组，可获得系统的近似稳态响应概率密度.
1.2 迭代过程
滞迟系统属于强非线性系统，仅使用一次加权残值法可能得不到具有足够精度的解.因此，采用文献[30-31]中的迭代技术来提高解的精度.令p(k)为迭代k次后得到的近似稳态概率密度函数，取代式(9) 中的pm，再应用加权残值法求出下一个近似稳态概率密度函数p(k+1).重复此步骤直到满足以下收敛条件
$$\begin{align} & \varepsilon =\sqrt{\iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{({{{\bar{p}}}^{\left( k \right)}}-{{{\bar{p}}}^{(k+1)}})}^{2}}\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}}}\approx \\ & \qquad \sqrt{\sum\limits_{i=0}^{{{N}_{1}}}{\sum\limits_{j=0}^{{{N}_{2}}}{\sum\limits_{k=0}^{{{N}_{3}}}{{{({{{\bar{p}}}^{(k)}}-{{{\bar{p}}}^{(k+1)}})}^{2}}\Delta x\Delta y\Delta z}}}}\le {{\varepsilon }_{0}} \\ \end{align}$$ (11) 其中，ε0为预设误差，N1, N2与N3为常数，分别表示状态空间离散的数目，∆x, ∆y和∆z表示离散积分步长.
为了估计解p(k)的精度∆pk，引入以下标准
$$\begin{align} & \Delta {{p}_{k}}=\sqrt{\iiint\limits_{{{\mathit{\Omega }}_{\rm{s}}}}{{{({{{\bar{p}}}^{\left( k \right)}}-{{p}_{R}})}^{2}}\rm{d}\mathit{x}\rm{d}\mathit{y}\rm{d}\mathit{z}}}\approx \\ & \qquad \sqrt{\sum\limits_{i=0}^{{{N}_{1}}}{\sum\limits_{j=0}^{{{N}_{2}}}{\sum\limits_{k=0}^{{{N}_{3}}}{{{({{{\bar{p}}}^{(k)}}-{{p}_{R}})}^{2}}\Delta x\Delta y\Delta z}}}} \\ \end{align}$$ (12) 其中，pR(x, y, z)表示Monte Carlo模拟得到的稳态概率密度函数.
需要指出的是, 当pm与真实解相差较大时，上述迭代过程可能不会收敛，特别是当系统为强非线性的情况.文献[31]提出了一种渐近迭代的方案，即先应用加权残值法获得系统弱非线性情形时或弱阻尼弱激励情形时的近似稳态概率密度函数，然后以此构造权函数，再在非线性参数空间渐近迭代.例如，针对滞迟系统(1)，为了获得激励强度D =0.2情形时的近似稳态响应概率密度函数, 首先采用迭代加权残值法获得弱随机激励D =0.1情形时的近似稳态概率函数pm1；然后将pm1代入式(9)，以此构造新的权函数，再次利用迭代加权残值法，求得随机激励D =0.15情形时的近似稳态概率函数pm2；最后将pm2代入式(9)，构造新的权函数，利用迭代加权残值法求解D =0.2时的概率密度函数.这种渐近迭代的方法可有效避免选取的初始值pm与真实解相差较大时导致的不收敛情况，提高了求解的效率.
2. 算例
考虑如图 1所示的水平随机地震作用下单自由度钢纤维陶粒混凝土框架滞迟系统，其平衡条件为
$$m{{{\dot{x}}}_{t}}+c\dot{x}+g(x,\dot{x})=0$$ (13) ![]()
图 1 地震地面运动下单自由度钢钎维陶粒混凝土框架结构Figure 1. Steel fiber reinforced ceramsite concrete (SFRCC) single degree of freedom frame subjected to earthquake ground motion其中，x(t)是集中质量的相对位移；xt(t)是集中质量的总位移；$g(x,\dot{x})$是恢复力.
由图 1中可得
$${{x}_{t}}=x+{{x}_{g}}$$ (14) $${{{\ddot{x}}}_{t}}=\ddot{x}+{{{\ddot{x}}}_{g}}$$ (15) 式中，xg和$\ddot {x}_g $分别表示地面的位移和加速度, 将方程(15) 代入方程(13) 可得
$$m\ddot{x}+m{{{\ddot{x}}}_{g}}+c\dot{x}+g(x,\dot{x})=0$$ (16) 或
$$m\ddot{x}+c\dot{x}+g(x,\dot{x})=F(t)$$ (17) 其中，F(t)为图 2所示导致水平地面加速度的等效载荷.该框架结构被简化为如图 2所示的等效系统.
将方程(17) 无量纲化，并把等效载荷理想化为独立的高斯白噪声.该系统的动力学方程化为
$$\ddot{x}+2\xi \dot{x}+\alpha x+(1-\alpha )z={{W}_{1}}(t)$$ (18) 其中，x为无量纲化的位移, ξ 为黏滞阻尼率, α是屈服前后的刚度比, k1为常数, W1(t)是强度为2D1的高斯白噪声, z为恢复力中依赖于时间的滞迟力，由以下公式表示
$$\dot{z}=\lambda \dot{x}-\beta z\left| {\dot{x}} \right|{{\left| z \right|}^{n-1}}-\gamma \dot{x}{{\left| z \right|}^{n}}$$ (19) 本文基于钢纤维陶粒混凝土柱受水平低周反复加载试验所得的数据，经最小二乘法辨识获得滞迟系统的参数: n =1, λ= 1.103 7, β=0.202和γ=0.314 7.理论曲线和试验曲线对比如图 3所示，二者吻合较好.系统的其他参数为: ξ=0.025, α=0.2, D=0.1, ε0 =0.001.
![]()
图 3 钢纤维陶粒混凝土柱完全塑性阶段试验曲线和理论曲线(实线表示理论结果，虚线表示试验结果)Figure 3. Theoretical curve and practical curve of SFRCC column in fully plastic stage (the solid line denotes the theoretical curve and the dotted line denotes the practical curve)最终可得如下形式的运动方程
$$\left. \begin{align} &\frac{\text{d}x}{\text{d}t}=y \\ &\frac{\text{d}y}{\text{d}t}=-0.05y-0.2x-0.8z+W(t) \\ &\frac{\text{d}z}{\text{d}t}=1.1037\dot{x}-0.202\left| {\dot{x}} \right|z-0.3147\dot{x}\left| z \right| \\ \end{align} \right\}$$ (20) 本文首先运用等效线性化法获得了系统的近似稳态响应概率密度函数，表达式如下
$$\begin{align} &{{p}_{0}}(x,y,z)=0.403380833\exp (2.559927296yz+ \\ &\qquad 13.60167586xz-2xy-1.026369574{{y}^{2}}- \\ &\qquad 6.00480473{{x}^{2}}-8.704825325{{z}^{2}}) \\ \end{align}$$ (21) 然后以式(21) 构造权函数，运用加权残值法获得了新的近似稳态响应概率密度函数，即
$$\begin{align} &{{p}_{1}}(x,y,z)=0.0880074433\exp (-0.1191584596{{x}^{2}}+ \\ &\qquad 0.3875850584xy-0.6653184665xz- \\ &\qquad 0.126015634{{y}^{2}}-0.804881587yz+ \\ &\qquad 1.597915482{{z}^{2}}-0.1550420718{{y}^{4}}+ \\ &\qquad 0.5391715902{{y}^{3}}z-1.288822698{{y}^{2}}{{z}^{2}}+ \\ &\qquad 2.259549091y{{z}^{3}}-3.723961542{{z}^{4}}- \\ &\qquad 0.8416906504{{x}^{3}}y-0.4595036571{{x}^{2}}{{y}^{2}}- \\ &\qquad 0.287403815x{{y}^{3}}-0.8905259629{{x}^{4}}+ \\ &\qquad 4.535946201{{x}^{3}}z-8.929948156{{x}^{2}}{{z}^{2}}+ \\ &\qquad 8.374445343x{{z}^{3}}+3.217731206{{x}^{2}}yz+ \\ &\qquad 1.285962464x{{y}^{2}}z-4.307610477xy{{z}^{2}}) \\ \end{align}$$ (22) 根据误差公式(11)，式(21) 与式(22) 的误差分别为0.163和0.249.显然，式(22) 的精度低于式(21).因此，为了进一步提高精度，现以解(22) 构造权函数，经过1次迭代，获得了误差为0.035的解
$$\begin{align} &{{p}_{4}}(x,y,z)=0.1996305707\exp (-3.377531344{{x}^{2}}- \\ &\qquad 0.5008616187xy+5.763301455xz- \\ &\qquad 0.3803141420{{y}^{2}}+0.03928836384yz- \\ &\qquad 1.703843746{{z}^{2}}-0.1739958276{{y}^{4}}+ \\ &\qquad 0.6570314190{{y}^{3}}z-1.070008642{{y}^{2}}{{z}^{2}}+ \\ &\qquad 1.084530978y{{z}^{3}}-2.724243965{{z}^{4}}- \\ &\qquad 0.2460616857{{x}^{3}}y-0.04611532747{{x}^{2}}{{y}^{2}}- \\ &\qquad 0.2945573801x{{y}^{3}}-0.3305917948{{x}^{4}}+ \\ &\qquad 0.5676808287x{{y}^{2}}z-1.088846865xy{{z}^{2}}) \\ \end{align}$$ (23) 图 4给出了关于p(x, z)的稳态边缘概率密度函数.由图 4和图 5可知，等效线性化的结果与Monte Carlo模拟结果相差甚远；加权残值法的计算结果精度尚未令人满意，但已经表现出系统的一些非线性特征；迭代加权残值法得到结果与Monte Carlo模拟的结果吻合较好.
![]()
图 4 D =0.1情形时关于p(x, z)的稳态边缘概率密度函数Figure 4. The steady-state marginal probability density function p(x, z) in case of D =0.1![]()
图 5 D =0.1情形时的稳态边缘概率密度函数p1(x)和p2(z) ($\circ, *$表示Monte Carlo模拟数据)Figure 5. The steady-state marginal probability density function p1(x) and p2(z) in case of D =0.1 ($\circ, *$ represent the Monte Carlo simulation data)另外，以式(23) 构造权函数，采用渐进迭代加权残值法，经过3次迭代获得了较强随机激励D= 0.2情形时系统的近似稳态概率密度函数闭合解
$$\begin{align} &{{p}_{3}}(x,y,z)=0.0634854305\exp (-0.715892171{{x}^{2}}- \\ &\qquad 0.0367187464xy+0.9442267208xz- \\ &\qquad 0.1393982600{{y}^{2}}-0.2161024989yz+ \\ &\qquad 0.2541890718{{z}^{2}}-0.1374801508{{y}^{4}}+ \\ &\qquad 0.1909194619{{y}^{3}}z-0.3314499778{{y}^{2}}{{z}^{2}}+ \\ &\qquad 0.5168474132y{{z}^{3}}-0.8904109436{{z}^{4}}- \\ &\qquad 0.1193061758{{x}^{3}}y-0.0310110389{{x}^{2}}{{y}^{2}}- \\ &\qquad 0.0860371179x{{y}^{3}}-0.1041557739{{x}^{4}}+ \\ &\qquad 0.514458620{{x}^{3}}z-0.9276493809{{x}^{2}}{{z}^{2}}+ \\ &\qquad 1.043816700x{{z}^{3}}+0.4036708712{{x}^{2}}yz+ \\ &\qquad 0.2010678510x{{y}^{2}}z-0.6387509013xy{{z}^{2}}) \\ \end{align}$$ (24) 式(24) 的精度为0.040. 图 6给出D=0.2情形时的稳态边缘概率密度函数. 图 6(a)与6(b)表示关于p(x, z)的稳态边缘概率密度函数. 图 6(c)与图 6 (d)分别表示关于p1(x)及p2(z)的稳态边缘概率密度及均方差.符号$(\circ, *)$表示样本数为40 000的Monte Carlo模拟结果.由图 6可知，理论解析解与模拟结果在较强随机激励情形时仍吻合得较好.
![]()
图 6 D =0.2情形时稳态边缘概率密度函数($\circ, *$表示Monte Carlo模拟结果)Figure 6. The steady-state marginal probability density function in case of D =0.2 ($\circ, *$ represent results from Monte Carlo simulation)3. 结论
本文首次获得Bouc-Wen滞迟系统的稳态响应概率密度函数的近似闭合解.首先，利用等效线性化法获得的高斯概率密度函数构造权函数，然后应用加权残值法获得系统指数多项式形式的非高斯概率密度函数，最后引入迭代的办法提高了加权残值法计算所得结果的精度.对较强随机激励情形，以弱随机激励情形时的稳态响应概率密度函数构造权函，再在随机激励的参数空间内渐近迭代.作为算例，本文研究了随机激励下钢纤维陶粒混凝土结构的稳态响应，其中Bouc-Wen系统的参数是基于拟静力学试验由最小二乘法辨识获得.对比Monte Carlo模拟结果，等效线性化的结果精度较差，加权残值法的结果虽表现出了系统的非线性特征但精度尚未令人满意，迭代加权残值法得到结果与Monte Carlo模拟的结果吻合得较好.对于较强随机激励情形，基于渐近迭代加权残值法具有较高的效率，且所得结果也具有较好的精度.本文所获得的闭合解，不仅对于土木工程领域具有重要的实际应用价值，还可作为一个标准，用来检验其他非线性系统随机响应预测方法的精度.
-
![]()
图 1 地震地面运动下单自由度钢钎维陶粒混凝土框架结构
Figure 1. Steel fiber reinforced ceramsite concrete (SFRCC) single degree of freedom frame subjected to earthquake ground motion
![]()
图 3 钢纤维陶粒混凝土柱完全塑性阶段试验曲线和理论曲线(实线表示理论结果，虚线表示试验结果)
Figure 3. Theoretical curve and practical curve of SFRCC column in fully plastic stage (the solid line denotes the theoretical curve and the dotted line denotes the practical curve)
![]()
图 4 D =0.1情形时关于p(x, z)的稳态边缘概率密度函数
Figure 4. The steady-state marginal probability density function p(x, z) in case of D =0.1
![]()
图 5 D =0.1情形时的稳态边缘概率密度函数p1(x)和p2(z) ($\circ, *$表示Monte Carlo模拟数据)
Figure 5. The steady-state marginal probability density function p1(x) and p2(z) in case of D =0.1 ($\circ, *$ represent the Monte Carlo simulation data)
-
[1] 万征, 姚仰平, 孟达.复杂加载下混凝土的弹塑性本构模型.力学学报, 2016, 48(5):1159-1171 http://lxxb.cstam.org.cn/CN/abstract/abstract146013.shtml Wan Zheng, Yao Yangping, Meng Da. An elastoplastic constitutive model of concrete of complicated load. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5):1159-1171 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract146013.shtml
[2] 李龙彪.纤维增强陶瓷基复合材料疲劳迟滞回线模型研究.力学学报, 2014, 46(5):710-729 doi: 10.11776/cjam.31.05.A022 Li Longbiao. Investigation on fatigue hysteresis loops models of fibre-reinforced ceramic-matrix composites. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(5):710-729 (in Chinese) doi: 10.11776/cjam.31.05.A022
[3] 郭洪宝, 贾普荣, 王波等.基于迟滞行为的2D-SiC/SiC复合材料组份力学性能分析.力学学报, 2015, 47(2):260-269 doi: 10.6052/0459-1879-14-247 Guo Hongbao, Jia Purong, Wang Bo et al. Study on constituent properties of a 2D-SiC/SiC composite by hysteresis measurements. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(2):260-269 (in Chinese) doi: 10.6052/0459-1879-14-247
[4] 万强, 陈常青, 沈亚鹏.铁电陶瓷PZT53复杂力电耦合行为的实验研究.力学学报, 2005, 37(4):413-420 http://lxxb.cstam.org.cn/CN/abstract/abstract141342.shtml Wan Qiang, Chen Changqing, Shen Yapeng. An experimental investigation into the complex electromechanical behavior of PZT53. Chinese Journal of Theoretical and Applied Mechanics, 2005, 37(4):413-420 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract141342.shtml
[5] Mayergoyz I. Mathematical models of hysteresis. IEEE Transactions on Magnetic, 1986, 22(5):603-608 doi: 10.1109/TMAG.1986.1064347
[6] Wen YK. Methods of random vibrationfor inelastic structures. Applied Mechanics Review, 1989, 42(2):39-52 doi: 10.1115/1.3152420
[7] Bouc R. Forced vibration of mechanical systems with hysteresis//Proceedings of the Fourth Conference on Non-Linear Oscillation, Prague, Czechoslovakia, 1967
[8] Caughey T. Random excitation of a system with bilinear hysteresis. ASME Journal of Applied Mechanics, 1960, 27(4):649-652 doi: 10.1115/1.3644077
[9] Suzuki Y, Minai R. Application of stochastic differential equations to seismic reliability analysis of hysteretic structures. Probabilistic Engineering Mechanics, 1988, 3(1):43-52 doi: 10.1016/0266-8920(88)90007-0
[10] Jennings PC. Periodic response of a general yielding structure. Journal of Engineering Mechanics-ASCE Division, 1964, 90(2):131-166
[11] Bhatti M, Pister K. A dual criteria approach for optimal-design for earthquake-resistant structural systems. Earthquake Engineering Structure Design, 1981, 9(6):557-572 https://www.researchgate.net/publication/290799999_DUAL_CRITERIA_APPROACH_FOR_OPTIMAL_DESIGN_OF_EARTHQUAKE-RESISTANT_STRUCTURAL_SYSTEMS
[12] Beck JL, Jayakumar P. Class of Masing models for plastic hysteresis in structures//Proceedings of 14th ASCE Structures Congress, 1996:1083-1090
[13] Visintin A. Differential Models of Hysteresis Ⅲ. Springer Science & Business Media, 2013
[14] Mayergoyz I. Mathematical models of hysteresis. IEEE Transactions on Magnetic, 1986, 22(5):603-608 doi: 10.1109/TMAG.1986.1064347
[15] Iwan WD. A distributed-element model for hysteresis and its steadystate dynamic response. ASME Journal of Applied Mechanics, 1966, 33:893 doi: 10.1115/1.3625199
[16] Roberts JB. The response of an oscillator with bilinear hysteresis to stationary random excitation. ASME Journal of Applied Mechanics, 1978, 45(4):923-928 doi: 10.1115/1.3424442
[17] Roberts JB. Application of averaging methods to randomly excited hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems, Berlin, Springer, 1988
[18] Wen YK. Equivalent linearization for hysteretic systems under random excitation. ASME Journal of Applied Mechanics, 1980, 47(1):150-154 doi: 10.1115/1.3153594
[19] Zhu WQ, Lei Y. Stochastic averaging of energy envelope of bilinear hysteretic systems//Nonlinear Stochastic Dynamic Engineering Systems, Berlin:Springer, 1988
[20] Lin YK, Cai GQ. Probabilistic Structural Dynamics:Advanced Theory and Applications. New York:McGraw-Hill, 1995
[21] ZG Ying, WQ Zhu, Ni YQ, et al. Stochastic averaging of Duhem hysteretic systems. Journal of Sound and Vibration, 2002, 254(1):91-104 doi: 10.1006/jsvi.2002.4086
[22] Mayergoyz ID, Korman CE. The Preisach model with stochastic input as a model for aftereffect. Journal of Applied Physics, 1994, 75(10):5478-5480 doi: 10.1063/1.355712
[23] Ni YQ, Ying ZG, Ko JM. Random response analysis of Preisach hysteretic systems with symmetric weight distribution. ASME Journal of Applied Mechanics, 2002, 69(2):171-178 doi: 10.1115/1.1428333
[24] Spanos PD, Cacciola P, Muscolino G. Stochastic averaging of Preisach hysteretic systems. ASCE Journal of Engineering Mechanics, 2004, 130(11):1257-1267 doi: 10.1061/(ASCE)0733-9399(2004)130:11(1257)
[25] Wang Y, Ying ZQ, Zhu WQ. Stochastic averaging of energy envelope of Preisach hysteretic systems. Journal of Sound and Vibration, 2009, 321(3-5):976-993 doi: 10.1016/j.jsv.2008.10.021
[26] Jin XL, Wang Y, Huang ZL. Analysis and control for transient responses of seismic-excited hysteretic structures. Soil Dynamics and Earthquake Engineering, 2015, 73(6):58-65 http://www.sciencedirect.com/science/article/pii/S0267726115000524
[27] Er GK. Exponential closure method for some randomly excited non-linear systems. International Journal of Nonlinear Mechanics, 2000, 35(1):69-78 doi: 10.1016/S0020-7462(98)00088-2
[28] Er GK. A consistent method for the solution to reduced FPK equation in statistical mechanics. Physica A, 1999, 262(1):118-128 http://www.sciencedirect.com/science/article/pii/S0378437198003628
[29] Er GK. An improved closure method for analysis of nonlinear stochastic systems. Nonlinear Dynamics, 1998, 17(3):285-297 doi: 10.1023/A:1008346204836
[30] Di Paola M, Sofi A. Approximate solution of the Fokker-PlanckKolmogorov equation. Probabilistic Engineering Mechanics, 2002, 17(4):369-384 doi: 10.1016/S0266-8920(02)00034-6
[31] Chen LC, Sun JQ. The closed-form solution of the reduced FokkerPlanck-Kolmogorov equation for nonlinear systems. Communications in Nonlinear Science & Numerical Simulation, 2016, 41(12):1-10 http://www.sciencedirect.com/science/article/pii/S1007570416300879
[32] Chen LC, Sun JQ, The closed-form steady-state probability density function of van der pol oscillator under random excitations. Journal of Applied Nonlinear Dynamics, 2016, 5(4):495-502 doi: 10.5890/JAND.2016.12.009
-
期刊类型引用(9)
1. Lincong CHEN，Zi YUAN，Jiamin QIAN，J.Q.SUN. Random vibration of hysteretic systems under Poisson white noise excitations. Applied Mathematics and Mechanics(English Edition). 2023(02): 207-220 . 
必应学术
2. 周碧柳，靳艳飞. 高斯色噪声和谐波激励共同作用下耦合SD振子的混沌研究. 力学学报. 2022(07): 2030-2040 . 本站查看
3. 徐严钢，朱海涛，刘治国. 三维路径积分法的积分区间以及子区间数目选取方法研究. 应用力学学报. 2021(04): 1358-1365 . 百度学术
4. 郭秀秀，李长宇，史庆轩. 基于改进Bouc-Wen模型的非线性结构非平稳随机地震响应分析. 振动与冲击. 2020(18): 248-254+268 . 百度学术
5. 张婉洁，牛江川，申永军，杨绍普，王丽. 一类阻尼控制半主动隔振系统的解析研究. 力学学报. 2020(06): 1743-1754 . 本站查看
6. 芮珍梅，陈建兵. 加性非平稳激励下结构速度响应的FPK方程降维. 力学学报. 2019(03): 922-931 . 本站查看
7. 石晟，杜东升，王曙光，李威威. 概率密度演化方程TVD格式的自适应时间步长技术及其初值条件改进. 力学学报. 2019(04): 1223-1234 . 本站查看
8. 邢子康，申永军，李向红. 接地式三要素型动力吸振器性能分析. 力学学报. 2019(05): 1466-1475 . 本站查看
9. 邹广平，张冰，唱忠良，刘松. 弹簧-金属丝网橡胶组合减振器迟滞力学模型及实验研究. 力学学报. 2018(05): 1125-1134 . 本站查看
其他类型引用(5)
下载:
计量
- 文章访问数: 1371
- HTML全文浏览量: 232
- PDF下载量: 840
- 被引次数: 14
