第21卷第1期2008年2月
振 动 工 程 学 报
JournalofVibrationEngineeringVol.21No.1
Feb.2008
流固耦合谐振分析的局部变分原理及其杂交元算式
苏海东1,黄玉盈2
(1.长江科学院材料与结构研究所,湖北武汉 430010;2.华中科技大学土木工程与力学学院,湖北武汉 430074)摘要:针对埋入二维无限域流场中一般结构的谐耦振问题,提出了相应的局部变分原理,并用圆形人工边界将无限域分成两部分,在人工边界内的结构及流场采用有限元数值解,人工边界外无穷大域则采用满足Sommerfeld辐射条件的解析解。通过构造泛函,使得提出的变分方程和给定谐耦振的边值问题完全等价。在此基础上进一步建立了杂交元算式,它能保证在圆形人工边界上数值解和解析解两种场函数及其导数的相容性。算例检验表明本文方法结果正确,而且具有很高的计算效率。通过基于大型软件的二次开发,将新的理论方法融入到常规的流固耦合分析中,极大地简化了编程和计算,为类似的理论研究提供了简便的数值计算实现途径。关键词:流固耦合振动;局部变分原理;杂交元
中图分类号:O422;O343 文献标识码:A 文章编号:1004-4523(2008)01-0091-05
度良好,比有限元法省时[11]。
引 言
对于无限域流场的Helmholtz或Laplace方程的外问题,数值计算不可能将区域离散到无穷远。通常的做法是,首先用人工边界截出一个有界域,并在该边界上提出合适的边界条件,然后在该区域内使用有限元法。但是怎样保证它的计算精度,其中如何引入适当的人工边界和合理的边界条件,是有限元法求解此类问题的关键。Givoli和Keller以及Harari等曾提出Dirichlet-to-Neumann(DtN)映射方法
;文献[4～6]则分别采用了不同的近似满足Sommerfeld辐射条件的无反射边界条件;Bettess
[7]
[1～3]
本文将其推广到分析埋入无限域流场中二维(或三维)结构的谐耦振问题,导出了耦振更一般形式的变分方程,再借助离散手段得到了相应的杂交
元算式。最后通过基于MARC软件的二次开发,方便实施了全部数值求解,并针对不可压缩流场中二维结构谐耦振的计算验证了方法的有效性。不难看出文献[11]研究一维结构的结果是本文的一个特例。
1 控制方程及边界条件
如图1所示,无限域流场中的二维结构K的表面作用谐激励分布力Teikt,其中k为激励频率,i=-1。设结构内表面为Γs,流固交界面为Γo。现用
等人还发展了一种无限元法(IEM),并经Burnett
等人的改进[8,9],能满足Sommerfeld条件,有效地提高了求解精度,但需要计算大量的中间数据。处理无限域问题的另一种思路是建立局部变分原理,它也是借助简单外形的人工边界将整个域分为两部分,然后通过构造合适的泛函,使提出的变分方程能够保证内域有限元解和外域解析解在交界面上的相容性。这种方法具有严格的理论基础,由于解析解的级数收敛快,通过杂交元(HybridElement)实现数值求解,计算量小。黄玉盈等研究了变水深环境下中厚度浮板的局部变分原理,收到了较好的效果。最近作者也应用这类局部变分原理解决了水中弹性环(一维结构)的声辐射问题,又一次表明精
[10]
图1 无限域流场中的结构
收稿日期:2007-02-26;修订日期:2007-09-30
基金项目:高校博士点基金项目(20040487013)资助92
振 动 工 程 学 报第21卷
f将无限域外流场划分为一适当半径R的圆周Γ
K1+K2。
1-2dk∫
K
1
(pf1,i,i+kpf1)Wpf1dK1-
2
不计结构阻尼,假设流场无粘、无旋,具有可压缩性。结构K内的应力和位移分别表示为e和u,流
场K1和K2内的流体压力分别为pf1和pf2,结构和流体的质量密度分别为m和d。系统的控制方程及边界条件为
eij,j+mkui=0,
ijnj-Ti=e0,eijnj+pf1ni=0,
2
∫Γ
o
pf12
pf1dΓo+-dkunW n pf2
n
dΓf+
∫[p
Γ
f
f2
-pf1]
W
∫Γ
f
pf1 pf2
-Wpf1dΓf+
n n
(10)
inK (1)
s (2)onΓ
onΓo (3)o (4)onΓ
∫Γ
f
1 pf2 Wpf2
Wpf2-pf2dΓf=0
2 n n
式中 考虑到流体Γo边界的法线n和固体Γo边界的法线n方向相反,即n=-n。
由方程(10)可知,只要最后一项积分为0,则由 pf2
的任意性,所有的控制方程和边 n
界条件将得到满足。方程(10)之所以被称为局部变Wui,Wpf1和W分原理,是因为其中不涉及人工圆边界以外K2的积分运算,这是它的重要优点。
现在剩下的问题是要进一步证明
1 pf2 Wpf2
f2-pf2dΓf=0WpΓ2 n nf
pf12
=dkun, n
pf1,i,i+kpf1=0,
22
inK1 (5)
2 (6)pf2,i,i+kpf2=0,inK
pf2
limrf2=0, r→∞ r+ikp222 r=x+y,onΓ∞ (7)
pf1 pf2
pf1=pf2, =,onΓf (8)
n n
其中,式(1),(5)和(6)分别为结构和流体的平衡微
∫(11)
在极坐标系下,满足条件(7)的Helmholtz方程(6)的解为
∞
分方程(按指标形式给出,对于平面问题,i=1,2,
kj=1,2),k=,c为流场中波速;式(2)中Ti为结构
c
所受分布力T的坐标分量;式(3)和(4)为流固耦合边界条件,式(8)反映了圆形人工边界上的数值解pf1和解析解pf2两种场函数及其导数的相容性,n代表边界面的法向;式(7)为Sommerfeld辐射条件,r以圆心为原点。
pf2(r,U)=
∑
H(n2)(kr)[ancos(nU)+bnsin(nU)]
(12)
n=0
式中 Hn(2)(kr)为第二类n阶Hankel函数;an,bn为级数系数。
对于K2,根据Green恒等式有
pf2 Wpf2
dΓ= Wpf2-pf2
Γ n n2 无限域流场中一般结构谐耦振的局
部变分原理
根据系统的总势能构造泛函如下
112
J=KeijXijdK-kmuiuidK-K22
1TiuidΓs+pf1undΓo+2 ΓΓdkso
∫∫[Wp p
K
f2
2
2
f2
-pf2 Wpf2]dK2
2
(13)
∫∫∫∫1ppdK-k∫pdK2∫1 pp-pdΓ∫2 n
K
f1,i
f1,i
2
1
222
式中 Γ=Γf+Γ∞。另由 pf2=-kpf2及 Wpf2=-k2Wpf2,可得
∫
f
(Γ+Γ)
∞
pf2 Wpf2
Wpf2-pf2dΓ=0(14) n n
r- pf2
=-ik r
rpf2,因此
将式(7)改写成 pf2
=W r(9)
Wpf2
= r
K
2
f1
1
+
1
ikWpf2,从而可得
Γ
f2f1
f2
f
f
∫
Γ
∞
pf2 Wpf2
f2-pf2dΓ=0,代入式(14)后Wp n n 下面将证明,上述泛函变分后,由WJ=0将得到
前一节所有的控制方程和边界条件。对式(9)中的逐项进行变分,得WJ=
立即看出式(11)成立。
∫(en
Γ
ij
s
j
s+-Ti)WuidΓ
∫(en+
Γ
ij
j
o
3 对应的杂交元矩阵算式
pf1ni)
将式(12)的级数截断后,它可表示为
pf2=Cq(15a)
o-WuidΓ
∫(e
Kij,j
+mkui)WuidK+
2
第1期
(2)
苏海东,等:流固耦合谐振分析的局部变分原理及其杂交元算式
(2)
(2)
93
式中 C=[H0(kr)H1(kr)cos(U)H1(kr)·
(2)(2)
sin(U)…Hn(kr)cos(nU) Hn(kr)sin(nU)](15b)
T
q=a0a1b1…anbn(15c)
将K和K1剖分有限单元,则每个单元内{u}e=Ns{d}e,pf1=Nf{p}e,其中,Ns和Nf分别为结构单元和流体单元的形函数,{d}e和{p}e分别为结构单元节点位移和流体单元节点压力。
将式(9)变分后用矩阵表示(均是单元矩阵的组集)
Wd[(K-kM)d-f-Gp]+
1{WpT[(H-k2E)p-2dk2TTdkGd-Lq]+Wq[Aq-Lp]}=0
式中 K=
T
2
T
e
这样,式(19)就可简化为
K-kM-G
2
-G
T
2
(H+H1)/dkd
=pf
(22)0求解方程(19)或(22)后,要分析的流固谐耦振问题就迎刃而解。
4 基于MARC软件二次开发实施数值求解
由方程(19)可知,要得到系统的响应,需形成结构的单元刚度矩阵和质量矩阵、流体单元矩阵和流固耦合矩阵,进行组集后再进行方程求解,这意味着要编制一套完整的流固耦振分析程序,工作量很大。实际上,现有的大型有限元分析软件,如MARC,ANSYS等均具有考虑流体压缩性的流固耦合分析功能。应用这些软件并加以二次开发,可极大地简化以上工作量。
考察方程(19)发现,本文推导的方程与常规有限元流固耦合方程不同之处仅仅在于含流体矩阵部分的附加项H1(式(20)),因此应用MARC软件的接口技术,不难将该附加项组集到流体矩阵中,从而就可通过常规的流固耦合分析步骤获得系统的响应,极大地简化了编程工作。具体做法如下:
(1)在MARC软件中按通常的步骤将K和K1分别划分为固体单元和流体单元,选中与Γf相关的流体单元的边,定义边界条件,输出计算文件。(2)自编程序,读入计算文件中的Γf边界条件信息,通过式(17)计算A,L矩阵,并形成附加项H1。
(3)式(20)中的H1为组集后的流体矩阵附加项,而MARC软件只提供对单元矩阵的操作功能——用户单元userelement。本文将Γf上的流体单元边界节点两两相连,定义为一维用户单元,并赋予其H1矩阵相应位置上的系数。对于已经用过的系数,令其为0,使得所有一维单元的系数组集后正好等于H1。
(4)返回到MARC计算程序中,通过用户子程序uselem读入上述H1,按通常的步骤计算系统响应,最后由图形界面显示计算结果。
除上述方法之外,还有另外一种基于MARC软件的二次开发技术:不使用用户单元(Userelement)接口,而是在自编程序中读入由MARC软件生成的每个单元矩阵及耦合矩阵,并计算H1矩阵,然后按照节点顺序对号入座,组集成式(19)或式(22)的总体矩阵,最后调用Fortran函数库求解方程组。这种方法的通用性甚至比前一种方法更强。(16)
∫
s
s
BDBdK(B和D分别为应变矩阵和K
T
弹性矩阵), M=
f=
T
Γ
∫
sfi
T
mNsNsdK,K
fT
T
s
o
∫NTdΓ, G=∫NnNdΓ,
N NH=∫dK,E=∫NNdK,
x x1C+ CCdA=∫C Γ,2 n n CL=∫NdΓ(17)
n
Γ
o
f
T
K
1
i
1
K
f
T
f1
1
T
T
Γ
f
f
T
Γ
ff
f
由Wd,Wp,Wq的任意性,经整理后得2T
K-kM-G0df-dkG
2
TTT
H-kE
2
-L
T
0-LA
这就是基于上述局部变分原理导出的杂交元算式,待求的列向量是由结构的节点位移、流体节点压力以及级数(12)的展开系数混合组成,杂交元正是由此而得名。由于圆形人工边界外域不必离散,故待求未知量大为减少。
根据方程(18)的最后一个等式消去q还可使方程缩聚,由此得
2
K-kM
-G式中
H1=-LTA-1L
流场域Laplace方程pf2,i,i=0的解为
∞
p=
q00(18)
-GT
2
(H-k2E+H1)/dkd
=pf0(19)(20)
若不考虑流体的可压缩性(k=0),满足无限远
pf2(r,U)=
∑
n=1
r
-n
[ancos(nU)+bnsin(nU)]
(21)94
振 动 工 程 学 报第21卷
5 不可压流场谐耦振的算例
由式(12)可知,pf2表达式中的H(kr)涉及较复杂的复数运算。为简便起见,本文采用不可压缩流体的算式(21)及(22)来验证局部变分原理。这时由于pf2函数形式简单,式(20)中的A,L矩阵均可以解析表达。
现在计算一个椭圆厚环(二维结构)的算例,它不同于一维薄环。图2所示为该椭圆环及其周围的流体网格。椭圆环厚度为0.2m,长半轴外半径为1m,短半轴外半径为0.8m,按需要剖分了两层网格。材料参数为:弹性模量105kN/m2,泊松比0.1,质量密度1000kg/m。在椭圆厚环水平方向的左右内侧1,2两点施加全约束,在3点处作用水平集中力激励,幅值为f0=100kN,频率为10rad/s。流体密度为500kg/m3。
3
(2)n
方案计算:
(1)令R=1.5m处节点的流体压力p=0,相当于固定边界条件;
(2)在R=1.5m处不做任何处理,相当于节点
p的流体压力法向导数=0,即反射边界条件;
n
(3)根据局部变分原理,在R=1.5m处的Γf
上引入解析解。
另外,为了对比,还增加了方案(4),采用图3所示9m厚范围的流体网格,经试算认为该网格已经离散到足够远了。
采用MARC软件计算上述几种方案(其中方案(3)应用了基于局部变分原理的二次开发技术),其结果如表1所示,表中给出了集中力作用点3两个方向(x,y)的位移幅值u,v。由表可见,方案(3)的计算结果最接近剖分很多流场网格的方案(4),这说明采用局部变分原理引入的解析解能完全反映远场效应。
表1 椭圆厚环的计算结果对比方案(1)固定条件
u/m
v/m
0.07650.0174
方案(2)反射条件0.09020.0212
方案(3)加入解析解0.08260.0192
方案(4)远场0.08170.0188
下面进一步讨论级数(21)的收敛性。表2给出了采用不同级数项数的计算结果,其中项数为0的情况就是方案(2)。结果表明收敛很快,级数取到5项(最多10项)所得精度就已足够。而从图3可看出,
图2 椭圆厚环网格及其周围0.5～0.7m范围的流体网格
用作对比的计算方案(4)在Γf以外的自由度相当多,用有限元计算十分费时,这说明局部变分原理的计算效率很高。
表2 级数项数与收敛性
级数项数
0
123451020
u/m0.090150.087040.083850.082700.0820.082620.082600.08260
v/m0.021150.020320.019490.019180.019170.019160.019160.01916
图3 椭圆厚环及其周围9m范围的流体网格
6 结束语
本文扩展了文[11]的思路,针对二维无限域流场中一般结构的谐耦振问题,提出了一个局部变分原理,用圆形人工边界将无限域分成两部分,在人工圆形人工边界的半径R=1.5m,这样,图2中的椭圆厚环附近的流体范围仅为0.5～0.7m,剖分了4层网格。在R=1.5m处的边界条件采用以下3种 第1期苏海东,等:流固耦合谐振分析的局部变分原理及其杂交元算式
95
边界内的二维结构及流场采用有限元解,人工边界外流场采用满足Sommerfeld辐射条件的解析解。在此基础上建立了杂交元算式,保证了在圆形人工边界上两种场函数及其导数的相容性。算例表明本方法精度好,计算效率高。其次,为实现本文方法,还提出了基于大型软件的二次开发技术,极大地简化了编程,为类似的理论研究提供了简便的数值计算实现途径。
值得指出,这里的算例验证虽然没有考虑流体的可压缩性,但要考虑这个因素的影响也没有实质性的困难,只是涉及复数运算,处理上稍微复杂一些。另外,本方法可进一步推广到一般的三维耦合场问题,此时只需用封闭球面替代圆形人工边界,用三维Helmholtz或Laplace方程的解析解代替文中二维解析解(12),其它步骤基本不变。因此,本文方法为将来水下结构的声-弹耦合问题的研究打下了基础。参考文献:
[1] GivoliD,KellerJB.Afiniteelementmethodfor
largedomains[J].Comput.Methods.Appl.Mech.
:41—66.Eng.,19,76
[2] KellerJB,GivoliD.Exactnon-reflectingBoundary
conditions[J].J.Comput.Phys.,19,82:172—192.
[3] HarariI,HughesTJR.Studiesofdomain-basedfor-
mulationsforcomputingexteriorproblemsofacous-tics[J].Int.J.Numer.MethodsEng.,1994,37:2935—2950.
[4] AssaadJ,DecarpignyJN,BruneelC.Applicationof
thefiniteelementmethodtotwo-dimensionalradiationproblem[J].J.Acoust.Soc.Am.,1993,94(1):562—573.
[5] KallivokasLF,BielakJ.
Time-domainanalysisof
transientstructuralacousticsproblemsbasedonthefiniteelementmethodandanovelabsorbingboundaryelement[J].J.
3480—3492.
Acoust.Soc.Am.,1993,94(6):
[6] GivoliD,KellerJB.Specialfiniteelementsforuse
withhigh-orderboundaryconditions[J].
Comput.
:199—213.Methods.Appl.Mech.Eng.,1994,119
[7] BettessP.Infiniteelements[J].Int.J.Numer.Meth-odsEng.,1977,11:53—.
[8] BurnettDS.Athree-dimensionalacousticinfiniteel-ementbasedonaprolatespheroidalmultipoleexpan-sion[J].J.Acoust.Soc.Am.,1994,96(5):2798—2816.[9] BurnettDS,
Holford
RL.
Prolateandoblate
Comput.
spheroidalacousticinfiniteelements[J].
Methods.Appl.Mech.Eng.,1998,158:117—141.[10]黄玉盈.变水深环境下中厚度浮板耦振问题的一个变
分解[J].力学学报,1988,20(5):401—410.
[11]HuangYuying,SuHaidong,XiangYu.Avariational
solutionof2-Dsound-structureinteractionproblems[J].ActaMech.Sinica,2006,22(2):127—131.
Alocalizedvariationalprincipleandhybridelementformulas
-structurecouplingvibrationsforfluid
12SUHai-dong,HUANGYu-ying
(1.ChangjiangRiverScientificResearchInstitute,Wuhan430010,China;2.SchoolofCivilEng.&Mechanics,HUST,Wuhan430074,China)
Abstract:Alocalizedvariationalprincipleispresentedforanalyzingfluid-structureinteractionproblemsunderharmonicexci-tationintwo-dimensionalinfiniteliquidfield.Theinfiniteliquidfieldisdividedintotwopartsbyanauxiliarycircle,insidewhichthestructureanditsneighboringliquidarecomputedbyFEM,andoutsidewhichananalyticalsolutionsatisfyingSom-merfeldconditionisused.Bymeansofthisvariationalprinciple,allgoverningequationsandboundaryconditionsaresatisfiedautomatically.Hybridelementformulasarealsoderivedbasedonthisprinciple.Viatheseformulas,thenumericalsolutionsareincorporatedwiththeanalyticalsolutionontheauxiliarycircle.Agivenexampleforincompressibleliquidcasehasdemon-stratedthevalidityandhighefficiencyoftheproposedapproach.Inaddition,interfaceofMARC(ageneral-purposeFEM
software)isalsoadoptedtomergethenewapproachintothetraditionalcouplinganalysis,whichimprovestheefficiencyofprogrammingandcomputationgreatly.Thepresentedapproachcanalsobeappliedinfuturestudiesonstructure-acousticproblemsofsubmergedstructures.
:fluid-structurecouplingvibration;localizedvariationalprinciple;hybridelementKeywords
作者简介:苏海东(1968—),男,工学博士,高级工程师。电话:(027)82829754;E-mail:suhd@mail.crsri.cn
因篇幅问题不能全部显示，请点此查看更多更全内容