第 34 卷 第 5 期 204 年 5 月 北京理工大学学报 TransactionsofBeijingInstituteofTechnology Vol.34 No.5 May204 基于结构势能原理的动态载荷等效静态转化方法 张艳岗, 苏铁熊, 李伟环 2, 毛虎平, 郭支明, 李建军 3 (. 中北大学机械与动力工程学院, 山西, 太原 03005;2. 北方信息控制集团, 山西, 太原 030006; 3. 中北大学理学院, 山西, 太原 03005) 摘要 : 针对结构动态响应优化过程中动态分析的复杂性和高耗时性, 以及直接进行动态优化的不可行性等问题, 提出了一种基于结构势能的动态载荷等效静态转化方法. 该方法通过分析结构势能的本构方程, 应用动态应力解空间谱元离散的关键点识别方法得到了结构最危险的时刻, 然后建立了基于结构势能原理的动态载荷等效静态转化数学模型, 应用区域细分算法求解, 最终实现了动态载荷等效静态转化, 并以 0 杆桁架结构 24 杆桁架结构以及两层桁架结构为例进行验证. 结果表明 : 文中方法得到的位移比文献方法更接近实际位移的 6.74%,4.57%, 8.50%, 说明了文中方法的可行性和有效性. 关键词 : 结构势能 ; 动态载荷 ; 等效静态转化 ; 关键点中图分类号 :TH22 文献标志码 :A 文章编号 :00-0645(204)05-0454-06 Transformation MethodofEquivalentStaticLoadsfrom Dynamic LoadsBasedonStructuralPotentialPrinciple ZHANG Yaṉgang, SU Tie-xiong, LIWei-huan 2, MAO Hu-ping, GUOZhi-ming, LIJiaṉjun 3 (.ColegeofMechanicalandElectricalEngineering,NorthUniversityofChina,Taiyuan,Shanxi03005,China; 2.NorthernofChinaInformationControlGroupCo.,Ltd.Taiyuan,Shanxi030006,China; 3.SchoolofScienceNorthUniversityofChina,Taiyuan,Shanxi03005,China) Abstract:Duetotheproblem ofcomplexityandhigẖcostcharacterofdynamicanalysisinthe processingofstructuraldynamicoptimization,andinfeasibilitytodirectdynamicoptimization,a transformationmethodofequivalentstaticloadsfrom dynamicloadsbasedonvolumepotential energywasproposed.analyzedtheconstitutiveequationofvolumepotentialenergy,themost dangeroustimepointsofthestructuredynamicresponsewerefoundbythecriticaltimepoints identification method for solution space of global dynamic stress based spectral element interpolation.a mathematicalmodelforthetransformationmethodofequivalentstaticloadswas constructedanditwassolvedbydirectalgorithm.last,0bartruss,24-memberplanetruss and-bay2-storeyframe weretakenasexamplestoverifythefeasibilityandvalidityofthe proposedmethod.resultsshowthatthedisplacementusedproposed methodis morerealistic thantheliteraturedisplacementof6.74%,4.57%,8.5%,indicatingthattheproposed method isfeasibleandefective,andlayafoundationforthestructuraldynamicresponseoptimization. Key words:structural potential;dynamicloads;transformation of equivalent staticloads; criticalpoints 收稿日期 :203 09 02 基金项目 : 国家自然科学基金资助项目 (5275489); 山西省基础研究计划 - 煤层气联合研究基金项目 (20202003); 山西省重点实验室 项目 作者简介 : 张艳岗 (98 ), 男, 博士生, 讲师,E-mail:zyg324@63.com. 通信作者 : 毛虎平 (974 ), 男, 博士, 讲师,E-mail:maohp8@gmail.com.
第 5 期 张艳岗等 : 基于结构势能原理的动态载荷等效静态转化方法 455 在真实的物理世界中, 大多数结构都是在动态 载荷作用下工作的, 由于动态响应计算以及优化的 复杂性与高耗时性 [], 使得结构动态响应优化设计 难以进行. 韩国学者 G.J.Park [26] 已经对动态载荷 等效转化进行总结和发展, 提出基于有限元模态位 移将动态载荷等效为静态载荷的方法, 等效原理是 结构等效前后位移场相同 ; 提出基于应力等效转化 方法, 并且与位移场等效方法进行比较 ; 应用位移等 效静态载荷的等效方法进行了动态响应优化. 在 G.J.Park 文献中, 虽然说明了静态载荷等 效的关键时间点, 即位移最大时间点作为等效载荷 转化的时间点, 然而结构在动态载荷的作用下, 最大 位移点与最大应力点有一个时间差 ; 另外, 结构在动 态载荷的作用下, 除了在动态力作用方向位移比较 大之外, 其他方向都有位移, 如果只关心发生最大位 移的方向而不考虑其他方向, 从一定意义上来说不 合理 ; 如果只考虑应力也不合理. 总之, 采用应力等效存在以下不足 : 最大应 力仅仅说明结构的局部受力状况, 不能反映应力的 分布情况 ;2 由于动态分析采用有限元法计算, 而 结构最大应力与单元类型及单元尺寸有关, 尤其在 应力集中的局部位置, 由于应力变化比较大, 因此单 元类型与尺寸对其影响非常大 ;3 大量计算表明, 一个结构在不同的单元类型与尺寸下获得的最大应 力值差别非常明显. 同理, 采用位移等效方法也存 在类似的问题. 因此本文提出基于结构势能动态载 荷等效方法. 结构势能 结构在动态载荷作用下, 若忽略结构阻尼, 能量 满足如下关系 : 2 U T M U + 2 UT KU + 2 UT SU = F(t)u(t)dt, t () 式中 : 2 U T M U 为动能 ; 2 UT KU + 2 UT SU 为结构 势能 ( 其中 2 UT KU 为体积应变能, 2 UT SU 为应力 刚度能, 也称为外力势能 ; tf(t)u(t)dt 为外力功 ;S 为应力刚度矩阵或应力刚化矩阵 (stressstifness matrix);u 为位移向量 ; U 为速度向量 ;M 为质量矩 阵 ;K 为刚度矩阵 ;F(t) 为动态载荷 ;u(t) 为结构在 力方向的位移. 由于结构线性振动时, 动能比较小, 因此将式 () 近似为 2 UT KU + 2 UT SU F(t)u(t)dt. (2) t 2 关键点识别 Courant [7] 在 943 年的工作就体现了有限元法的思想,960 年 Clough [8] 第一次提出了有限元法 的概念. 有限元法为计算分析结构动力学问题提供 了较为有效的方法, 然而, 结构动态分析过程非常复 杂且耗时较大. 因此, 为有效解决动态响应优化问 题, 考虑用静态载荷等效代替动态载荷, 关键问题是 首先要确定在哪个时刻进行等效, 即在动态载荷作 [9] 用下结构关键点的识别. 关键点是指动态载荷作用下的结构最危险的时 刻. 一直以来, 学术界普遍认为关键点即结构在动 态载荷作用下, 结构位移绝对值最大的时刻, 或者部 分研究中将动态载荷最大的时刻作为关键点. 作者 在前期的研究中发现, 结构在动态载荷作用下位移 达到最大值的时刻与应力达到最大值的时刻有一个 时间差, 另一方面, 动态载荷最大的时刻, 结构位移 与应力都远远没有达到最大值, 因此不论以位移绝 对值最大或者动态载荷最大的时刻作为关键点都不 合适. 因此, 作者提出基于动态应力解空间谱元离 散的关键点识别方法 [0], 并在此基础上进行动态载 荷等效静态转化. 3 等效静态载荷转化 结构在动态载荷作用下, 基于结构势能进行动 态载荷等效静态转化方法如下. 步骤 建立有限元模型. 运动微分方程为 M q+ċq+kq=f, (3) 式中 :M = Φ T mφ,c = Φ T cφ,k = Φ T kφ,f = Φ T f(t),q 为模态位移向量 ; q 为模态加速度向量 ; q 为模态速度向量 ;f(t) 为时间的任意函数 ;Φ 为模态 矩阵 ;m 为质量系数 ;c 为阻尼系数 ;k 为刚度系数. 在本文中假设 c=0. 步骤 2 瞬态动力学分析. 采用模态叠加法. 先求解式 (3), 获得 q, 然后 得到结构的节点位移 x =qφ, 再由公式 σ k (t)= DB kqkφ 获得第 k 阶模态应力, 通过公式 σ i (t)= m σ k(t) i 得到结构动态应力. k=
456 北京理工大学学报第 34 卷 步骤 3 识别关键时间点. 获得结构动态应力解空间, 将其在高斯勒让德 罗巴托点谱元离散, 构造时间点与其对应的动态应 m max ì í s.t. î i= σ i (t k )Φ k+(t) i (t-t i )Φ i ( t k+ k ), 式中 :m 为解空间离散的单元数 ;k 为每一个单元离 散 GLL 点数 ;σ i (t k ) 为第 i 个单元的第 k 个 GLL 点 对应的动态应力 ;t 0 为仿真开始时间, 通常为 0; t 为仿真结束时间. 步骤 4 布 结构势能. 计算关键时间点对应的结构位移分 根据步骤 3 中得到的关键时间点, 通过公式 () 计算结构在这一点的结构势能及位移分布. 模型 式中 : 步骤 5 建立结构动态载荷等效静态转化数学 n min Fi, 2 i= s j d j,j=,2,,m, E g s =E g d, E ge s =E ge d. (5) E g s = 2 UT sk+ 2 UT ssu s ; E g d = 2 UT dk+ 2 UT dsu d ; E e s = 2 UeT s K e + 2 UeT ss e U e s; 力解空间矩阵, 应用 Lagrange 插值, 得到高精度的 近似函数, 最后识别关键时间点. 关键时间点识别 的数学模型 : ìφ i k+(t)=(t-t )(t-t 2 ) (t-t k )(t-t k+ ) t k t (-t 2 ) d k í { 2 k k! dt k[( t 2 -) k ]= 0 }. î t 0 t t (4) 69GPa, 密度为 270kg/m 3, 截面积为.94 0-4 m 2, 泊松比为 0.3, 结构及载荷如图 所示. 图 动态载荷作用的 0 杆桁架结构 Fig. The0bartrussstructurewithdynamicloads 4.. 关键点识别 基于前面的关键点识别方法进行该桁架结构的 时间关键点识别, 图 2 为采用基于应力的谱元离散插 值的关键点识别方法获得的识别结果. 其中有 2 个 关键点, 可以看出 0.033967436002547s 时刻应力 最大, 达到 -98639.2087022Pa, 见表. E e d = 2 UeT d K e + 2 UeT ds e U e d. 式中 :E g s 为静态力作用下整体结构势能 ;E g d 为动态力作用下整体结构势能 ;E e s 为静态力作用下某单元结构势能 ;E e d 为动态力作用下某单元结构势能. 步骤 6 求解. 由于结构等效载荷是未知的, 其特性极为复杂, 一般为非线性多峰值函数, 因此采用全局优化算法 区域细分算法 (DIRECT) 进行搜索结构等效静态载荷. 4 方法应用 4. 0 杆桁架结构本算例的动态载荷为常数, 结构弹性模量为 图 2 0 杆桁架结构关键点识别结果 Fig.2 Criticalpointidentificationof0bartrussstructure 表 0 杆桁架结构关键点识别结果 Tab. Criticalpointidentificationresultsof0 bartrus structure 关键时间点 /s 应力值 /Pa 0.033967436002547-98639.20870220 0.065027948029393 237738.03083364 4..2 关键点等效载荷转化表 2 为本文方法获得 0 杆桁架结构等效静态
第 5 期 张艳岗等 : 基于结构势能原理的动态载荷等效静态转化方法 457 载荷与文献 [2] 方法获得的等效静态载荷. 表 2 等效静态载荷转化 Tab.2 Transformationofequivalentstaticloads 载荷类型 载荷 /N 节点 2 节点 3 节点 6 节点 5 动态 0 0-445.00-445.00 静态 ( 本文 ) 9.02000-58.87-370.37-772.97 静态 ( 文献 [2]) 33.25358-633.29-507.05-535.29 4..3 施加等效载荷计算结果通过本文方法对 0 杆桁架结构等效转化, 获得 了表 3 的结果. 从表 3 中, 可以看出本文方法比文献方法更加合理, 因为本文方法获得的结构节点位移最大相差 2.7%, 而文献 [2] 方法相差 9.440%; 特别是文献 [2] 中没有考虑结构势能的变化, 而本文兼顾了位移与结构势能, 结果很明显, 本文方法获得的结构势能与动态作用下相差.293%, 而文献 [2] 方法相差 3.839%. 载荷类型 表 3 施加等效载荷计算结果对比 Tab.3 Calculationresultscontrastofequivalentload 位移 /m 单元结构势能 /J 节点 2 节点 3 节点 6 节点 5 单元 单元 3 单元 5 单元 0 结构势能 /J 动态 -0.0028326-0.0064780-0.0065344-0.002957.0027 0.0023337.2340 0.00427 4.996 静态 ( 本文 ) -0.0027563-0.0066060-0.006567-0.0029790.0559 0.0058666.609 0.036476 4.2539 静态 ( 文献 [2]) -0.0028326-0.0070894-0.0070340-0.002998.2667 0.0022554.3482 0.08625 4.7808 4.2 24 杆组成的平面桁架为增加对比性, 本算例采用与文献 [2] 中相同的算例. 如图 3 所示, 该桁架结构有 49 个铰链,94 个自由度, 弹性模量 E=207GPa, 泊松比 γ=0.3, 密度 ρ=7850kg/m 3, 杆的截面积 0.645 0-4 m 2. 动态载荷为半正弦函数. 在节点,20,9,8,7, 6,5 的 x 正方向作用同样大的动态载荷, 在节点,2,3,4,5 的 y 负方向也作用同样大的动态载荷. 表 4 为应用基于应力的谱元离散插值的关键点识别方法获得的 4 个关键点, 其中 0.0652572s 时应力最大, 达到 226099.4Pa, 如图 4 所示. 表 4 24 杆桁架结构关键点识别结果 Tab.4 Criticalpointidentification of24-member plane trus 序号 关键时间点 /s σ/pa 0.034390489-523974.8 2 0.0652572 226099.4 3 0.9259952-29976.6 4 0.27604573 4932887.7 表 5 为本文方法获得 24 杆桁架结构等效静态载荷与文献 [2] 方法获得的等效静态载荷. 从其中可以看出, 静态载荷的作用方向除了与动态载荷方向一致之外, 还有与动态载荷方向相反的静态载荷. 这就说明传统的动态因子法不合适, 其仅仅用与动 态载荷方向一致的最大幅值乘以一个系数. 表 6 是应用本文方法获得等效静态载荷, 并与文 献 [2] 的方法进行比较, 可以看出本文方法的优势. 图 3 24 杆平面桁架及动态载荷 Fig.3 The24-memberplanetrussstructurewithdynamicloads
458 北京理工大学学报第 34 卷 Fig.4 图 4 24 杆桁架结构关键点识别结果 Criticalpointidentificationresultsof24-memberplane trussstructure 态载荷作用下最大相差 4.7%, 结构势能相差.20%; 而文献 [2] 方法 x 方向的位移与动态载荷作用下最大相差 6.75%,y 方向的位移与动态载荷作用下最大相差 3.2%, 结构势能相差 9.423%. 4.3 均布力与集中力作用的平面桁架如图 5 所示, 均匀分布力与集中力共同作用的平面桁架结构. 弹性模量 E=207GPa, 密度 ρ= 7850kg/m 3, 竖直杆横截面积为 0.03m 2, 水平杆横截面积为 0.0068m 3. 半正弦均匀分布力 F (t) 作用在水平杆上, 半正弦集中力 F 2 (t) 作用在节点 2 和 3 处. Tab.5 x 方向 y 方向 节点 表 5 24 杆桁架结构等效静态载荷转化 Transformation ofequivalentstaticloadsto24- memberplanetrusstructure 等效静态载荷 /N 本文文献 [2] 7567.4 8039.82 5 8492.2 8598.88 6 8029.30 7734.02 7 5886.48 606.88 8 4040.4 5997.95 9 53.89 3347.6 20 4648.6-270.40 3033.57 5225.4 2 2952.75-70.97 3-6630.83-336.70 4-32.44-5406.72 5-933.84-665.64 表 6 24 杆桁架结构等效静态载荷作用下计算结果对比 Tab.6 Calculationresultscontrastofequivalentloadto24- x 方向 y 方向 memberplanetrusstructure 节点 5 6 7 8 9 20 2 3 4 5 关键时间点的位移 /m 等效静态力作用下的位移 ( 本文 )/m 文献 [2] 分析位移 /m 0.2552 0.2202 0.23007 0.08098 0.0806 0.082 0.073 0.0085 0.0460 0.257 0.2558 0.303 0.4779 0.4926 0.5473 0.730 0.7386 0.8022 0.9368 0.9773 0.20566 0.07269 0.07376 0.07650 0.03689 0.03622 0.036 0.00205 0.00208 0.00209-0.03273-0.03579-0.03886-0.0686-0.0784-0.07765 本文方法及文献 [2] 中结构势能分别为 3593.83,3888.9J. 本文方法 x 方向的位移与动 态载荷作用下最大相差 2.8%,y 方向的位移与动 图 5 均布力与集中力作用的平面桁架 Fig.5 Planetrusssubjecttouniformandconcentratedforce 通过关键识别得到 : 关键时间点为 0.02303560s, 应力为 4905597.27Pa, 见表 7 及 图 6. 表 7 Tab.7 均布力与集中力作用的平面桁架关键点识别结果 Criticalpointidentification resultsofplanetrus subjecttouniformandconcentratedforce 序号关键时间点 /s 单元 0 的应力 /Pa 0.02303560 4905597.27 2 0.043959-48939233.6 图 6 均布力与集中力作用的平面桁架单元 0 的动态应力 Fig.6 Dynamicstressofelement0ofplanetrusssubjectto uniformandconcentratedforce 将动态力作用的点作为静态力的作用点, 这样 就有 44 个静态力, 并将这些点的与静态力作用方向 一致的位移作为约束, 再加上结构势能相等的约束,
第 5 期 张艳岗等 : 基于结构势能原理的动态载荷等效静态转化方法 459 求解结果是等效静态力作用下的结构势能为 464.35J, 动态力作用下为 46.59J, 位移场最大 相差 0.00m. 等效静态载荷为图 7 所示, 从图中进 一步说明了等效静态载荷法比动态因子法更加 合理. 图 7 集中力与分布力共同作用下的等效静态载荷 Fig.7 Transformationofequivalentstaticloadstoplanetruss subjecttouniformandconcentratedforce 5 结论 从结构势能的角度出发, 建立数学模型, 应用区 域细分算法获得等效静态载荷. 如果仅仅从应力 位移角度出发进行等效转化, 无论从应力还是位移 的角度, 都具有一定的片面性, 都没有综合考虑应力 和位移, 基于这一背景下, 本文研究基于结构势能的 结构动态载荷等效静态转化方法, 获得如下结论 : 本文研究的前提是弹性结构及忽略结构阻 尼, 本文方法能够综合体现等效转化时各个物理量 的影响. 结构势能不仅可以反映应力 应变的影响 因素, 而且反映了体积积分的区域, 并反映了位移 弹性模量的影响因素 ; 2 只通过应力等效或者位移等效, 不能全面体 现结构整体的响应, 相反, 通过结构势能等效能包含 动态载荷作用下结构各方面的物理因素, 结构势能 中包含了应力与位移这两种响应 ; 3 静态等效的结果更接近动态载荷作用下的 结构响应, 将本方法应用于结构动态响应优化设计 分析, 可降低动态响应计算的复杂性和高消耗, 同时 有效提高动态响应优化的效率和精度. 参考文献 : []James M L. Vibration of mechanicaland structural systems [M ]. 2nd ed. New York: Harper Colins,994. [2]ChoiW S,ParkGJ.Transformationofdynamicloads intoequivalentstaticloadsbasedon modalanalysis[j]. International Journal for Numerical Methods in Engineering,999,46:29 43. [3]Choi W S, Park K B, Park G J. Calculation of equivalentstaticloadsanditsapplication[j].nuclear EngineeringandDesign,2005,235:2337 2348. [4]KangBS,ChoiW S,ParkGJ.Structuraloptimization underequivalentstaticloadstransformedfrom dynamic loads based on displacement [J]. Computers & Structures,200,79:45 54. [5]Park K J, Lee J N, Park G J.Structural shape optimization using equivalentstaticloadstransformed from dynamic loads [J].International Journal for Numerical Methods in Engineering,2005,63 (4): 589 602. [6]Park G J.Analytic methodsfordesignpractice[m]. New York:Springer,2006. [7]CourantR.Variationalmethodforsolutionsofproblems ofequilibrium andvibrations[j].bulam MathSoc, 943,49: 23. [8]CloughR W.Thefiniteelementmethodinplanestress analysis [C] ASCE Preceeding of 2nd ASCE Conference on Electronic Computation. Pitsburgh, USA:ASCE,960:345 378. [9]GrandhiR V, Haftka R T, Watson L T.Desigṉ oriented identification of critical times in transient response[j].aiaajournal,986,24(4):649 656. [0] 毛虎平, 吴义忠, 陈立平. 基于时间谱元法的动态响应优化 [J]. 机械工程学报,200,46(6):79 87. Mao Huping, Wu Yizhong,Chen Liping.Dynamic responseoptimizationbasedontimespectralelement method[j].journalofmechanicalengineering,200, 46(6):79 87.(inChinese) ( 责任编辑 : 刘雨 )