基于回传射线矩阵法的非均匀ig梁多轴响应

基于回传射线矩阵法的非均匀ig梁多轴响应

回传波束矩阵法是一种新型发展趋势的分析方法，用于解决空间和时间变量的线性偏差方程和复杂边界条件。在分析平面结构的瞬态响应时，ho等人提出了该方法，并在结构和层的瞬态响应分析中成功应用于结构和层的瞬态响应分析。tian等人讨论了梁缝的检测问题，并进一步介绍了该方法在实际工程中的应用价值。pao等人进一步改进了回传波束矩阵法，并与传统的传递矩阵法进行了详细比较，在估计静态和初始瞬态响应的情况下，决定了其准确性、数值稳定性和列式统一等优点。这是一种理想的弹性动力学分析。回传射线矩阵法的主要思路是:采用连续单元模型和一套对偶局部坐标系统,应用积分变换手段将控制偏微分方程组转换为变换域内的常微分方程组;以波幅为基本未知量,并根据物理意义将弹性波分成离开波(出射波)和到达波(入射波),引入节点处力平衡条件和位移协调条件,建立散射关系,从数学上清晰地表示出节点处的出射波和入射波波幅之间的关系;联合根据对偶局部坐标系的特点建立的相位关系,求解波幅未知量;最后利用积分逆变换方法得到所关心的真实物理空间的解.迄今为止,回传射线矩阵法列式仅能考虑有限个集中力荷载作用情况,无法处理其他荷载形式,如分布荷载和移动荷载.梁结构在分布荷载作用下的动力响应是实际工程中经常出现的问题,具有很强的现实意义.Timoshenko等人以Euler梁为对象,利用模态叠加法求解了均匀简支梁在不同分布荷载下的动力响应.Tong等人基于传递矩阵法对简谐分布荷载作用下的非均匀Timoshenko梁进行了分析.Doyle提出了频谱单元法,将分布荷载等效为节点荷载.本文从回传射线矩阵法的基本思路出发,以非均匀Timoshenko梁为例,对任意分布荷载作用下结构的动力响应进行分析.本文不仅在理论上对回传射线矩阵法进行了深化,而且拓展了其应用范围,使其能更好地为实际工程服务.1d4vd4+2密度7结构参数考虑一根受竖向分布荷载q(x,t)作用的非均匀梁,长度为L,如图1(a)所示.梁的材料密度ρ、弹性模量E和横截面面积A可沿梁的长度方向任意变化.对于这样一根非均匀梁,直接求解往往非常困难,很多情况下难以求得解析解.为此,首先将梁沿长度方向划分为N段,当每段长度足够小时,可近似认为各梁段内的材料参数均为常数,由此就建立了分段均匀梁模型,如图1(b)所示.视两端边界或相邻梁段的交界点为节点,每个梁段用其两端的节点编号表示,与该梁段相关的各个物理量采用节点编号作上标以示区别,如lIJ表示梁段IJ的长度.在回传射线矩阵法中,对每个梁段IJ引入一对对偶局部坐标系(xIJ,yIJ)和(xJI,yJI),坐标原点分别位于两端节点处,如图1(c)所示.采用Timoshenko梁理论,考虑剪切和转动惯量的影响,局部坐标系下每个梁段(忽略上标IJ或JI)的运动方程为EΙz∂2ϕ∂x2+κAG(∂v∂x-ϕ)=ρΙz∂2ϕ∂t2,κAG(∂2v∂x2-∂ϕ∂x)=ρA∂2v∂t2-q(x,t).}(1)式中:E为弹性模量,G为剪切模量,ρ为密度,A为横截面面积,Iz为横截面惯性矩,κ为剪切系数,v为挠度,ϕ为转角.假设L0、A0、E0、ρ0为相应的长度、面积、模量和密度的参考量,令ξ=x/L0,τ=t/t0,ˉv=v/L0,ˉq(ξ,τ)=qL0/(E0A0),其中t0=L0/c1,c1=√E0/ρ0为纵波波速,则方程(1)可无量纲化为∂2ˉϕ∂ξ2+mn[∂ˉv∂ξ-ˉϕ]=mEmρ∂2ˉϕ∂τ2,∂2ˉv∂ξ2-∂ˉϕ∂ξ=nmEmρ∂2ˉv∂τ2-nmAmEˉq.}(2)式中:符号上一横表示无量纲量;n=E/(κG)为弹性模量与修正剪切模量之比;m=AL20/I是与横截面尺寸有关的无量纲系数;mΑ=A0/A,mΕ=E0/E,mρ=ρ0/ρ分别梁段面积、弹性模量、密度的比值.假设梁受到简谐荷载ˉq(ξ,τ)=ˉq0(ξ,ˉω)eiˉωτ的作用,式中ˉω=ωt0为无量纲频率.不失一般性,考虑梁段材料参数相同(mΕ=mρ=1)情形.设ˉv=V(ξ)eiˉωτ,ˉϕ=Φ(ξ)eiˉωτ,方程(2)可简化为d2Φdξ2+mn(dVdξ-Φ)=-ˉω2Φ,d2Vdξ2-dΦdξ=-nˉω2V-nmAˉq0(ξ).}(3)从式(3)中消去变量Φ可得d4Vdξ4+ˉω2(n+1)d2Vdξ2+(nˉω2-m)ˉω2V=(m-nˉω2)mAˉq0-nmAd2ˉq0dξ2.易得上式的齐次解为V(ξ)=a2eik2ξ+d2e-ik2ξ+a3eik3ξ+d3e-ik3ξ.(4)式中:a2、a3为待定入射波波幅;d2、d3为待定出射波波幅;k2、k3为波数,按下式确定:k2,3=√[(n+1)ˉω2±√(n+1)2ˉω4-4ˉω2(nˉω2-m)]/2.(5)同理可求得Φ(ξ)=g2a2eik2ξ-g2d2e-ik2ξ+g3a3eik3ξ-g3d3e-ik3ξ.(6)其中:g2,3=imk2,3m+n(k22,3-ˉω2).(7)相应于分布荷载的非齐次解可采用Fourier变换求得V(ξ)=-A2k2∫ξ0ˉq0(s)sin(k2(ξ-s))ds+A3k3∫ξ0ˉq0(s)sin(k3(ξ-s))ds,(8)Φ(ξ)=-A4∫ξ00(s)cos(k2(ξ-s))ds+A4∫ξ00(s)cos(k3(ξ-s))ds.(9)式中:A2,3=(m-nˉω2+nk22,3)mAk22-k23,A4=m[JX-*2/7].[JX*2/7]mAk22-k23.最后可得Timoshenko梁在任意简谐均布荷载作用下位移ˉv(ξ,τ)和转角ˉϕ(ξ,τ)的通解为ˉv(ξ,τ)=[a2eik2ξ+d2e-ik2ξ+a3eik3ξ+d3e-ik3ξ-A2k2∫ξ0ˉq0(s)sin(k2(ξ-s))ds+A3k3∫ξ0ˉq0(s)sin(k3(ξ-s))ds]eiˉωτ,(10)ˉϕ(ξ,τ)=[g2a2eik2ξ-g2d2e-ik2ξ+g3a3eik3ξ-g3d3e-ik3ξ-A4∫ξ0ˉq0(s)cos(k2(ξ-s))ds+A4∫ξ0ˉq0(s)cos(k3(ξ-s))ds]eiˉωτ.(11)梁内的剪力和弯矩可相应表示为Μ(ξ,τ)=EΙzL0⋅dˉϕ(ξ,τ)dξ,Q(ξ,τ)=κAG⋅[dˉv(ξ,τ)dξ-ˉϕ(ξ,τ)].}(12)1.1局部散射的方程下面根据回传射线矩阵法的基本思路对任意简谐分布荷载作用下的非均匀Timoshenko梁进行分析.首先,在对偶局部坐标系下的所有节点建立力的平衡和位移协调关系;以梁单元IJ和JK的共同节点J为例,有ˉvJΙ(0,τ)=-ˉvJΚ(0,τ),ˉϕJΙ(0,τ)=ˉϕJΚ(0,τ),ΜJΙ(0,τ)+ΜJΚ(0,τ)=0,QJΙ(0,τ)-QJΚ(0,τ)=0.}(13)令式(10)～(12)中ξ=0并将其代入上式可得dJ=SJaJ.(14)式中:dJ=[dJΙ2,dJΙ3,dJΚ2,dJΚ3]T为入射波波幅向量;aJ=[aJΙ2,aJΙ3,aJΚ2,aJΚ3]T为出射波波幅向量;SJ为4×4的系数矩阵,称为节点J处的局部散射矩阵,其元素不难由式(13)得出.在两端边界节点处(节点0和N),局部散射矩阵的维数将变为2×2,由相应的边界条件得到d0=S0a0,dN=SNaN.(15)式中:d0=[d012,d013]T;dN=[dΝ(Ν-1)2,dΝ(Ν-1)3]T;a0=[a012,a013]T;aN=[aΝ(Ν-1)2,aΝ(Ν-1)3]T.将式(14)所示的局部散射关系按节点顺序排列组集,可得到分段均匀梁的整体散射关系为d=Sa.(16)式中:d4N×1=[(d0)T,(d1)T,…,(dN-1)T,(dN)T]T和a4N×1=[(a0)T,(a1)T,…,(aN-1)T,(aN)T]T分别为整体出射波和入射波波幅向量;S=[S002×4⋯02×402×204×2S1⋯04×404×2⋮⋮⋱SΝ-104×202×202×4⋯02×4SΝ]为整体散射矩阵.向量d和a一共含有8N个未知量,而式(16)只提供4N个方程.为求得所有的未知量,须再补充一组方程.1.2相位矩阵的建立同一梁段某点处的物理量可以在一对对偶局部坐标系下同时表示,因此两个表达式必须等价.以位移为例,在同一梁段IJ内某点处存在以下关系:IJ(ξIJ,τ)=-JI(IJ-ξIJ,τ).(17)其中IJ为梁段IJ的无量纲长度.对同一梁单元,易知kΙJi=kJΙi,AΙJi=AJΙi(i=2,3),ΙJ0=-JΙ0.经过一定的数学推导,最后可得到以下关系:aΙJ2=-dJΙ2e-ikΙJ2ˉlΙJ+qΙJ2,aΙJ3=-dJΙ3e-ikΙJ3ˉlΙJ+qΙJ3,aJΙ2=-dΙJ2e-ikΙJ2ˉlΙJ+qJΙ2,aJΙ3=-dΙJ3e-ikΙJ3ˉlΙJ+qJΙ3.}(18)式中:qΙJ2,3=±AΙJ2,32ikΙJ2,3∫ˉlΙJ0ˉqΙJ0(s)e-ikΙJ2,3ˉlΙJds,qJΙ2,3=∓AΙJ2,3e-ikΙJ2,3ˉlΙJ2ikΙJ2,3∫ˉlΙJ0ˉqΙJ0(s)eikΙJ2,3ˉlΙJds.}(19)由式(18)可以得到如下相位关系:a=PUd+q.(20)式中:P为整体相位矩阵(维数4N×4N),定义为Ρ=[Ρ0104×4⋯04×404×4Ρ12⋯04×4⋮⋮⋱⋮04×404×4⋯Ρ(Ν-1)Ν],(21)ΡΙJ=[-e-ik2ΙJl¯ΙJ0000-e-ik3ΙJl¯ΙJ0000-e-ik2ΙJl¯ΙJ0000-e-ik3ΙJl¯ΙJ],(22)U是整体置换矩阵,表示为U=[U004×4⋯04×404×4U0⋯04×4⋮⋮⋱⋮04×404×4⋯U0],(23)其中U0={{0,0,1,0},{0,0,1,0},{1,0,0,0},{0,1,0,0}};q=[q201,q301,q210,q310,…,q2(Ν-1)Ν,q3(Ν-1)Ν,q2Ν(Ν-1),q3Ν(Ν-1)]T,与均布荷载有关.将式(20)代入式(16)有d=(I-R)-1s.(24)式中:R=SPU称为回传射线矩阵(thereverberation-raymatrix),s=Sq称为整体源向量,I为单位矩阵.将求解出的d的表达式回代入式(20)或者式(16),即可求得a.由此未知波幅向量a和d完全确定,将其代入位移、转角或者内力的表达式中,即可求得梁上任意一点的各物理量随时间的变化.可以看出,与集中力作用下的回传射线矩阵法列式不同的是,本文中与外荷载有关的源向量出现在相位关系里面,而不是出现在散射关系中.当谐荷载为零,即s=0时,为得到d的非平凡解,必须满足|Ι-R|=0.(25)上式即为分段均匀梁的自由振动频率方程.若梁上施加非简谐荷载,则运用傅里叶变换:g(τ)=12π∫-∞∞g^(ω)eiωτdω,g^(ω)=∫-∞∞g(τ)e-iωτdτ.}(26)易知在材料参数不变的情况下,式(2)经变换后形式与式(3)相同,只需将式(3)中的V(ξ)、Φ(ξ)、0(ξ)分别替换为v^(ξ,ω)、ϕ^(ξ,ω)、q^(ξ,ω)即可.前面的推导步骤可保持不变,但最后给出的是频域中的表达式.为得到时域解,需要将响应作傅里叶逆变换,以挠度为例,有v¯(ξ,τ)=12π∫-∞∞v^(ξ,ω)eiωτdτ.(27)在v^(ξ,ω)中包含了入射波波幅和出射波波幅,其中d须由式(24)式确定.但是当ω=i时(i为第i阶无量纲自振频率),|Ι-R|=0,不能直接求逆.为解决这个问题,在回传射线矩阵法中采用下述Neumann级数展开(I-R)-1=I+R+R2+…+RN+….(28)这不仅消除了瞬态响应积分计算中出现的极点,而且展开式中的每一项都有明确的物理意义,对应不同的波动过程.关于该级数的收敛性以及展开项数的讨论详见文献.2材料及边界条件考虑一根长为3m的矩形简支梁,截面高0.1m,宽0.05m,梁的弹性模量为206GPa,剪切模量为79GPa,密度7850kg/m3,剪切系数κ取为π2/12.若梁上施加有大小为2000N/m、频率为102rad/s的简谐均布荷载,用本文方法计算梁的稳态响应.计算时长度单位用L0=1m无量纲化,其余参考量皆取为梁自身的材料常数.图2为本文方法与精确解在0.25L¯处的竖向位移比较曲线,可看出两者结果是一致的.3非均匀简支梁瞬态响应曲线算例1非均匀简支梁的稳态响应.无量纲参考量及材料性质均与上例同.假设梁的无量纲长、宽分别为L¯=1、b¯=0.1;截面高度沿梁的长度方向呈指数变化h¯=0.2e-ξ/4;梁上作用有无量纲简谐分布荷载q¯=1×10-3eiτ/50.在计算前先将梁分段均匀化,此处考虑了3种情况,即梁的划分段数分别为N=5,10,20.图3、4分别显示了用本文方法计算所得的L¯/4、L¯/2处挠度响应和截面最大弯曲应变响应曲线,可以看出随段数N值的增加计算结果很快收敛.算例2非均匀简支梁的瞬态响应.假设作用在Timoshenko梁上的荷载为如下矩形脉冲荷载:q¯(x,t)={1×10-6,0≤t≤6t0;0,t>6t0.其他参数与算例1相同.在用回传射线矩阵法进行计算时,梁的划分段数取为20,Neumann级数展开项数取为80,计算所得的瞬态响应如图5、6所示,其中图6中的h1、h2分别为L¯/40和L¯/4处的截面高度.为验证结果的正确性,将其与Ansys软件的计算结果进行了对比,其中位移符合良好.需指出的是,在用Ansys软件计算瞬态响应时,梁单元必须划