弹性动力学问题一种新的时空域边界积分方程_姚振汉
- 格式:pdf
- 大小:993.88 KB
- 文档页数:9
冲击动力学冲击动力学分为四章。
第一章包括两章:弹性波和弹塑性波。
第二部分介绍了不同应变率下的动态力学实验技术,总结了高应变率下材料的本构关系。
第三章着重分析了刚塑性梁板的动力响应,第五章介绍了惯性效应和塑性铰,第六章分析了悬臂梁的动力响应,第七章讨论了轴力和剪力对梁动力性能的影响,第八章介绍了模态分析技术、极限定理和刚塑性模型的适用性,第九章介绍了刚塑性板的动力响应分析。
第四章研究材料和结构的能量吸收,其中第10章讨论了材料和结构吸能的一般特征,第11章介绍了典型的吸能结构和材料。
”“碰撞动力学”着重阐述了碰撞动力学的基本概念、基本模型和基本方法。
文中还介绍了动态实验方法以及冲击动力学在冲击防护问题中的应用。
每章附有练习和主要参考文献,供教学和科研参考。
以冲击动力学为教材,可用于40门课程的研究生课程,为固体力学、航空航天、汽车工程、防护工程和国防工程研究生等前沿科学领域的冲击动力学及相关研究方法打下基础。
为他们进行相关的科学研究。
同时,也可供教师、科研人员、工程技术人员和相关专业大四学生自学参考。
作者简介余同希英国剑桥大学哲学博士、科学博士。
曾任北京大学力学系教授、博士生导师;英国曼彻斯特理工大学机械工程系教授。
1995年加入香港科技大学,先后任工学院副院长、机械工程系系主任、协理副校长、霍英东研究院院长等职。
研究主要集中于冲击动力学、塑性力学、结构与材料的能量吸收、复合材料与多胞材料等领域,擅长对工程问题建立力学模型并由此揭示其变形和失效机理。
已发表论文300余篇,担任《国际冲击工程学报》副主编、《国际机械工程学报》副主编,以及十余种学术刊物的编委。
目录绪论第一篇固体中的应力波第1章弹性波1.1 圆杆中的弹性波1.2 弹性波的分类1.3 波的反射和相互作用思考题习题第2章弹塑性波2.1 一维弹塑性波2.2 有限长度杆在高速冲击下的大变形2.2.1 taylor模型2.2.2 用能量法求解taylor杆问题。
第一章.波动方程§1 方程的导出。
定解条件1.细杆(或弹簧)受某种外界原因而产生纵向振动,以u(x,t)表示静止时在x 点处的点在时刻t 离开原来位置的偏移,假设振动过程发生的张力服从虎克定律,试证明),(t x u 满足方程()⎪⎭⎫⎝⎛∂∂∂∂=⎪⎭⎫ ⎝⎛∂∂∂∂x u E x t u x t ρ 其中ρ为杆的密度,E 为杨氏模量。
证:在杆上任取一段,其中两端于静止时的坐标分别为 x 与+x x ∆。
现在计算这段杆在时刻t 的相对伸长。
在时刻t 这段杆两端的坐标分别为:),();,(t x x u x x t x u x ∆++∆++其相对伸长等于 ),()],([)],([t x x u xxt x u x t x x u x x x ∆+=∆∆-+-∆++∆+θ令→∆x ,取极限得在点x 的相对伸长为x u ),(t x 。
由虎克定律,张力),(t x T 等于),()(),(t x u x E t x T x =其中)(x E 是在点x 的杨氏模量。
设杆的横截面面积为),(x S 则作用在杆段),(x x x ∆+两端的力分别为x u x S x E )()(x u x x S x x E t x )()();,(∆+∆+).,(t x x ∆+于是得运动方程 tt u x x s x ⋅∆⋅)()(ρxESu t x =),(x x x x x ESu x x |)(|)(-∆+∆+利用微分中值定理,消去x ∆,再令0→∆x 得u x s x )()(ρx ∂∂=x ESu ()若=)(x s 常量,则得22)(tu x ∂∂ρ=))((xu x E x∂∂∂∂即得所证。
2.在杆纵向振动时,假设(1)端点固定,(2)端点自由,(3)端点固定在弹性支承上,试分别导出这三种情况下所对应的边界条件。
解:(1)杆的两端被固定在l x x ==,0两点则相应的边界条件为.0),(,0),0(==t l u t u(2)若l x =为自由端,则杆在l x =的张力xu x E t l T ∂∂=)(),(|l x =等于零,因此相应的边界条件为xu ∂∂|l x ==0同理,若0=x 为自由端,则相应的边界条件为xu∂∂∣00==x(3)若l x =端固定在弹性支承上,而弹性支承固定于某点,且该点离开原来位置的偏移由函数)(t v 给出,则在l x =端支承的伸长为)(),(t v t l u -。
第9卷㊀第2期2024年3月气体物理PHYSICSOFGASESVol.9㊀No.2Mar.2024㊀㊀DOI:10.19527/j.cnki.2096 ̄1642.1098Mach数和壁面温度对HyTRV边界层转捩的影响章录兴ꎬ㊀王光学ꎬ㊀杜㊀磊ꎬ㊀余发源ꎬ㊀张怀宝(中山大学航空航天学院ꎬ广东深圳518107)EffectsofMachNumberandWallTemperatureonHyTRVBoundaryLayerTransitionZHANGLuxingꎬ㊀WANGGuangxueꎬ㊀DULeiꎬ㊀YUFayuanꎬ㊀ZHANGHuaibao(SchoolofAeronauticsandAstronauticsꎬSunYat ̄senUniversityꎬShenzhen518107ꎬChina)摘㊀要:典型的高超声速飞行器流场存在着复杂的转捩现象ꎬ其对飞行器的性能有着显著的影响ꎮ针对HyTRV这款接近真实高超声速飞行器的升力体模型ꎬ采用数值模拟方法ꎬ研究Mach数和壁面温度对HyTRV转捩的影响规律ꎮ采用课题组自研软件开展数值计算ꎬMach数的范围为3~8ꎬ壁面温度的范围为150~900Kꎮ首先对γ ̄Re~θt转捩模型和SST湍流模型进行了高超声速修正:将压力梯度系数修正㊁高速横流修正引入到γ ̄Re~θt转捩模型ꎬ并对SST湍流模型闭合系数β∗和β进行可压缩修正ꎻ然后开展了网格无关性验证ꎬ通过与实验结果对比ꎬ确认了修正后的数值方法和软件平台ꎻ最终开展Mach数和壁面温度对HyTRV边界层转捩规律的影响研究ꎮ计算结果表明ꎬ转捩区域主要集中在上表面两侧㊁下表面中心线两侧ꎻ增大来流Mach数ꎬ上下表面转捩起始位置均大幅后移ꎬ湍流区大幅缩小ꎬ但仍会存在ꎬ同时上表面层流区摩阻系数不断增大ꎬ下表面湍流区摩阻系数不断减小ꎻ升高壁面温度ꎬ上下表面转捩起始位置先前移ꎬ然后快速后移ꎬ最终湍流区先后几乎消失ꎮ关键词:转捩ꎻHyTRVꎻ摩阻ꎻMach数ꎻ壁面温度㊀㊀㊀收稿日期:2023 ̄12 ̄13ꎻ修回日期:2024 ̄01 ̄02基金项目:国家重大项目(GJXM92579)ꎻ广东省自然科学基金-面上项目(2023A1515010036)ꎻ中山大学中央高校基本科研业务费专项资金(22qntd0705)第一作者简介:章录兴(1998 )㊀男ꎬ硕士ꎬ主要研究方向为高超声速空气动力学ꎮE ̄mail:184****8082@163.com通信作者简介:张怀宝(1985 )㊀男ꎬ副教授ꎬ主要研究方向为空气动力学ꎮE ̄mail:zhanghb28@mail.sysu.edu.cn中图分类号:V211ꎻV411㊀㊀文献标志码:AAbstract:Thereisacomplextransitionphenomenonintheflowfieldofatypicalhypersonicvehicleꎬwhichhasasignifi ̄cantimpactontheperformanceofthevehicle.TheeffectsofMachnumberandwalltemperatureonthetransitionofHyTRVwerestudiedbynumericalsimulationmethods.Theself ̄developedsoftwareoftheresearchgroupwasusedtocarryoutnu ̄mericalcalculations.TherangeofMachnumberwas3~8ꎬandtherangeofwalltemperaturewas150~900K.Firstlyꎬthehypersoniccorrectionsoftheγ ̄Re~θttransitionmodelandtheSSTturbulencemodelwerecarriedout.Thepressuregradientcoefficientcorrectionandthehigh ̄speedcross ̄flowcorrectionwereintroducedintotheγ ̄Re~θttransitionmodelꎬandthecom ̄pressibilitycorrectionsoftheclosurecoefficientsβ∗andβoftheSSTturbulencemodelwerecarriedout.Thenꎬthegridin ̄dependenceverificationwascarriedoutꎬandthemodifiednumericalmethodandsoftwareplatformwereconfirmedbycom ̄paringwithexperimentalresults.FinallyꎬtheeffectsofMachnumberandwalltemperatureonthetransitionlawoftheHyTRVboundarylayerwerestudied.Theresultsshowthatthetransitionareaismainlyconcentratedonbothsidesoftheuppersurfaceandthecenterlineofthelowersurface.WiththeincreaseoftheincomingMachnumberꎬthestartingpositionoftransitionontheupperandlowersurfacesisgreatlybackwardꎬandtheturbulentzoneisgreatlyreducedꎬbutitstillex ̄ists.Atthesametimeꎬthefrictioncoefficientofthelaminarflowzoneontheuppersurfaceincreasescontinuouslyꎬandthefrictioncoefficientoftheturbulentzoneonthelowersurfacedecreases.Asthewalltemperatureincreasesꎬthestartingposi ̄tionoftransitionontheupperandlowersurfacesshiftsforwardꎬthenrapidlyshiftsbackwardꎬandfinallytheturbulentzonealmostdisappears.气体物理2024年㊀第9卷Keywords:transitionꎻHyTRVꎻfrictionꎻMachnumberꎻwalltemperature引㊀言高超声速飞行器具有突防能力强㊁打击范围广㊁响应迅速等显著优势ꎬ正逐渐成为各国空天竞争的热点[1]ꎮ高超声速飞行器边界层转捩是该类飞行器气动设计中的重要问题[2]ꎮ在边界层转捩过程中ꎬ流态由层流转变为湍流ꎬ飞行器的表面摩阻急剧增大到层流时的3~5倍ꎬ严重影响飞行器的气动性能与热防护系统ꎬ转捩还会导致飞行器壁面烧蚀㊁颤振加剧㊁飞行姿态控制难度大等一系列问题ꎬ对飞行器的飞行安全构成严重的威胁[3 ̄5]ꎬ开展高超声速飞行器边界层转捩研究具有十分重要的意义ꎮ影响边界层转捩的因素很多ꎬ例如ꎬMach数㊁Reynolds数㊁湍流强度㊁表面传导热等ꎮ在高超声速流动条件下ꎬ强激波㊁强逆压梯度㊁熵层等高超声速现象及其相互作用ꎬ会使得转捩流动的预测和研究难度进一步增大[6]ꎮ目前高超声速飞行器转捩数值模拟方法主要有直接数值模拟(DNS)㊁大涡模拟(LES)和基于Reynolds平均Navier ̄Stokes(RANS)的转捩模型方法ꎬ由于前两种计算量巨大ꎬ难以推广到工程应用ꎬ基于Reynolds平均Navier ̄Stokes的转捩模型在工程实践中应用最为广泛ꎬ其中γ ̄Re~θt转捩模型基于局部变量ꎬ与现代CFD方法良好兼容ꎬ目前已经有多项研究尝试从一般性的流动问题拓展到高超声速流动转捩模拟[6 ̄9]ꎮ目前高超声速流动转捩的研究对象主要是结构相对简单的构型ꎮMcDaniel等[10]研究了扩口直锥在高超声速流动条件下的转捩现象ꎮPapp等[11]研究了圆锥在高超声速流动条件下的转捩特性ꎮ美国和澳大利亚组织联合实施的HIFiRE计划[12]ꎬ研究了圆锥形状的HIFiRE1和椭圆锥形的HIFiRE5的转捩问题ꎮ杨云军等[13]采用数值模拟方法ꎬ分析了椭圆锥的转捩影响机制ꎬ并研究了Reynolds数对转捩特性的影响规律ꎮ另外ꎬ袁先旭等[14]于2015年成功实施了圆锥体MF ̄1航天模型飞行试验ꎮ以上对高超声速流动的转捩研究ꎬ都取得了比较理想的结果ꎬ然而所采用的模型都是圆锥㊁椭圆锥等简单几何外形ꎬ这与真实高超声速飞行器有较大差异ꎬ较难反映真实的转捩特性ꎮ为了有效促进对真实高超声速飞行器的转捩问题研究ꎬ中国空气动力研究与发展中心提出并设计了一款接近真实飞行器的升力体模型ꎬ即高超声速转捩研究飞行器(hypersonictransitionresearchvehicleꎬHyTRV)[15]ꎬ模型详细的参数见参考文献[16]ꎮHyTRV外形如图1所示ꎬ其整体外形较为复杂ꎬ不同区域发生转捩的情况也不尽相同ꎮ对HyTRV的转捩问题研究能够显著提高对真实高超声速飞行器转捩特性的认识水平ꎮLiu等[17]采用理论分析㊁数值模拟和风洞实验3种方法对HyTRV的转捩特性进行了研究ꎻ陈坚强等[15]分析了HyTRV的边界层失稳特征ꎻChen等[18]对HyTRV进行了多维线性稳定性分析ꎻQi等[19]在来流Mach数6㊁攻角0ʎ的条件下对HyTRV进行了直接数值模拟ꎻ万兵兵等[20]结合风洞实验与飞行试验ꎬ利用eN方法预测了HyTRV升力体横流区的转捩阵面形状ꎮ目前ꎬ相关研究主要集中在HyTRV的稳定性特征及转捩预测两个方面ꎬ而对若干关键参数ꎬ特别是Mach数和壁面温度对转捩的影响研究还比较少ꎮ(a)Frontview(b)Sideview㊀㊀㊀图1㊀HyTRV外形Fig.1㊀ShapeofHyTRV基于此ꎬ本文采用数值模拟方法ꎬ应用课题组自研软件开展Mach数和壁面温度对HyTRV转捩流动的影响规律研究ꎮ1㊀数值方法1.1㊀控制方程和数值方法控制方程为三维可压缩RANS方程ꎬ采用结构网格技术和有限体积方法ꎬ变量插值方法采用2阶MUSCL格式ꎬ通量计算采用低耗散的通量向量差分Roe格式ꎬ黏性项离散采用中心格式ꎬ时间推进方法采用LU ̄SGS格式ꎮ壁面采用等温㊁无滑移壁面条件ꎬ入口采用Riemann远场边界条件ꎬ出口采用零梯度外推边界条件ꎮ1.2㊀γ ̄Re~θt转捩模型γ ̄Re~θt转捩模型是Menter等[21ꎬ22]于2004年提01第2期章录兴ꎬ等:Mach数和壁面温度对HyTRV边界层转捩的影响出的一种基于拟合公式的间歇因子转捩模型ꎬ在2009年公布了完整的拟合公式及相关参数[23]ꎮ许多学者也开发了相应的程序ꎬ并进行了大量的算例验证[24 ̄28]ꎬ证明了该模型具有较好的转捩预测能力ꎬ预测精度较高ꎻ通过合适的标定ꎬγ ̄Re~θt转捩模型可以适用于多种情况下的转捩模拟ꎮ该模型构建了关于间歇因子γ的输运方程和关于转捩动量厚度Reynolds数Re~θt的输运方程ꎮ具体来说ꎬγ表示该位置是湍流流动的概率ꎬ取值范围为0<γ<1ꎮ关于γ的控制方程为Ə(ργ)Ət+Ə(ρujγ)Əxj=Pγ-Eγ+ƏƏxjμ+μtσfæèçöø÷ƏγƏxjéëêêùûúú其中ꎬPγ为生成项ꎬEγ为破坏项ꎮ关于Re~θt的输运方程为Ə(ρRe~θt)Ət+Ə(ρujRe~θt)Əxj=Pθt+ƏƏxjσθt(μ+μt)ƏRe~θtƏxjéëêêùûúú其中ꎬPθt为源项ꎬ其作用是使边界层外部的Re~θt等于Reθtꎬ定义式为Pθt=cθtρt(Reθt-Re~θt)(1.0-Fθt)Reθt采用以下经验公式Reθt=1173.51-589 428Tu+0.2196Tu2æèçöø÷F(λθ)ꎬTuɤ0.3Reθt=331.50(Tu-0.5658)-0.671F(λθ)ꎬTu>0.3ìîíïïïïF(λθ)=1+(12.986λθ+123.66λ2θ+405.689λ3θ)e-(Tu1.5)1.5ꎬ㊀λθɤ0F(λθ)=1+0.275(1-e-35.0λθ)e-(Tu0.5)ꎬλθ>0ìîíïïïï在实际计算中ꎬ通过γ ̄Re~θt转捩模型获得间歇因子ꎬ再通过间歇因子来控制SSTk ̄ω湍流模型中湍动能的生成ꎮγ ̄Re~θt转捩模型与SSTk ̄ω湍流模型耦合为Ə(ρk)Ət+Ə(ρujk)Əxj=γeffτijƏuiƏxj-min(max(γeffꎬ0.1)ꎬ1.0)ρβ∗kω+ƏƏxjμ+μtσkæèçöø÷ƏkƏxjéëêêùûúúƏ(ρω)Ət+Ə(ρujω)Əxj=γvtτijƏuiƏxj-βρω2+ƏƏxj(μ+σωμt)ƏωƏxjéëêêùûúú+2ρ(1-F1)σω21ωƏkƏxjƏωƏxj模型中具体参数定义见文献[23]ꎮ1.3㊀高超声速修正原始SST湍流模型及γ ̄Re~θt转捩模型都是基于不可压缩流动发展的ꎬ为了更好地预测高超声速流动转捩ꎬ本节引入了3种重要的高超声速修正方法ꎮ1.3.1㊀压力梯度修正压力梯度对边界层转捩的影响较大ꎬ在高Mach数情况下ꎬ边界层厚度较大ꎬ进而影响压力梯度的大小ꎬ因此在模拟高超声速流动时应该考虑Mach数对压力梯度的影响ꎮ本文采用张毅峰等[29]提出的压力梯度修正方法ꎬ具体修正形式如下λᶄθ=λθ1+γᶄ-12Maeæèçöø÷其中ꎬMae为边界层外缘Mach数ꎬγᶄ为比热比ꎮ1.3.2㊀高速横流修正在原始γ ̄Re~θt转捩模型中ꎬ没有考虑横流不稳定性对转捩的影响ꎬ对于横流模态主导的转捩ꎬ原始转捩模型计算的结果并不理想ꎮLangtry等[30]在2015年对γ ̄Re~θt转捩模型进行了低速横流修正ꎬ向星皓等[9]在Langtry低速横流修正的基础上ꎬ对高超声速椭圆锥转捩DNS数据进行了拓展ꎬ提出了高速横流转捩判据ꎬ本文直接采用向星皓提出的高速横流转捩方法ꎮLangtry将横流强度引入转捩发生动量厚度Reynolds数输运方程中Ə(ρRe~θt)Ət+Ə(ρujRe~θt)Əxj=Pθt+DSCF+ƏƏxjσθt(μ+μt)ƏRe~θtƏxjéëêêùûúú式中ꎬDSCF为横流源项ꎬLangtry低速横流修正为DSCF=cθtρtccrossflowmin(ReSCF-Re~θtꎬ0.0)Fθt2其中ꎬReSCF为低速横流判据ReSCF=θtρUlocal0.82æèçöø÷μ=-35.088lnhθtæèçöø÷+319.51+f(+ΔHcrossflow)-f(-ΔHcrossflow)其中ꎬh为壁面粗糙度高度ꎬθt为动量厚度ꎬ11气体物理2024年㊀第9卷ΔHcrossflow是横流强度抬升项ꎮ向星皓提出的高速横流转捩判据ꎬ其中高速横流源项DSCF ̄H为DSCF ̄H=cCFρmin(ReSCF ̄H-Re~θtꎬ0)FθtReSCF ̄H=CCF ̄1lnhlμ+CCF ̄2+(Hcrossflow)其中ꎬCCF ̄1=-9.618ꎬCCF ̄2=128.33ꎻlμ为粗糙度参考高度ꎬlμ=1μmꎻf(Hcrossflow)为抬升函数f(Hcrossflow)=60000.1066-ΔHcrossflow+50000(0.1066-ΔHcrossflow)2其中ꎬΔHcrossflow与Langtry低速横流修正中保持一致ꎮ1.3.3㊀SST可压缩修正高超声速流动具有强可压缩性ꎬ所以在进行高超声速计算时ꎬ应该对湍流模型进行可压缩修正ꎮSarkar[31]提出了膨胀耗散修正ꎬ对SST湍流模型中的闭合系数β∗ꎬβ进行了可压缩修正ꎬWilcox[32]在Sarkar修正的基础上考虑了可压缩生成项产生时的延迟效应ꎬ使得可压缩修正在湍流Mach数较小的近壁面关闭ꎬ在湍流Mach数较大的自由剪切层打开ꎬ本文采用Wilcox提出的可压缩性修正β∗=β∗0[1+ξ∗F(Mat)]β=β0-β∗0ξ∗F(Mat)其中ꎬβ0ꎬβ∗均为原始模型中的系数ꎬξ∗=1.5ꎮF(Mat)=[Mat-Mat0]H(Mat-Mat0)Mat0=1/4ꎬH(x)=0ꎬxɤ01ꎬx>0{其中ꎬMat=2k/a为湍流Mach数ꎬa为当地声速ꎮ2㊀网格无关性验证及数值方法确认2.1㊀网格无关性验证计算采用3套网格ꎬ考虑到HyTRV的几何对称性ꎬ生成3套半模网格ꎬ第1层网格高度为1ˑ10-6mꎬ确保y+<1ꎬ流向ˑ法向ˑ周向的网格数分别为:网格1是301ˑ201ˑ201ꎬ网格2是301ˑ301ˑ201ꎬ网格3是401ˑ381ˑ281ꎮ全模下表面如图2所示ꎬ选取y/L=0中心线和x/L=0.5处ꎬ对比3套网格的表面摩阻系数ꎬ计算结果如图3所示ꎮ采用网格1时ꎬ表面摩阻系数分布与另外两个结果存在明显差异ꎻ而采用网格2和网格3时ꎬ表面摩阻系数曲线基本重合ꎬ表明在流向㊁法向和周向均满足网格无关性ꎬ后续数值计算采用网格2ꎮ图2㊀截取位置示意图Fig.2㊀Schematicdiagramoftheinterceptionlocation(a)Surfacefrictionaty/L=0(b)Surfacefrictionatx/L=0.5图3㊀采用3套网格计算得到的摩阻对比Fig.3㊀Comparisonofthefrictiondragcalculatedusingthreesetsofgrids2.2㊀数值方法和自研软件的确认采用修正后的转捩模型对HyTRV开展计算ꎬ计算工况为Ma=6ꎬ来流温度Tɕ=97Kꎬ单位21第2期章录兴ꎬ等:Mach数和壁面温度对HyTRV边界层转捩的影响Reynolds数为Re=1.1ˑ107/mꎬ攻角α=0ʎꎬ来流湍流度FSTI=0.8%ꎬ壁面温度T=300Kꎮ为方便对比分析ꎬ计算结果与参考结果均采用上下对称形式布置ꎬ例如ꎬ图4是模型下表面计算结果与实验结果对比:对于下表面两侧转捩的起始位置ꎬ高超声速修正前的转捩位置在x=0.68m附近ꎬ高超声速修正后的计算结果与实验结果吻合良好ꎬ均在x=0.60m附近ꎬ并且湍流边界层区域形状基本一致ꎬ说明修正后的转捩模型能够较好地预测HyTRV转捩的位置ꎮ(a)Calculationofthefrictiondistribution(beforehypersoniccorrection)(b)Calculationofthefrictiondistribution(afterhypersoniccorrection)(c)Experimentalresultsoftheheatfluxdistribution[17]图4㊀下表面计算结果和实验结果对比Fig.4㊀Comparisonofthecalculatedandexperimentalresultsonthelowersurface3㊀HyTRV转捩的基本流动特性计算工况采用Ma=6ꎬ攻角α=0ʎꎬ来流湍流度FSTI=0.6%ꎬ分析HyTRV转捩的基本流动特性ꎮ从图5可以看出ꎬ模型两侧和顶端均出现高压区ꎬ高压区之间为低压区ꎬ横截面上存在周向压力梯度ꎬ流动从高压区向低压区汇集ꎬ从而在下表面中心线附近和上表面两侧腰部区域均形成流向涡结构(见图6)ꎬ沿流动方向ꎬ高压区域逐渐扩大ꎬ流向涡结构的影响范围也越大ꎮ在流向涡结构的边缘位置ꎬ壁面附近的低速流体被抬升到外壁面区域ꎬ外壁面区域的高速流体又被带入到近壁面区域ꎬ进而导致流向涡结构边缘处壁面的摩阻显著增加ꎬ最终诱发转捩ꎬ这些流动特征与文献[15]的结果一致ꎮ图7显示了上下表面摩阻的分布情况ꎬ其中上表面两侧区域在x/L=0.80附近ꎬ摩阻显著增加ꎬ出现明显的转捩现象ꎬ转捩区域分布在两侧边缘位置ꎻ而下表面两侧区域在x/L=0.75附近ꎬ也出现明显的转捩ꎬ转捩区域相对集中在中心线两侧ꎮ图5㊀不同截面位置处的压力云图Fig.5㊀Pressurecontoursatdifferentcross ̄sectionlocations图6㊀不同截面位置处的流向速度云图Fig.6㊀Streamwisevelocitycontoursatdifferentcross ̄sectionlocations31气体物理2024年㊀第9卷(a)Uppersurface㊀㊀㊀㊀㊀(b)Lowersurface图7㊀上下表面摩阻分布云图Fig.7㊀Frictioncoefficientcontoursontheupperandlowersurfaces4㊀不同Mach数对HyTRV转捩的影响保持来流湍流度FSTI=0.6%不变ꎬMach数变化范围为3~8ꎮ图8是不同Mach数条件下HyTRV上下表面的摩阻分布云图ꎬ从图中可知ꎬ随着Mach数的增加ꎬ上下表面的湍流区域均逐渐减少ꎬ其中上表面两侧转捩起始位置由x/L=0.56附近后移至x/L=0.92附近ꎬ下表面两侧转捩起始位置由x/L=0.48附近后移至x/L=0.99附近ꎬ上下表面两侧转捩起始位置均大幅后移ꎬ说明Mach数对HyTRV转捩的影响很大ꎮuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(a)Ma=3uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(b)Ma=4uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(c)Ma=541第2期章录兴ꎬ等:Mach数和壁面温度对HyTRV边界层转捩的影响uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(d)Ma=6uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(e)Ma=7uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(f)Ma=8图8㊀不同Mach数条件下摩阻系数分布云图Fig.8㊀FrictioncoefficientcontoursatdifferentMachnumbers上表面选取图7中z/L=0.12的位置ꎬ下表面选取z/L=0.10的位置进行分析ꎮ从图9中可以分析出ꎬ随着Mach数的增加ꎬ上表面转捩起始位置不断后移ꎬ当Mach数增加到7时ꎬ由于湍流区的缩小ꎬ此处位置不再发生转捩ꎬ此外ꎬMach数越高层流区摩阻系数越大ꎻ下表面转捩起始位置也不断后移ꎬ当Mach数增加到8时ꎬ此处位置不再发生转捩ꎬ此外ꎬMach数越高ꎬ湍流区的摩阻系数越小ꎬ这些结论与关于来流Mach数对转捩位置影响的普遍研究结论一致ꎮ(a)Uppersurface㊀㊀㊀㊀㊀(b)Lowersurface图9㊀不同位置摩阻系数随Mach数的变化Fig.9㊀VariationoffrictioncoefficientwithMachnumberatdifferentlocations51气体物理2024年㊀第9卷5㊀不同壁面温度对HyTRV转捩的影响保持来流湍流度FSTI=0.6%及Ma=6不变ꎬ壁面温度的变化范围为150~900Kꎮ图10是不同壁面温度条件下HyTRV上下表面的摩阻分布云图ꎬ可以看出随着壁面温度的增加ꎬ上表面两侧湍流区域先是缓慢扩大ꎬ在壁面温度为500K时湍流区域快速缩小ꎬ增加到900K时ꎬ已无明显湍流区域ꎻ下表面两侧湍流区域先是无明显变化ꎬ同样当壁面温度升高到500K时ꎬ湍流区域快速缩小ꎬ当壁面温度升高到700K时ꎬ两侧已经无明显的湍流区域ꎬ相比上表面两侧湍流区域ꎬ下表面湍流区域消失得更早ꎮ由此可以得出壁面温度对转捩的产生有较大的影响ꎬ壁面温度增加到一定程度将导致HyTRV没有明显的转捩现象ꎮuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(a)T=150Kuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(b)T=200Kuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(c)T=300Kuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(d)T=500K61第2期章录兴ꎬ等:Mach数和壁面温度对HyTRV边界层转捩的影响uppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(e)T=700Kuppersurface㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀lowersurface(f)T=900K图10㊀不同壁面温度条件下摩阻系数分布云图Fig.10㊀Frictioncoefficientcontoursatdifferentwalltemperatureconditions上表面选取z/L=0.125的位置ꎬ下表面选取z/L=0.100的位置进行分析ꎮ从图11中可以分析出ꎬ随着壁面温度的增加ꎬ上表面转捩起始位置先前移ꎬ当壁面温度增加到500K时ꎬ转捩起始位置后移ꎬ转捩区长度逐渐增加ꎬ层流区域的摩阻系数逐渐增加ꎬ当壁面温度增加到700K时ꎬ该位置已不再出现转捩ꎻ下表面转捩起始位置先小幅后移ꎬ当壁面温度增加到300K时ꎬ转捩起始位置开始后移ꎬ当壁面温度增加到700K时ꎬ由于湍流区域的减小ꎬ该位置不再发生转捩ꎮ(a)Uppersurface㊀㊀㊀㊀㊀(b)Lowersurface图11㊀不同位置摩阻系数随壁面温度的变化Fig.11㊀Variationoffrictioncoefficientwithwalltemperatureatdifferentlocations为进一步分析壁面温度的影响ꎬ本文分别在上下表面湍流区选取一点(0.9ꎬ0.029ꎬ0.14)ꎬ(0.97ꎬ-0.34ꎬ0.12)ꎬ分析边界层湍动能剖面ꎬ结果如图12所示ꎮ从图中可以看到ꎬ随着壁面温度升高ꎬ边界层厚度先略微变厚ꎬ再变薄ꎬ当壁面温度升高到700K时ꎬ边界层厚度迅速降低ꎮ这些结果与转捩位置先前移再后移的结论相符合ꎬ因为边界层厚度会影响不稳定波的时间和空间尺度ꎬ边界层厚度低时ꎬ不稳定波增长速度变慢ꎬ延迟转捩发生ꎮ需要指出的是ꎬ仅采用当前使用的方法ꎬ无法从更深层71气体物理2024年㊀第9卷次揭示转捩反转的流动机理ꎬ而须另外借助稳定性分析方法ꎬ例如ꎬ使用eN方法开展基于模态的稳定性研究ꎮ文献[33]采用该手段研究了大掠角平板钝三角翼随壁温比变化出现转捩反转的内在机理:壁温比升高促进横流模态和第1模态扰动增长ꎬ抑制第2模态发展ꎬ在第1㊁2模态联合作用影响下ꎬ出现转捩反转现象ꎮ我们将在后续开展进一步研究ꎮ(a)Uppersurface(b)Lowersurface图12㊀不同位置湍动能剖面随壁面温度的变化Fig.12㊀Variationofturbulentkineticenergywithwalltemperatureatdifferentlocations6㊀结论针对HyTRV转捩问题ꎬ在Mach数Ma=3~8ꎬ壁面温度T=150~900K的条件下ꎬ基于课题组自研软件ꎬ对γ ̄Re~θt转捩模型和SST湍流模型进行了高超声速修正ꎬ研究了Mach数和壁面温度对HyTRV转捩的影响ꎬ得出以下结论:1)经过高超声速修正后的γ ̄Re~θt转捩模型和SST湍流模型能够较为准确地预测HyTRV转捩位置ꎬ并且湍流边界层区域形状与实验结果基本一致ꎻHyTRV存在多个不同的转捩区域ꎬ上表面两侧转捩区域分布在两侧边缘位置ꎬ下表面两侧转捩区域分布在中心线两侧ꎮ2)Mach数的增加会导致上下表面转捩起始位置均大幅后移ꎬ湍流区大幅缩小ꎬ但当Mach数增加到8时ꎬ湍流区仍然存在ꎬ并没有消失ꎻ上表面层流区摩阻不断增加ꎬ下表面湍流区摩阻不断减小ꎮ3)壁面温度的增加会导致上下表面转捩起始位置先前移ꎬ再后移ꎬ这与边界层厚度变化规律一致ꎬ当壁面温度增加到700K时ꎬ下表面湍流区已经基本消失ꎬ当壁面温度增加到900K时ꎬ上表面湍流区也基本消失ꎻ上表面在层流区域的摩阻系数逐渐增大ꎬ在湍流区的摩阻系数逐渐减小ꎮ致谢㊀感谢中国空气动力研究与发展中心和空天飞行空气动力科学与技术全国重点实验室提供的HyTRV模型数据和实验数据ꎮ参考文献(References)[1]㊀OberingIIIHꎬHeinrichsRL.Missiledefenseforgreatpowerconflict:outmaneuveringtheChinathreat[J].Stra ̄tegicStudiesQuarterlyꎬ2019ꎬ3(4):37 ̄56. [2]ChengCꎬWuJHꎬZhangYLꎬetal.Aerodynamicsanddynamicstabilityofmicro ̄air ̄vehiclewithfourflappingwingsinhoveringflight[J].AdvancesinAerodynamicsꎬ2020ꎬ2(3):5.[3]BertinJJꎬCummingsRM.Fiftyyearsofhypersonics:whereweᶄvebeenꎬwhereweᶄregoing[J].ProgressinAerospaceSciencesꎬ2003ꎬ39(6/7):511 ̄536. [4]陈坚强ꎬ涂国华ꎬ张毅锋ꎬ等.高超声速边界层转捩研究现状与发展趋势[J].空气动力学学报ꎬ2017ꎬ35(3):311 ̄337.ChenJQꎬTuGHꎬZhangYFꎬetal.Hypersnonicboundarylayertransition:whatweknowꎬwhereshallwego[J].ActaAerodynamicaSinicaꎬ2017ꎬ35(3):311 ̄337(inChinese).[5]段毅ꎬ姚世勇ꎬ李思怡ꎬ等.高超声速边界层转捩的若干问题及工程应用研究进展综述[J].空气动力学学报ꎬ2020ꎬ38(2):391 ̄403.DuanYꎬYaoSYꎬLiSYꎬetal.Reviewofprogressinsomeissuesandengineeringapplicationofhypersonicboundarylayertransition[J].ActaAerodynamicaSinicaꎬ2020ꎬ38(2):391 ̄403(inChinese).[6]ZhangYFꎬZhangYRꎬChenJQꎬetal.Numericalsi ̄81第2期章录兴ꎬ等:Mach数和壁面温度对HyTRV边界层转捩的影响mulationsofhypersonicboundarylayertransitionbasedontheflowsolverchant2.0[R].AIAA2017 ̄2409ꎬ2017. [7]KrauseMꎬBehrMꎬBallmannJ.Modelingoftransitioneffectsinhypersonicintakeflowsusingacorrelation ̄basedintermittencymodel[R].AIAA2008 ̄2598ꎬ2008. [8]YiMRꎬZhaoHYꎬLeJL.Hypersonicnaturalandforcedtransitionsimulationbycorrelation ̄basedintermit ̄tency[R].AIAA2017 ̄2337ꎬ2017.[9]向星皓ꎬ张毅锋ꎬ袁先旭ꎬ等.C ̄γ ̄Reθ高超声速三维边界层转捩预测模型[J].航空学报ꎬ2021ꎬ42(9):625711.XiangXHꎬZhangYFꎬYuanXXꎬetal.C ̄γ ̄Reθmodelforhypersonicthree ̄dimensionalboundarylayertransitionprediction[J].ActaAeronauticaetAstronauticaSinicaꎬ2021ꎬ42(9):625711(inChinese).[10]McDanielRDꎬNanceRPꎬHassanHA.Transitiononsetpredictionforhigh ̄speedflow[J].JournalofSpacecraftandRocketsꎬ2000ꎬ37(3):304 ̄309.[11]PappJLꎬDashSM.Rapidengineeringapproachtomodelinghypersoniclaminar ̄to ̄turbulenttransitionalflows[J].JournalofSpacecraftandRocketsꎬ2005ꎬ42(3):467 ̄475.[12]JulianoTJꎬSchneiderSP.InstabilityandtransitionontheHIFiRE ̄5inaMach ̄6quiettunnel[R].AIAA2010 ̄5004ꎬ2010.[13]杨云军ꎬ马汉东ꎬ周伟江.高超声速流动转捩的数值研究[J].宇航学报ꎬ2006ꎬ27(1):85 ̄88.YangYJꎬMaHDꎬZhouWJ.Numericalresearchonsupersonicflowtransition[J].JournalofAstronauticsꎬ2006ꎬ27(1):85 ̄88(inChinese).[14]袁先旭ꎬ何琨ꎬ陈坚强ꎬ等.MF ̄1模型飞行试验转捩结果初步分析[J].空气动力学学报ꎬ2018ꎬ36(2):286 ̄293.YuanXXꎬHeKꎬChenJQꎬetal.PreliminarytransitionresearchanalysisofMF ̄1[J].ActaAerodynamicaSinicaꎬ2018ꎬ36(2):286 ̄293(inChinese).[15]陈坚强ꎬ涂国华ꎬ万兵兵ꎬ等.HyTRV流场特征与边界层稳定性特征分析[J].航空学报ꎬ2021ꎬ42(6):124317.ChenJQꎬTuGHꎬWanBBꎬetal.Characteristicsofflowfieldandboundary ̄layerstabilityofHyTRV[J].ActaAeronauticaetAstronauticaSinicaꎬ2021ꎬ42(6):124317(inChinese).[16]陈坚强ꎬ刘深深ꎬ刘智勇ꎬ等.用于高超声速边界层转捩研究的标模气动布局及设计方法.中国:109969374B[P].2021 ̄05 ̄18.ChenJQꎬLiuSSꎬLiuZYꎬetal.Standardmodelaero ̄dynamiclayoutanddesignmethodforhypersonicboundarylayertransitionresearch.CNꎬ109969374B[P].2021 ̄05 ̄18(inChinese).[17]LiuSSꎬYuanXXꎬLiuZYꎬetal.Designandtransitioncharacteristicsofastandardmodelforhypersonicboundarylayertransitionresearch[J].ActaMechanicaSinicaꎬ2021ꎬ37(11):1637 ̄1647.[18]ChenXꎬDongSWꎬTuGHꎬetal.Boundarylayertran ̄sitionandlinearmodalinstabilitiesofhypersonicflowoveraliftingbody[J].JournalofFluidMechanicsꎬ2022ꎬ938(408):A8.[19]QiHꎬLiXLꎬYuCPꎬetal.Directnumericalsimulationofhypersonicboundarylayertransitionoveralifting ̄bodymodelHyTRV[J].AdvancesinAerodynamicsꎬ2021ꎬ3(1):31.[20]万兵兵ꎬ陈曦ꎬ陈坚强ꎬ等.三维边界层转捩预测HyTEN软件在高超声速典型标模中的应用[J].空天技术ꎬ2023(1):150 ̄158.WanBBꎬChenXꎬChenJQꎬetal.ApplicationsofHyTENsoftwareforpredictingthree ̄dimensionalboundary ̄layertransitionintypicalhypersonicmodels[J].AerospaceTechnologyꎬ2023(1):150 ̄158(inChinese). [21]MenterFRꎬLangtryRBꎬLikkiSRꎬetal.Acorrelation ̄basedtransitionmodelusinglocalvariables PartⅠ:modelformulation[J].JournalofTurbomachineryꎬ2006ꎬ128(3):413 ̄422.[22]LangtryRBꎬMenterFRꎬLikkiSRꎬetal.Acorrelation ̄basedtransitionmodelusinglocalvariables PartⅡ:testcasesandindustrialapplications[J].JournalofTur ̄bomachineryꎬ2006ꎬ128(3):423 ̄434.[23]LangtryRBꎬMenterFR.Correlation ̄basedtransitionmodelingforunstructuredparallelizedcomputationalfluiddynamicscodes[J].AIAAJournalꎬ2009ꎬ47(12):2894 ̄2906.[24]孟德虹ꎬ张玉伦ꎬ王光学ꎬ等.γ ̄Reθ转捩模型在二维低速问题中的应用[J].航空学报ꎬ2011ꎬ32(5):792 ̄801.MengDHꎬZhangYLꎬWangGXꎬetal.Applicationofγ ̄Reθtransitionmodeltotwo ̄dimensionallowspeedflows[J].ActaAeronauticaetAstronauticaSinicaꎬ2011ꎬ32(5):792 ̄801(inChinese).[25]牟斌ꎬ江雄ꎬ肖中云ꎬ等.γ ̄Reθ转捩模型的标定与应用[J].空气动力学学报ꎬ2013ꎬ31(1):103 ̄109.MouBꎬJiangXꎬXiaoZYꎬetal.Implementationandcaliberationofγ ̄Reθtransitionmodel[J].ActaAerody ̄namicaSinicaꎬ2013ꎬ31(1):103 ̄109(inChinese). [26]郭隽ꎬ刘丽平ꎬ徐晶磊ꎬ等.γ ̄Re~θt转捩模型在跨声速涡轮叶栅中的应用[J].推进技术ꎬ2018ꎬ39(9):1994 ̄2001.91气体物理2024年㊀第9卷GuoJꎬLiuLPꎬXuJLꎬetal.Applicationofγ ̄Re~θttran ̄sitionmodelintransonicturbinecascades[J].JournalofPropulsionTechnologyꎬ2018ꎬ39(9):1994 ̄2001(inChinese).[27]郑赟ꎬ李虹杨ꎬ刘大响.γ ̄Reθ转捩模型在高超声速下的应用及分析[J].推进技术ꎬ2014ꎬ35(3):296 ̄304.ZhengYꎬLiHYꎬLiuDX.Applicationandanalysisofγ ̄Reθtransitionmodelinhypersonicflow[J].JournalofPropulsionTechnologyꎬ2014ꎬ35(3):296 ̄304(inChi ̄nese).[28]孔维萱ꎬ阎超ꎬ赵瑞.γ ̄Reθ模式应用于高速边界层转捩的研究[J].空气动力学学报ꎬ2013ꎬ31(1):120 ̄126.KongWXꎬYanCꎬZhaoR.γ ̄Reθmodelresearchforhigh ̄speedboundarylayertransition[J].ActaAerody ̄namicaSinicaꎬ2013ꎬ31(1):120 ̄126(inChinese). [29]张毅锋ꎬ何琨ꎬ张益荣ꎬ等.Menter转捩模型在高超声速流动模拟中的改进及验证[J].宇航学报ꎬ2016ꎬ37(4):397 ̄402.ZhangYFꎬHeKꎬZhangYRꎬetal.ImprovementandvalidationofMenterᶄstransitionmodelforhypersonicflowsimulation[J].JournalofAstronauticsꎬ2016ꎬ37(4):397 ̄402(inChinese).[30]LangtryRBꎬSenguptaKꎬYehDTꎬetal.Extendingtheγ ̄Reθtcorrelationbasedtransitionmodelforcrossfloweffects(Invited)[R].AIAA2015 ̄2474ꎬ2015. [31]SarkarS.Thepressure ̄dilatationcorrelationincompressi ̄bleflows[J].PhysicsofFluidsAꎬ1992ꎬ4(12):2674 ̄2682.[32]WilcoxDC.Dilatation ̄dissipationcorrectionsforadvancedturbulencemodels[J].AIAAJournalꎬ1992ꎬ30(11):2639 ̄2646.[33]马祎蕾ꎬ余平ꎬ姚世勇.壁温对钝三角翼边界层稳定性及转捩影响[J].空气动力学学报ꎬ2020ꎬ38(6):1017 ̄1026.MaYLꎬYuPꎬYaoSY.Effectofwalltemperatureonstabilityandtransitionofhypersonicboundarylayeronabluntdeltawing[J].ActaAerodynamicaSinicaꎬ2020ꎬ38(6):1017 ̄1026(inChinese).02。
第 39 卷第 6 期2023 年12 月结构工程师Structural Engineers Vol. 39 , No. 6Dec. 2023含正交各向异性夹层的Kirsch问题的弹性解AHEHEHINNOU OUGBE ANSELME1,2高梦园1,2,*(1.浙江大学土木工程系,杭州 310058; 2.浙江大学平衡建筑研究中心,杭州 310007)摘要针对含有正交各向异性圆环的Kirsch问题,利用级数展开求解材料内部的弹性控制方程,获得了两种不同荷载下各向异性圆环嵌入无限基体的位移/应力解析表达式,其中圆环材料的形态分别考虑为径向各向异性、环向各向异性、横观各向同性。
在该理论退化验证良好的基础上,对各向异性材料进行力学分析。
结果表明:不同形态的圆环对孔口周围的应力分布有较大的影响。
另外,随着各向异性圆环的厚度增大,圆环和基体之间界面处的应力值也增大。
关键词拓展Kirsch问题,正交各向异性材料,二维弹性解,应力集中Elastic Solutions of Kirsch Problem with OrthogonalAnisotropic InterlayerAHEHEHINNOU OUGBE ANSELME1,2GAO Mengyuan1,2,*(1.School of Civil Engineering,Zhejiang University, Hangzhou 310058, China;2.Center for Balance Architecture,Zhejiang University, Hangzhou 310007, China)Abstract The elastic governing equations within the material are solved through the utilization of a series expansion for the Kirsch problem,incorporating an orthogonal anisotropic hollow layer. Analytical expressions for the displacement/stress of anisotropic hollow layers embedded in an infinite matrix are obtained for two different loads.The material morphology of the hollow layer is considered to be radially anisotropic,circumferentially anisotropic,and transversely isotropic.A mechanical analysis of the anisotropic materials is performed based on the well-verified theoretical degradation. The results show that the various configurations of the hollow layers have a great effect on the stress distribution around the orifice.In addition,as the thickness of the anisotropic circular hole increases,the stress value at the interface between the circular hole and the substrate increases.Keywords expanding the kirsch problem,orthotropic materials,two-dimensional elastic solution,stress concentration0 引言钻井工程中井壁失稳问题是比较棘手的井下事故,井壁失稳通常表现为井眼漏失和井壁坍塌[1-2]。
1172009全国结构动力学学术研讨会 安徽省安庆市,2009.10.28-31 中国振动工程学会结构动力学专业委员会
弹性动力学问题一种新的时空域边界积分方程i 姚振汉 清华大学航天航空学院工程力学系, 北京, 100084 Email: demyzh@tsinghua.edu.cn 摘要:弹性动力学问题传统的时空域边界积分方程采用含时间的基本解基于动力学互等定理来建立。弹性动力学含时间的基本解是在无限弹性空间某点于某瞬时作用单位集中力脉冲的解,其中不仅含有压力波、剪切波,还有波速介于两者之间的Laplace波。本文采用加权余量格式由弹性动力学偏微分方程初边值问题出发导出一种新的时空域边界积分方程。方程中只分别利用于某瞬时会聚于弹性体边界某点的球面会聚压力波和剪切波作为核函数,从而使方程显著简化。由此建立的边界元法将比传统方法具有更高的计算效率。
关键词:弹性动力学,边界积分方程,边界元法,球面会聚压力波,球面会聚剪切波
引言 众所周知,边界元法是比有限元法稍晚几年发展起来的,最早可以看到关于间接法的一系列工作,其中求解的边界未知量并不是原问题未知场变量的边界值,而是为求解而引进的辅助变量。最早的间接法边界积分方程方法的文献可追溯到1958年(Smith 和 Pierce用于位势问题)。直接法边界积分方程方法的文献出现得稍晚一些,1963年 Jaswon将其用于位势问题。1967年Rizzo发表了关于弹性静力学问题直接法边界积分方程方法的论文,我国从事固体力学边界元法研究的一些作者曾经把它作为边界元法的第一篇文献。1968年Cruse和Rizzo就发表了弹性动力学问题直接法边界积分方程方法的文章[1, 2]。弹性动力学问题在重大工程问题中广泛存在,因此弹性动力学是固体力学边界元法中最重要的研究领域之一。在近年Aliabadi的边界元法专著[3]中也有专门的一章。 上述最早的弹性动力学边界积分方程方法的文献将边界元结合Laplace变换,然后求解变换域中的椭圆型方程,后来Manolis和Beskos对其做了一些改进[4]。时间-空间域边界元描述最早是由Cole、Kosloff和Minster于1978年对反平面问题给出的[5],后来Niwa、Kobayashi和Kitahara给出了一般形式的描述[6]。进一步的改进还可见于Antes[7],Karabalis和Beskos[8],以及Mansur等的文献[9]。 基于弹性动力学方程的问题除弹性波问题之外还有弹性体振动问题。主要对于后者,Nardini和Brebbia基于弹性静力学描述导出了质量阵和刚度阵[10],后来发展成为双重互易法,用于将惯性力的域内积分化为边界积分。弹性动力学边界元法在广泛的应用中受到重视,还进一步发展了用于土壤-结构相互作用和动态断裂力学的方法。 弹性动力学传统的时空域边界积分方程采用含时间的基本解、基于动力学互等定理来建立。该基本解是在无限弹性空间某点于某瞬时作用单位集中力脉冲的解,其中不仅含有压力波、剪切波,还有波速介于两者之间的Laplace波。
i 此项研究得到国家自然科学基金资助(10602029) 118
本文采用建立边界积分方程的一般方法,由加权余量格式从弹性动力学偏微分方程初边值问题出发,导出一种新的时空域边界积分方程。方程中只利用某瞬时会聚于弹性体边界某点的球面会聚压力波和剪切波作为核函数,使方程显著简化,从而使建立的边界元法具有更高的计算效率。 弹性动力学传统的时空域边界积分方程
弹性动力学传统的时空域边界积分方程采用弹性动力学方程含时间的基本解、基于弹性动力学的互等定理来建立。
1. 弹性动力学含时间的基本解 弹性动力学传统的时空域边界积分方程所采用的含时间基本解所满足的方程是 ()),(Q)(P,,22,2221tuucucckiskisjjkisijkjτδρρρ∆∆−=−+− (1)
其中()s
P, ; Q, kjutτ
的物理意义是τ瞬时作用在无限弹性体中一点P沿kx方向的单位集中力脉冲所引起
的弹性体任意一点Q在t瞬时jx
方向的位移分量;
122, GGccλ
ρρ+
==,分别为压力波和剪切波的波速,ρ为弹性材料的密度,
kiδ是Kroneckerδ,
(), tτ∆,()P, Q∆
为Dirac∆函数,
()()()()()(), 0 , , d, d1P, Q0 PQ, P, QdQ1Vttttt
Vτττττ
∞∞
−∞−∞∆=∀≠∆=∆=
∆=∀≠∆=∫∫
∫
该基本解的具体公式为 ()()s
212
222122122
1P, Q, 3,,HH411 +,,, , , kikikikikitrrutrrtt
rccr
rrrrrttt
cccccc
τδπρ
δ
⎧⎡⎤⎛⎞⎛⎞′
⎪
′′;=−−−−
⎢⎥⎨
⎜⎟⎜⎟
⎢⎥⎝⎠⎝⎠⎪
⎣⎦⎩
⎫⎡⎤⎛⎞⎛⎞⎛⎞
⎪
′′′∆−∆+∆
⎢⎥⎬⎜⎟⎜⎟⎜⎟
⎢⎥⎝⎠⎝⎠⎝⎠⎪⎣⎦⎭
(2)
式中ttτ′=−,H为Heaviside函数, ()()()1 H, H, d0 ttatatatatta
′
−∞
′∀>
⎧
′′=−=∆
⎨
′∀<
⎩
∫。
与上述基本解相应的面力基本解为 ()2122
1P, q, ,2,2,,4 2,,22,Skikiikkikikikirrttrnrnrr
rrnrn
cr
rrrn
rnrrrc
ψχχτδ
π
χψχχ
⎡∂∂∂⎛⎞⎛⎞⎛⎞
;=−+−−⎜⎟⎜⎟⎜⎟⎢
∂∂∂⎝⎠⎝⎠⎝⎠⎣
⎤⎛⎞∂∂∂∂⎛⎞
−+−−−⎥⎜⎟
⎜⎟
∂∂∂∂⎝⎠⎥
⎝⎠⎦
(3)
其中 223212
222211
1HH,
213, ,
crrr
tttt
ccrcr
crr
tt
rcrcc
ψ∆χψ∆∆
⎡⎤⎛⎞⎛⎞⎛⎞
′′′′=−−−+
⎢⎥⎜⎟⎜⎟⎜⎟
⎢⎥⎝⎠⎝⎠⎝⎠⎣⎦
⎛⎞⎛⎞′′=−−
⎜⎟⎜⎟
⎝⎠⎝⎠
⎥⎦⎤⎢⎣⎡⎟⎟⎠⎞⎜⎜⎝⎛′∆+⎟⎟⎠⎞⎜⎜⎝⎛′∆+⎥⎦⎤⎢⎣
⎡
⎟⎟⎠⎞⎜⎜⎝⎛′∆+⎟⎟⎠⎞⎜⎜
⎝
⎛′∆−−=
∂
∂
⎥⎦⎤⎢⎣
⎡
⎟⎟
⎠
⎞
⎜⎜⎝⎛′∆+⎟⎟
⎠
⎞
⎜⎜
⎝
⎛′∆−−=
∂
∂
1112212222222222
,,1,,13
,,1
crtcrcrtrcccrtcrcrtrrrcrtcrcrtrrr
χχχψ 119
以上公式中q表示弹性体任意的边界点,当基本解的源点,即集中力脉冲的作用点P趋于边界点时,也用小写字母p来表示。 这一组含时间的基本解中不仅包含弹性体中的压力波和剪切波,而且还含有波速介于两者之间的Laplace波。 2. 弹性动力学时空域边界积分方程 Graf [11]把弹性静力学的Betti功互等定理推广到了弹性动力学问题。考虑同一弹性体同一时间-空间域的两个独立的弹性动力学状态,()()()()()()()11101011, , , , ,
iiiiiiutfuuf
和
()()()()()()()22202022, , , , ,
iiiiiiutfuuf
,互等定理为
∫∫∫∫∫∫−+=−+
ViiViiSiiViiViiSii
QVtQutQuQVtQutQfqStqutqtQVtQutQuQVtQutQfqStqutqt)(d),(*),()(d),(*),()(d),(*),(
)(d),(*),()(d),(*),()(d),(*),(
)1()2()1()2()1()2()2()1()2()1()2()1(
ρρ
(4)
式中*号表示卷积积分。 若将待求的状态作为第1组状态,而将基本解状态作为第2组状态,则
)Q;(p,)Q,;(p,)(Q,)Q;(p,)Q,;(p,)(Q,),(Q)(p,)(Q,)q;(p,)q,;(p,)(q,)Q;(p,)Q,;(p,)(Q,)q;(p,)q,;(p,)(q,)2()2()2()2()2()2(τττττδττττττ−≡=−≡=∆∆=−≡=−≡=−≡=tututututututtftttttttutututututu
skiskiiskiskii
kiiskiskii
skiskiiskiskii
方程(4)中的卷积积分应表示为类似如下公式的形式: "",d)(q,)q;(p,)(q,*)(q,0)2()1(∫
−≡
t
tiskiiittututtτττ
于是方程(4)可以改写成
∫∫∫∫∫∫∫∫∫∫−−−+−=−−−+
VttiskiSttiskiVttiskiVttiskiSttiskiiki
VutuSttuVFtuVutuSutttuC00000(Q)dd)(Q,)Q;(p,(q)dd)(q,)q;(p,(Q)dd)(Q,)Q;(p,
(Q)dd)(Q,)Q;(p,(q)dd)(q,)q;(p,)(p,(p)
τττρττττττρτττρτττ
(5)
将方程两端的最后一个积分分别对时间进行分部积分,就可以最终得到如下时空域的位移边界积分方程:
∫∫∫∫∫∫∫∫−+−+−=−+−+
ViskiSttiskiVttiskiViskiSttiskiikiVtuttuSttuVFtuVtuttuSutttuC(Q)d)(Q,)Q;(p,(q)dd)(q,)q;(p,(Q)dd)(Q,)Q;(p,(Q)d)(Q,)Q;(p,(q)dd)(q,)q;(p,)(p,(p)0000
000
ρττττττρρτττ
(6)
这一边界积分方程可以用时空域边界元法求解。 弹性动力学一种新的时空域边界积分方程 为了得到一种新的时空域边界积分方程,本文采用加权余量格式直接从弹性动力学的偏微分方程初边值问题出发进行推导。
1. 由加权余量格式推导边界积分方程 弹性动力学方程是一个场方程,因此对该方程求散度和旋度均应等于零,即 ()[]()[]⎪⎩
⎪⎨
⎧
=−++−=−++−
0)(Q,)(Q,)(Q,)(Q,0)(Q,)(Q,)(Q,)(Q,
,,22,2221,,22,2221
mskiijjiijjkmi
iskiijjiijj
tutftuctucce
tutftuctucc
ρρρρρρ
(7)
这组方程是和原始的弹性动力学方程等价的。和它等价的加权余量格式的积分方程是