留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

实时子结构试验中显式算法对比分析

唐玉 覃晖

唐玉, 覃晖. 实时子结构试验中显式算法对比分析[J]. 工程力学, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
引用本文: 唐玉, 覃晖. 实时子结构试验中显式算法对比分析[J]. 工程力学, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
Yu TANG, Hui QIN. COMPARISONS OF MODEL-BASED EXPLICIT INTEGRATION ALGORITHMS IN REAL-TIME SUBSTRUCTURE TESTING[J]. Engineering Mechanics, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
Citation: Yu TANG, Hui QIN. COMPARISONS OF MODEL-BASED EXPLICIT INTEGRATION ALGORITHMS IN REAL-TIME SUBSTRUCTURE TESTING[J]. Engineering Mechanics, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001

实时子结构试验中显式算法对比分析

doi: 10.6052/j.issn.1000-4750.2019.04.S001
基金项目: 国家自然科学基金项目(41904095,51979027);中央高校基本科研业务费项目(DUT19JC23,DUT19RC(4)024)
详细信息
    作者简介:

    唐玉(1985−),女,辽宁人,讲师,博士,主要从事结构动力分析及子结构试验技术研究(E-mail: ytang@dlut.edu.cn)

    通讯作者: 覃晖(1985−),男,陕西人,讲师,博士,主要从事隧道工程动力分析及试验技术研究(E-mail: hqin@dlut.edu.cn)
  • 中图分类号: TU317

COMPARISONS OF MODEL-BASED EXPLICIT INTEGRATION ALGORITHMS IN REAL-TIME SUBSTRUCTURE TESTING

图(7)
计量
  • 文章访问数:  109
  • HTML全文浏览量:  35
  • PDF下载量:  71
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-04-30
  • 修回日期:  2019-12-31
  • 网络出版日期:  2020-06-01
  • 刊出日期:  2020-06-01

实时子结构试验中显式算法对比分析

doi: 10.6052/j.issn.1000-4750.2019.04.S001
    基金项目:  国家自然科学基金项目(41904095,51979027);中央高校基本科研业务费项目(DUT19JC23,DUT19RC(4)024)
    作者简介:

    唐玉(1985−),女,辽宁人,讲师,博士,主要从事结构动力分析及子结构试验技术研究(E-mail: ytang@dlut.edu.cn)

    通讯作者: 覃晖(1985−),男,陕西人,讲师,博士,主要从事隧道工程动力分析及试验技术研究(E-mail: hqin@dlut.edu.cn)
  • 中图分类号: TU317

摘要: 实时子结构试验技术中的一个关键问题是求解数值子结构的动力响应,而这一过程可选取合适的数值积分算法来实现。针对目前已建立的三种基于模型的显式积分算法(Chang法、CR法和实时子结构RST法),对比分析了各算法在线性系统和非线性系统中的数值特性。结果表明:三种算法在线性系统和具有刚度软化特性的非线性系统中均是无条件稳定的,在具有刚度硬化特性的非线性系统中变为有条件稳定。当结构阻尼比为零时,三种算法均无数值阻尼,且周期延长率结果完全一致,并随Ω的增加而增大;当结构阻尼比不为零时,三种算法均存在数值阻尼,且CR法的数值阻尼绝对值最小,而RST法的周期延长率最小。两个算例表明,RST法和Chang法的精度要优于CR法,但因Chang法是半显式的,因此RST法更适于实时子结构试验中数值子结构的仿真计算。

English Abstract

唐玉, 覃晖. 实时子结构试验中显式算法对比分析[J]. 工程力学, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
引用本文: 唐玉, 覃晖. 实时子结构试验中显式算法对比分析[J]. 工程力学, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
Yu TANG, Hui QIN. COMPARISONS OF MODEL-BASED EXPLICIT INTEGRATION ALGORITHMS IN REAL-TIME SUBSTRUCTURE TESTING[J]. Engineering Mechanics, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
Citation: Yu TANG, Hui QIN. COMPARISONS OF MODEL-BASED EXPLICIT INTEGRATION ALGORITHMS IN REAL-TIME SUBSTRUCTURE TESTING[J]. Engineering Mechanics, 2020, 37(S): 1-5, 12. doi: 10.6052/j.issn.1000-4750.2019.04.S001
  • 实时子结构试验作为近年来新兴的试验技术得到了广泛关注,该技术的核心思想是基于结构动力分析中的子结构概念,将整体结构划分为若干个子结构,其中难以模拟或者具有强非线性的构件作为试验子结构进行物理模型试验,而其余部分作为数值子结构进行数值仿真计算,通过子结构之间相互作用传递,获得整体结构的动力反应[1-2]。在求解数值子结构的运动方程时,通常需要选定合适的数值积分算法以保证实时子结构试验顺利开展,因此研究数值积分算法的精度和稳定性就显得尤为重要。

    实时子结构试验的数值积分算法有显式和隐式之分,当第i+1步的位移和速度均可由第i步及其以前的解表达时称之为显式算法,反之为隐式算法[3]。显式算法的最大缺点在于条件稳定性,特别是多自由度体系子结构试验的条件稳定性对时间步长要求十分严格;隐式算法的主要优点是无条件稳定和能量耗散性能高,但缺点是需要迭代[4]。近年来,不少学者都致力于开展无条件稳定显式积分算法的研究,提出了几种基于模型的积分算法,即算法的积分参数与结构模型有关。2002年,台湾学者Chang[5-6]最先提出了基于模型的无条件稳定半显式积分算法Chang法,其位移表达式为显式,但速度表达式因含有下一步的加速度而为隐式;2008年,Chen和Ricles[7]利用离散传递函数的概念构造了一种新算法−CR法,实现了位移和速度的显式表达;随后,部分学者[8-9]保留CR法的位移和速度表达式,在积分参数的表达式中引入控制参数,扩大了CR法的应用范围;Kolay和Ricles[10]则利用离散控制理论提出了具有可控数值阻尼的无条件稳定显式KR-α法,该算法可以为高阶模态提供数值阻尼,而忽略低阶模态的数值阻尼,但计算过程十分复杂;Tang和Lou[11-12]在CR法的基础上进一步改进,提出了一种实时子结构试验RST法,并验证了该算法在线性系统中的无条件稳定性,但是缺乏算法在非线性系统中稳定性和精度的研究。本文对比分析Chang法(CEM)、CR法(CRM)及RST法在线性系统中的精度问题及在非线性系统中的数值特性,以寻求更适用于实时子结构试验的数值积分算法。

    • 当数值子结构和试验子结构均为单自由度(SDOF)体系时,结构运动微分方程可以写成离散形式为:

      $${m_N}{\ddot u_{i + 1}} + {c_N}{\dot u_{i + 1}} + {k_N}{u_{i + 1}} + {R_E}\left( {{u_{i + 1}},{{\dot u}_{i + 1}}} \right) = {f_{i + 1}}$$ (1)

      式中: ${m_N}$ ${c_N}$ ${k_N}$ 分别为数值子结构的质量、阻尼和刚度; ${\ddot u_{i + 1}}$ ${\dot u_{i + 1}}$ ${u_{i + 1}}$ 分别为数值子结构在第i+1步的加速度、速度和位移; ${R_E}\left( {{u_{i + 1}},{{\dot u}_{i + 1}}} \right)$ 为试验子结构在第i+1步的反作用力; ${f_{i + 1}}$ 为第i+1步的外部激励。将位移和速度表达式统一写成:

      $$\left\{ \begin{array}{l} {u_{i + 1}} = {u_i} + {\alpha _1}\Delta t{{\dot u}_i} + {\alpha _2}\Delta {t^2}{{\ddot u}_i} \\ {{\dot u}_{i + 1}} = {{\dot u}_i} + {\alpha _3}\Delta t{{\ddot u}_i}{\rm{ + }}{\alpha _{\rm{4}}}\Delta t{{\ddot u}_{i{\rm{ + 1}}}} \end{array} \right.$$ (2)

      式中: $\Delta t$ 为时间步长;基于模型的显式积分算法中系数 ${\alpha _1}$ ${\alpha _2}$ ${\alpha _3}$ ${\alpha _{\rm{4}}}$ 分别为:

      $$\left\{ {\begin{array}{*{20}{l}} {{\alpha _1} = \dfrac{{4 + 4\xi \Omega }}{{4 + 4\xi \Omega + {\Omega ^2}}},}&{{\alpha _2} = \dfrac{2}{{4 + 4\xi \Omega + {\Omega ^2}}},}\\ {{\alpha _3} = \dfrac{1}{2},}&{{\alpha _4} = \dfrac{1}{2}}&\!\!{\left( {{\rm{CEM}}} \right)} \\ {{\alpha _1}{\rm{ = 1}},}&{{\alpha _2} = \dfrac{4}{{4 + 4\xi \Omega + {\Omega ^2}}},}\\ {{\alpha _3} = \dfrac{4}{{4 + 4\xi \Omega + {\Omega ^2}}},}&{{\alpha _4} = 0}&\!\!{\left( {{\rm{CRM}}} \right)} \\ {{\alpha _1} = \dfrac{4}{{4 + 4\xi \Omega + {\Omega ^2}}},}&{{\alpha _2} = \dfrac{{4 - 2\xi \Omega - 8{\xi ^2}}}{{4 + 4\xi \Omega + {\Omega ^2}}},}\\ {{\alpha _3} = 1,}&{{\alpha _4} = 0}&\!\!{\left( {{\rm{RST}}} \right)} \end{array}} \right.$$ (3)

      式中: $\Omega \!\! = \!\! \Delta t{\omega _n} \! = \! \Delta t\sqrt {\left( {{k_N} + {k_E}} \right)/m}$ ; ${\xi _N} \! = \! {c_N}/\left( {2m{\omega _n}} \right)$ ; ${\xi _E} = {c_E}/\left( {2m{\omega _n}} \right)$ $\xi {\rm{ = }}{\xi _N} + {\xi _E}$

    • 多自由度(MDOF)体系的结构运动微分方程可以写成矩阵的形式:

      $${{M}} \cdot {\ddot{ U}}\left( t \right) + {{C}} \cdot {\dot{ U}}\left( t \right) + {{K}} \cdot {{U}}\left( t \right) = {{F}}\left( t \right)$$ (4)

      式 中:MCK分别为多自由度体系的质量矩阵、阻尼矩阵和刚度矩阵; ${\ddot{ U}}$ ${\dot{ U}}$ ${{U}}$ 分别为加速度向量、速度向量和位移向量; ${{F}}$ 为力向量。为简单起见,这里假定矩阵C为比例阻尼矩阵。

      对于一个MDOF体系,各自由度之间的动力反应是相互耦合的,因此需要通过转换到模态坐标系进行解耦,考虑到模态的正交性,将式(2)和式(4)改写为:

      $${{{M}}^*} \cdot {\ddot{ Y}}\left( t \right) + {{{C}}^{\rm{*}}} \cdot {\dot{ Y}}\left( t \right) + {{{K}}^{\rm{*}}} \cdot {{Y}}\left( t \right) = {{{F}}^{\rm{*}}}\left( t \right)$$ (5)
      $$\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left\{ \begin{array}{l} {{{Y}}_{i + 1}} = {{{Y}}_i} + {{\alpha }}_1^{\rm{*}}\Delta t{{{\dot{ Y}}}_i} + {{\alpha }}_2^{\rm{*}}\Delta {t^2}{{{\ddot{ Y}}}_i} \\ {{{\dot{ Y}}}_{i + 1}} = {{{\dot{ Y}}}_i} + {{\alpha }}_3^{\rm{*}}\Delta {t^2}{{{\ddot{ Y}}}_i} \\ \end{array} \right.$$ (6)

      式(5)和式(6)均代表N个解耦的联立方程,其中 ${{{M}}^*} = {{{\varPhi }}^{\rm{T}} }{{M\varPhi }}$ ${{{C}}^*} = {{{\varPhi }}^{\rm{T}} }{{C\varPhi }}$ ${{{K}}^{\rm{*}}} = {{{\varPhi }}^{\rm{T}} }{{K\varPhi }}$ 分别为模态质量矩阵、模态阻尼矩阵和模态刚度矩阵,均为对角阵; ${{{F}}^*}\left( t \right) = {{{\varPhi }}^{\rm{T}} } \cdot {{F}}\left( t \right)$ 为模态力向量; ${{\varPhi }}$ N×N阶振型矩阵, ${{\varPhi }}{\rm{ = }}$ [ $\begin{array}{*{20}{c}} {{\phi _1}}&{{\phi _2}}& \cdots &{{\phi _N}} \end{array}$ ] ( $i =1, 2, \cdots, N $ )为第i阶振型向量,N为模态数; ${{Y}}\left( t \right)$ 为模态坐标向量,其与几何坐标向量 ${{U}}\left( t \right)$ 的关系为 ${{U}}\left( t \right) = {{\varPhi }} \cdot {{Y}}\left( t \right)$ ${{\alpha }}_1^{\rm{*}}{\rm{ = }}{{{\varPhi }}^{ - 1}}{\alpha _1}{{\varPhi }}$ ${{\alpha }}_2^{\rm{*}}{\rm{ = }}{{{\varPhi }}^{ - 1}}{\alpha _2}{{\varPhi }}$ ${{\alpha }}_3^{\rm{*}}{\rm{ = }}{{{\varPhi }}^{ - 1}}{\alpha _3}{{\varPhi }}$ 为对角积分参数矩阵。

      MDOF体系在模态坐标系下的动力反应是解耦的,因此对应于每阶模态的运动方程求解都可以采用同SDOF体系完全相同的步骤进行,然后再利用两个坐标系下各矩阵关系,最终得到几何坐标系下的解。

    • 研究积分算法的数值特性主要从稳定性和精度两个方面着手。首先将式(1)~式(3)写成递推矩阵的形式:

      $${{{X}}_{i + 1}} = {{A}}{{{X}}_i} + {{L}}{F_{i + 1}}$$ (7)

      式中: ${{{X}}_i} = {[{u_i},\Delta t{\dot u_i},\Delta {t^2}{\ddot u_i}]^{\rm{T}}}$ ${{L}}{\rm{ = [0,\;0,\;}}{{\Delta {t^2}} / m}{{\rm{]}}^{\rm{T}}}$ A为放大矩阵。

      以SDOF无阻尼系统为模型考核算法稳定性是有意义的,文献[4]中已经验证了各算法在线性系统中的无条件稳定性,因此不再赘述。此处主要讨论算法的精度问题,通常采用数值损耗和飘移来衡量,这两项指标是由算法的数值阻尼引起的,具体表现在计算中会出现振幅衰减(或增长)和周期延长(或缩短)。当算法无条件稳定时, ${\lambda _{1,2}}$ 是一对共轭复根,记作:

      $${\lambda _{1,2}}{\rm{ = }}\sigma \pm {\rm{i}}\varepsilon = \exp ( { - \bar \xi \cdot \bar \Omega \pm {\rm{i}}{{\bar \Omega }_d}} )$$ (8)

      式中: ${\rm{i}} = \sqrt { - 1} $ $\sigma = {A_1}$ $\varepsilon {\rm{ = }}\sqrt {{A_2} - A_1^2} $ $ {\bar \Omega _d}{\rm{ = }}{\bar \omega _d}\Delta t = $ $ {\rm{arctan}}\left( {\left| {{\varepsilon / \sigma }} \right|} \right)$ $ \bar \Omega = {\bar \omega _n}\Delta t = {{{{\bar \Omega }_d}} / {\sqrt {1 - {{\bar \xi }^2}} }}$ $\bar \xi = - {{\ln ( {{\sigma ^2} + } }}$ ${{\varepsilon ^2})/ {2\bar \Omega }}$ 为数值阻尼比。

      定义数值阻尼比与结构阻尼比的差值为算法阻尼比KSI,周期延长率PE为:

      $$ {\rm{KSI = }}\bar \xi - \xi , \;\; {{\rm{PE}} = \dfrac{{\bar T - T}}{T} = \dfrac{\omega }{{\bar \omega }} - 1} $$ (9)

      在线性系统中,三种算法具有相似的数值特性,图1给出了三种算法在结构阻尼比为0.1时,算法阻尼比随Ω的变化曲线,可见当结构阻尼比为0时,各算法的数值阻尼比始终为0,说明此时三种算法均无数值阻尼;当结构阻尼比不为零时,算法阻尼比的绝对值均随Ω的增加而增大,CR法和RST法的算法阻尼比重合,且绝对值小于Chang法。图2显示,三种算法的周期延长率在线性条件下保持一致,且其值随Ω的增加而增大,随ξ的增大而减小。

      图  1  算法阻尼比随Ω的变化曲线

      Figure 1.  Variation of numerical damping ratio with Ω

      图  2  周期延长率随Ω的变化曲线

      Figure 2.  Variation of period errors with Ω

    • 为了描述积分算法在非线性系统中的数值特性,需要定义结构在第i+1步结束时的瞬时刚度为 ${k_{i + 1}}$ ,其与结构初始刚度 ${k_0}$ 的比值称为瞬时非线性度 ${\delta _{i + 1}}$ 。当 ${\delta _{i + 1}}{\rm{ = }}1$ 时,表示第i+1步结束时的瞬时刚度等于初始刚度; $0 < {\delta _{i + 1}} < 1$ 时,表示第i+1步结束时发生了瞬时刚度软化;而 ${\delta _{i + 1}} > 1$ 时,表示第i+1步结束时发生了瞬时刚度强化。由稳定条件和放大矩阵的特征值可以得到各算法在特定瞬时非线性度 ${\delta _{i + 1}}$ 下的稳定性,结果可知三种算法对于刚度软化系统是无条件稳定的,对于刚度硬化系统有条件稳定,其中各算法在 ${\delta _{i + 1}}{\rm{ = }}2.0$ 时谱半径随Ω的变化曲线如图3所示。图3中可见,在同一ξ条件下Chang法的稳定界限最大,其次是RST法;Chang法和RST法的稳定界限随ξ的增加而增大,CR法的稳定界面与ξ无关,始终保持为2。图4图5分别给出了在ξ=0.1且δ分别等于0.5和2.0条件下,三种算法的算法阻尼比和周期延长率随Ω的变化规律。

      图  3  谱半径随Ω的变化曲线(δi+1>1)

      Figure 3.  Variation of upper stability limit with Ω (δ i+1>1)

      图  4  算法阻尼比随Ω的变化曲线(ξ=0.1)

      Figure 4.  Variation of numerical damping ratio with Ω (ξ=0.1)

      图  5  周期延长率随Ω的变化曲线(ξ=0.1)

      Figure 5.  Variation of period errors limit with Ω (ξ=0.1)

    • 本算例为一个三层隔振体系,其中主体结构为剪切型框架(数值子结构),底层为隔振层并设有一个阻尼器(试验子结构)。数值子结构的各层质量 $m = 160 \times {10^3}\;{\rm{kg}}$ ,刚度 ${k_N} = 360 \times {10^6}\;{\rm{N/m}}$ ,阻尼 ${c_N} = 0\;{\rm{N}} \cdot {\rm{s/m}}$ ;试验子结构的刚度 ${k_E} = 9 \times $ 106 N/m,阻尼 ${c_E} = 30 \times {10^5}\;{\rm{N}} \cdot {\rm{s/m}}$ 。初始条件为 ${u_0} $ =1 cm, ${\dot u_0} = 0\;{\rm{cm/s}}$ 。使用三种算法对结构进行自由振动下的结反应计算,积分时间间隔Δt=0.02 s,Newmark显式算法(Δt=0.001 s)的计算结果作为准确解用于对比。

      图6给出了结构作自由振动时顶层位移、速度和加速度的时程反应曲线。图6中可见,RST法和Chang法的精度优于CR法;当t=0时刻,由各数值积分算法得到的速度反应均小于原体系的速度,说明数值积分算法会改变原体系初始条件;在0.1 s内由各数值积分算法得到的顶层加速度反应与准确解的误差较大;0.1 s后与准确解吻合较好。

    • 本算例为一个具有刚度软化特性的两层剪切型框架结构,将底层框架作为试验子结构,而二层框架作为数值子结构。各层质量分别为 ${m_E} = {\rm{3}} \times $ $ {10^{\rm{3}}}\;{\rm{kg}}$ ${m_N} = {\rm{2}}{\rm{.5}} \times {10^{\rm{2}}}\;{\rm{kg}}$ ;各层初始刚度分别为 ${k_{E,{\rm{0}}}} = {10^{\rm{7}}}\;{\rm{N/m}}$ ${k_{N,0}} = {10^4}\;{\rm{N/m}}$ ;第i层的刚度表达式为 ${k_i} = {k_{i,0}}( {1 - \sqrt {\left| {\Delta {u_i}} \right|} } )$ ,其中 ${k_{i,0}}$ 为第i层的初始刚度, $\Delta {u_i}$ 为第i层的层间位移;各层阻尼分别为 ${c_E} = {c_N} = 0\;{{{\rm{N}} \cdot {\rm{s}}} / {\rm{m}}}$ 。计算可得前两阶初始频率分别为 ${\omega _1} = 6.32\;{{\rm{rad}} / {\rm{s}}}$ ${\omega _2} = 57.76\;{{\rm{rad}} / {\rm{s}}}$ 。由三种算法求解结构的动力反应,积分时间步长Δt=0.02;外荷载为 $\sin \left( {2t} \right)$ ,初始条件为 ${u_0} = 0\;{\rm{cm}}$ ${\dot u_0} = 0\;{\rm{cm/s}}$ 。Newmark显式算法(Δt=0.001 s)的计算结果作为准确解用于对比。图7给出了三种积分算法计算的结构顶层反应时程曲线,可以得到与精度分析相同的结论。

      图  6  三层隔振体系自由振动的顶层结构反应

      Figure 6.  Free vibration responses of the top story of a 3-story shear building with isolation

      图  7  两层剪切型框架结构受迫振动的顶层结构反应

      Figure 7.  Forced vibration responses of the top story of a 2-story shear building

    • 数值积分算法是求解结构运动方程的有效方法,也是实现实时子结构试验的关键技术之一。本文选取三种基于模型的显式积分算法:Chang法、CR法和RST法,对比分析了各算法在线性系统和非线性系统中的稳定性和精度。稳定性分析可知,三种算法在线性系统和具有刚度软化特性的非线性系统中均是无条件稳定的,在具有刚度硬化特性的非线性系统中变为有条件稳定,且Chang法和RST法的稳定界限随Ω的增加而增大,CR法的稳定界限始终为2。精度分析可知,当结构阻尼比为零时,三种算法均无数值阻尼;当结构阻尼比不为零时,三种算法均存在数值阻尼,并且在同一个Ω下CR法的数值阻尼比的绝对值最小。三种算法的周期延长率均随Ω的增加而增大,并且在同一个Ω下RST法的周期延长率最小。两个算例表明,RST法和Chang法的精度要优于CR法,但因Chang法是半显式的,因此RST法更适于实时子结构试验中数值子结构的仿真计算。

参考文献 (12)

目录

    /

    返回文章
    返回