突然開孔時孔口氣流動力特性參數(shù)的數(shù)值模擬

1、第23卷第10期 Vol.23 No.102006年 10 月 Oct. 2006文章編號：1000-4750(2006)10-0055-06 工 程 力 學 ENGINEERING MECHANICS 55突然開孔時孔口氣流動力特性參數(shù)的數(shù)值模擬盧 旦，*樓文娟(浙江大學建筑工程學院, 杭州 310027)摘 要：建筑物突然開孔時瞬態(tài)內(nèi)壓響應(yīng)的規(guī)律可以用一個二階非線性常微分方程來描述，由于開孔形狀和孔邊物理特征的復(fù)雜性，采用風洞試驗獲取方程中各項系數(shù)的方法往往比較困難。首先利用計算流體動力學的方法，模擬建筑物突然開孔瞬時的流場變化，獲得風致內(nèi)壓的Helmholtz頻率和孔口氣流的等效線性阻尼

2、。再由這兩個參數(shù)計算得到的內(nèi)壓增益并與風洞試驗結(jié)果相比較，二者結(jié)果吻合較好，表明數(shù)值計算能夠準確地模擬突然開孔結(jié)構(gòu)孔口處的氣體流動狀態(tài)。其次，利用孔口氣流振蕩曲線，采用參數(shù)擬合的方法進一步獲得了內(nèi)壓傳播方程中的各特征參數(shù)。關(guān)鍵詞：風洞試驗；風致內(nèi)壓；數(shù)值模擬；突然開孔；動力特征參數(shù)中圖分類號：TU311.3 文獻標識碼：ANUMERICAL SIMULATION OF FLOW DYNAMICAL PARAMETERS ATSUDDEN OPENINGLU Dan , *LOU Wen-juan(Department of Civil Engineering, Zhejiang Univers

3、ity, Hangzhou 310027, China)Abstract: The transient response of internal pressure following a sudden openning was described by a second-order, non-liner, ordinary differential equation. Because of the complexity of the geometry and the opening configuration, it is quite difficult for wind tunnel tes

4、t to correctly estimate the constant coefficients in the equation. The transient response of flow pattern at sudden opening is simulated by Computational Fluid Dynamics method. Both the Helmholtz frequency and equivalent viscous losses are obtained. The gain function obtained by using these paramete

5、rs coincides with the result of wind tunnel test fairly accurately. It is shown that the numerical computation can simulate the flow pattern at sudden opening correctly. Finally the flow dynamical parameters are derived by fitting the differential equation to the numerical response.Key words: wind t

6、unnel test; wind-induced internal pressure; numerical simulation; sudden openning;dynamical parameter傳統(tǒng)上人們對結(jié)構(gòu)抗風的設(shè)計往往只考慮建筑物的外表面所承受的風荷載。但在風災(zāi)天氣下，建筑物的門、窗、幕墻等圍護結(jié)構(gòu)特別容易被強風吹開，風從開孔突然涌入。無論是內(nèi)壓靜力效應(yīng)還是由于突然開孔所引起的內(nèi)壓超載和紊流引起的收稿日期：2004-12-03；修改日期：2005-04-06基金項目：國家自然科學基金(50378085)作者簡介：盧 旦(1978)，男，浙江上虞人，博士生，從事結(jié)構(gòu)風工程研究(E-

7、mail: bblu007)；56 工 程 力 學進出振蕩的頻率稱為Helmholtz頻率。研究內(nèi)壓響應(yīng)的瞬態(tài)過程，并由此獲得Helmholtz頻率及孔口氣流阻尼能夠方便地計算出內(nèi)壓響應(yīng)的增益，為進一步研究開孔結(jié)構(gòu)的風振響應(yīng)提供條件。隨著計算機硬件和軟件技術(shù)的日新月異，采用計算流體動力學(CFD)對鈍體繞流進行數(shù)值模擬分析，獲取建筑物表面風壓及風速流線的技術(shù)已經(jīng)比較成熟4,5。CFD技術(shù)不但能節(jié)約試驗成本、縮短試驗周期，而且與試驗相比可以獲得非常詳盡的資料，有助于對問題發(fā)生的機理進行研究。Sharma和Richards6采用數(shù)值模擬研究了結(jié)構(gòu)突然開孔后內(nèi)壓發(fā)生的靜脈收縮式振動現(xiàn)象，比較了紊流引

8、起的粘滯阻尼損失和內(nèi)外氣流相互作用引起的能量損耗。本文通過CFD數(shù)值模擬，獲得突然開孔結(jié)構(gòu)內(nèi)壓瞬態(tài)響應(yīng)方程中的Helmholtz頻率及孔口氣流的等效線性阻尼，計算出內(nèi)壓響應(yīng)的增益，使之與試驗結(jié)果相比較，研究利用流體動力學數(shù)值計算方法來模擬開孔結(jié)構(gòu)孔口處氣流的運動狀態(tài)，并由此獲得相關(guān)參數(shù)的可行性。方程(1)滿足伯努利方程，穩(wěn)態(tài)流中孔口處的氣流將會發(fā)生收縮，因此方程中出現(xiàn)孔口收縮系數(shù)c。Vickery和Bloxham研究認為在非穩(wěn)態(tài)流中將不出現(xiàn)孔口氣體收縮，因此方程中不應(yīng)該出現(xiàn)c，而代以另外兩個系數(shù)CI和CL：V0&&aLeaqV02&|C&+C=C (3) Cpi

9、+CL|CpipipipeA0pa2(A0pa)Helmholtz頻率的表達式也不包含流動系數(shù)cfHH(Vickery)=12Ap (4) V0aLe=L0+CIA0是氣流的有效長度，CI是慣性其中Le系數(shù)，CL是損失系數(shù)。文獻6建議，對于墻面較薄的開孔采用下式：LV&&qV2&|C&+ Cpi+CL|Cpipi2cA0pa2(A0pa)µeffrPLeV0&+C=C (5) CpipipecA0p0對于墻面較厚、孔道長度大于孔道直徑的開孔，建議采用如下方程：V0&&LqV02&|C&+ |CCpi+CLpipi

10、A0pa2(A0pa)µeffV0rPLe&+C=C (6) CpipipeA0p01 內(nèi)壓傳播方程理論迄今為止，各國學者已經(jīng)在內(nèi)壓響應(yīng)方程類型上達成共識，即認為可以用二階常微分方程來描述突然開孔結(jié)構(gòu)的內(nèi)壓響應(yīng)。Liu和Saathoff7采用非定常等熵伯努利方程對孔口處入射氣流進行分析后得到：aLeV0&&aqV02&|C&+C=C (1) Cpi+|CpipipipecA0pa2(cA0pa)式中，、a、pa分別是開孔周圍空氣的比熱、密度和壓強；A0開孔面積，V0是建筑物內(nèi)部容積；Cpi,Cpe分別為內(nèi)壓力系數(shù)和外壓力系數(shù)，q是參考點風壓，即

11、q=1/2avh，vh為參考點風速；c是氣µeff是考慮紊流的有效摩擦系數(shù)，r為流速從平均速度到達零速度時經(jīng)過的距離，P為孔口周長。盡管方程(5)、方程(6)已經(jīng)明確表達了瞬態(tài)內(nèi)壓響應(yīng)的規(guī)律，但方程各項系數(shù)的取值卻受到孔口物理性質(zhì)的影響。不同的問題對應(yīng)有不同的系數(shù)，而問題的關(guān)鍵在于開孔處流動狀態(tài)的區(qū)別。由于受到試驗條件的限制，通過風洞模型試驗獲取這些參數(shù)的方法往往會很困難。而數(shù)值模擬獲取參數(shù)的可行性和可靠性正是本文所要研究的主題。流的孔口收縮流動系數(shù)；Le=L0+A0/是孔口處空氣柱的有效長度，L0是開孔的實際深度，一般指墻體的厚度，系數(shù)依賴于開孔的形狀等其他因素需由試驗測得。上式中

12、的一至四項可分別被視作慣性項、阻尼項、彈性項和系統(tǒng)所受的外力。在無阻尼自由振動條件下(即令方程的第二項和第四項為零)求解(1)式便得到Helmholtz頻率：1cA0pafHH(Liu/S)= (2)2aLeV02 風洞試驗風洞試驗剛性模型選用優(yōu)質(zhì)木材做成，屋面為800mm×500mm的長方形，底裙高200mm，屋面板和底裙之間為固支，用膠帶紙密封。測壓管為外徑1mm的銅管，埋在屋面板上下表面中間預(yù)先設(shè)計的夾層內(nèi)。為了使從測壓點處引出的銅管不影響模型內(nèi)部風場，在屋面板下的模型內(nèi)設(shè)置空心柱工 程 力 學 57子，屋面的測壓管通過空心柱與底部的測壓模塊相連接。模型在迎風面墻板正中開孔，開

13、孔面積為A0=0.1m×0.1m。試驗時先將孔口用一擋板密封，在風洞中吹幾秒鐘后用預(yù)先連在擋板上的繩子突然拉開，模擬實際建筑門窗產(chǎn)生風致破壞的情形。模型如圖1所示。由于開孔前內(nèi)、外壓存在壓力差，因此開孔瞬間外壓的激勵類似于一個矩形脈沖。而從圖2所示的內(nèi)壓系數(shù)時程中可以看出，開孔瞬間內(nèi)壓突然增大，但由于試驗中使用的是湍流流場，實際內(nèi)壓響應(yīng)含有脈動風壓，從其時間歷程中不易看出Helmholtz振動；又由于本模型Helmholtz振動的能量較小，因此依據(jù)其功率譜也難以獲得Helmholtz 振動的信息。3 CFD數(shù)值模擬采用計算機進行數(shù)值模擬不但能對結(jié)構(gòu)進行準確地三維計算，采用非穩(wěn)態(tài)計算還

14、能實時記錄下流體的壓強、流速矢量等參數(shù)的變化，逼真地再現(xiàn)氣流的運動變化規(guī)律。而這些都是風洞試驗所難以圖1 突然開孔風洞試驗剛性模型Fig.1 Stiff model of the structure with a sudden opening或無法做到的。數(shù)值模擬對象為前文中圖1所示的風洞模型。計算流域為長×寬×高6m×3m×2m，符合風洞試驗阻塞比小于4%的要求，采用非結(jié)構(gòu)化網(wǎng)格劃分，孔口和建筑物附近網(wǎng)格加密。為提高計算效率，考慮到結(jié)構(gòu)的對稱性采用對稱建模，在計算過程中采用自適應(yīng)網(wǎng)格技術(shù)。計算網(wǎng)格模型如圖3所示。由于來流在屋檐處分離造成迎風屋面邊緣附近

15、形成高負壓區(qū)，這種高負壓區(qū)僅限于邊緣區(qū)域8，所以在屋面上表面的邊緣測壓點布置較密。考慮到室內(nèi)氣壓較均勻，所以在屋面下表面及內(nèi)墻面的測點進行均勻布置。由于模型具有對稱性，因此對模型屋蓋上、下表面的測點都采取1/2單邊布置(見圖1)，但在迎風墻面孔口附近重點研究位置則采取整個墻面對稱布置。本次風洞試驗在南京航空航天大學603研究所NH-2低速風洞中進行。試驗段截面為矩形帶小切角，高2.5m，寬3.0m，長6.0m。在試驗段的進風口設(shè)置格柵以形成一定的湍流度。數(shù)據(jù)的采集由多通道電子掃描脈動風壓測量系統(tǒng)完成。2圖3 開孔結(jié)構(gòu)模型的計算網(wǎng)格 Fig.3 Computing grid of the ope

16、ning model內(nèi)壓系數(shù)Cpi1.5因為是迎風面開孔，同時又考慮到孔道內(nèi)氣體流動具有強烈的各向異性，因此湍流模型選用可實現(xiàn)的k模型(Realizable k)。流場為均勻湍流場，參考風速14m/s，湍流度I=20%。由于脈動風對平均風的影響是通過湍流度計入運算過程，因此計算結(jié)果得到將是開孔瞬時平均風速的變化過程。計算分為兩步：首先，關(guān)閉(deactivate)結(jié)構(gòu)內(nèi)部網(wǎng)時間t/s10.50格進行穩(wěn)態(tài)計算(steady-state solution)，用來模擬開孔前的流場環(huán)境，計算收斂的結(jié)果保存作為下一步非穩(wěn)態(tài)(unsteady-state)計算的初始條件。第二步，激活(activate)結(jié)

17、構(gòu)內(nèi)部網(wǎng)格使之參與非穩(wěn)態(tài)計算。根圖2 內(nèi)壓系數(shù)時程Fig.2 Time history of internal pressure coefficient58 工 程 力 學據(jù)初步估算，結(jié)構(gòu)的Helmholtz頻率約為60Hz，由此設(shè)定計算時間步長為0.001s。值得指出的是，通過改變時間步長可以控制突然開孔的延續(xù)時間，這一點是風洞試驗無法精確做到的，同時計算時間步長的大小只會影響響應(yīng)峰值的大小(步長越小峰值越高)，但不會改變Helmholtz頻率和阻尼值。本次計算時間步總數(shù)為100步，為節(jié)省計算機資源采取自適應(yīng)迭代，每個時間步最大迭代次數(shù)為30次。與以往研究風荷載時考慮空氣為不可壓縮氣體不同的

18、是，此次問題必須把空氣視為可壓縮氣體來研究。即空氣的密度和壓強必須滿足以下公式：ppa=const (7) a內(nèi)孔口附近氣流仍是往內(nèi)回的。由此可見，出風時不但孔口有效面積較小，而且由于氣流的復(fù)雜流動，阻尼損耗也較大，整個氣流循環(huán)過程是一個嚴重的非線性過程。由于Helmholtz頻率和孔口氣流阻尼是開孔結(jié)構(gòu)內(nèi)壓響應(yīng)所固有的特性，它們只和孔口本身的物理性質(zhì)有關(guān)，因此這些參數(shù)同樣可以從進出孔口氣流的瞬態(tài)響應(yīng)中獲取。圖5是數(shù)值模擬得到的突然開孔后孔口氣流的振蕩曲線。風速v/(m/s)為提高非穩(wěn)態(tài)計算的精度，對可壓縮氣體方程的求解采用2階couple implicit格式算法。時間t/s室外 孔口室內(nèi)圖

19、5 數(shù)值模擬進出孔口氣流的振蕩曲線 Fig.5 Computational oscillating curve of flow followingsudden opening(a) 進風(0.01s)結(jié)構(gòu)開孔前孔口風速的平均值為零，在開孔之后的0.05s時間，由于局部風場環(huán)境發(fā)生了變化，風速也發(fā)生了變動。圖5可以看到，風速中明顯包含了Helmholtz振動。且氣流在開孔后0.01s時間內(nèi)是處于進風狀態(tài)，在之后的時間內(nèi)全都處于出風狀態(tài)，直至趨向于零風速。為計算Helmholtz頻率和孔口阻尼，還需對圖中原始曲線進行處理，首先采用HHT(HilbertHuang transform)方法9獲得原始

20、振蕩曲線的趨勢線(圖5中虛線表示)。該方法從本質(zhì)上講是對一個信號進行平穩(wěn)化處理，其結(jié)果是將信號中不同尺度的波動或趨勢逐級分解開來，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，每一個序列稱為一個本征模函數(shù)(Intrinsic Mode Function，IMF)分量。對經(jīng)EMD分解得到的各IMF分量再進行線使得變換后得到的結(jié)性、平穩(wěn)處理即Hilbert變換，果能夠反應(yīng)真實的物理過程。將原始數(shù)據(jù)減去所有室外 室內(nèi)(b) 出風(0.015s) 圖4 孔口附近風速矢量圖 Fig.4 Velocity vectors at the opening圖4(a)為開孔后0.01s時刻開孔位置處氣流從模型外部流入模型

21、內(nèi)部的矢量圖，圖4(b)為開孔后0.015s時刻開孔位置處氣流從模型內(nèi)部流出模型內(nèi)部的矢量圖。從圖4可見，進風和出風均有孔口收縮現(xiàn)象，并且進風時在室內(nèi)的孔口附近、出風時在室外的孔口附近均伴有渦流產(chǎn)生。此外，出風時室IMF分量的總和即為原始型號的趨勢曲線。原始曲線減去由HHT方法得到的趨勢曲線后，可以得到趨勢為零的自由衰減曲線(經(jīng)EMD分解和Hilbert變換后得到的所有IMF分量總和)。至此，工 程 力 學 59便可由常用的計算單自由度振動系統(tǒng)頻率和阻尼的方法計算Helmholtz頻率和孔口阻尼。其中，孔口氣流振蕩的等效線性阻尼及其他參數(shù)阻尼的確定可以分別通過等效線性方程和非線性常微分方程參數(shù)

22、擬合得到。(10)得出內(nèi)壓增益。圖6為分別為由CFD計算和風洞試驗得到的內(nèi)壓增益比較。試驗和數(shù)值模擬獲得的Helmholtz和等效線性阻尼系數(shù)基本一致，分別為fHH=43Hz，ceq=35Hz。從圖6中可以看出，CFD計算和風洞試驗到的內(nèi)壓增益吻合得非常好。說明利用CFD數(shù)值模擬能夠十分準確地獲得開孔結(jié)構(gòu)的Helmholtz振動和孔口氣流的阻尼。為了進一步準確獲得孔口氣流的非線性動力特性，計算孔口氣流運動的動力特征參數(shù)，也即方程(5)或方程(6)中的常系數(shù)值?？梢岳脠D5中的衰減曲線，采用多變量尋優(yōu)擬合的方法確定這些常系數(shù)。根據(jù)前文介紹，孔口氣流的動力特性只和孔口本身的物理性質(zhì)有關(guān)，因此氣流的

23、動力特性參數(shù)同樣也可以從進出孔口氣流的瞬態(tài)響應(yīng)中獲取。首先，利用圖5中的衰減曲線分別計算出進出口風速的一階和二階導數(shù)，再對非線性方程(5)或方程(6)的常系數(shù)進行參數(shù)擬合。由于系數(shù)中a、V0、pa、A0、q和P這幾個參數(shù)都可以由模型的幾何4 內(nèi)壓增益及參數(shù)確定增益函數(shù)|H()|是反映振動系統(tǒng)動力特性的最重要的函數(shù)之一，定義為響應(yīng)與激勵的傅立葉變換之比的模。它在振動理論中占有很重要的地位。本文討論的內(nèi)壓增益是指內(nèi)壓時程(內(nèi)壓系數(shù)時程)與孔口位置處外壓時程(外壓系數(shù)時程)的傅立葉變換之比的模。選擇內(nèi)壓增益作為數(shù)值模擬和風洞試驗的比較指標，可以有效分析數(shù)值模擬的準確性。對于風洞試驗只需分別測出結(jié)構(gòu)外

24、表面開孔部位附近的風壓系數(shù)和結(jié)構(gòu)內(nèi)部風壓系數(shù)的時程數(shù)據(jù)便可由下式獲得內(nèi)壓增益。fft(Cpi)|H()|= (8)fft(Cpe)內(nèi)壓增益|H(f)|尺寸和風場環(huán)境方便地獲得，因此方程中剩余的Le/c、CL、µeff/(rc)等系數(shù)便可隨之獲得。本文模型采用方程(5)，擬合所得常系數(shù)和對應(yīng)孔口動力特性參數(shù)為：aLeV0L=1.3×105e=0.19m (11)cA0pacaqV02CL=2.9×107CL=1.25 (12) 22(A0pa)µeffrPLeV0=5.1×105 22cA0pa頻率f/Hz圖6 CFD計算與風洞試驗的增益函數(shù)比較

25、 Fig.6 Comparison of gain functions obtained frommeasurements and CFDµeff對于CFD計算，首先將內(nèi)壓傳播方程線性化10&&+cC&+2C=2C (9) CpieqpiHHpiHHperc根據(jù)孔邊材料的物理性質(zhì)或相關(guān)材料試驗，例如取µeff等于層流粘性系數(shù)(1.8×10-5kg·m-1·s-1)，r=0.3mm。則由式(13)可計算得到本模型的c=0.5，再由式(11)可得Le=0.095m=L0+0.075m。=0.12kgm-2s-1 (13)式中

26、，ceq稱之為等效線性阻尼系數(shù)。再將方程(9)經(jīng)拉普拉斯變換得到內(nèi)壓增益的計算公式。22c2+eq|H()|=HHHHHH1/25 結(jié)論突然開孔結(jié)構(gòu)的瞬態(tài)內(nèi)壓傳遞函數(shù)已被公認為是一個二階非線性微分方程，方程各項系數(shù)隨孔口物理性質(zhì)的不同而不同。若能利用數(shù)值模擬獲得上述系數(shù)將大大節(jié)省研究成本。本文利用CFD數(shù)值計算模擬了突然開孔結(jié)構(gòu)的孔口氣流變化，并與(10)因此，只要CFD模型能準確計算出Helmholtz頻率HH=2fHH和等效線性阻尼系數(shù)，便能由式60 工 程 力 學on high buildings J. Journal of Xian Jiaotong University, 2001,

27、 35(5): 471474. (in Chinese)5 盧旦, 樓文娟, 陳勇, 唐錦春. 雙幕墻建筑通風性能的數(shù)值模擬研究J. 浙江大學學報, 2005, 39(1): 4650. Lu Dan, Lou Wenjuan, Chen Yong, Tang Jinchun. Study on ventilation in double-skin facades by numerical simulation J. Journal of Zhejiang University, 2005, 39(1): 4650. (in Chinese)6 Sharma R N, Richards P J

28、. Computational modeling ofthe transient response of building internal pressure to a sudden opening J. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 72: 149161.7 Liu H, Saathoff P J. Building internal pressure: suddenchang J. J. Eng. Mech. Div. ASCE, 1981, 107(Em2): 309321.8 Uematsu

29、 Y, Yamada M, Sasaki A. Wind-induced dynamicresponse and resultant load estimation for a flat long-span roof J. Journal of Wind Engineering and Industrial Aerodynamics, 1996, 65: 155166.9 Huang N E, Shen Z, Long S R. A new view of nonlinearwater waves: the Hilbert spectrum J. Ann Rev Fluid Mech, 199

30、9, 31: 34173457.10 Sharma R N, Richards P J. Computational modeling inthe prediction of building internal pressure gain function J. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 67 & 68: 815825.試驗結(jié)果相比較得出以下結(jié)論：(1) CFD數(shù)值模擬能對突然開孔結(jié)構(gòu)的孔口氣流變化進行準確地模擬。(2) 采用內(nèi)壓增益作為比較指標，可以直觀地檢驗數(shù)值模擬結(jié)果的正確性。(3

31、) 模擬得到的響應(yīng)曲線可以對瞬態(tài)內(nèi)壓傳遞方程中的各項系數(shù)進行尋優(yōu)擬合，從而獲得突然開孔時孔口氣流動力特性參數(shù)。 參考文獻：1 Liu H, Saathoff P J. Internal pressure and building safetyJ. J. Struct. Div. ASCE, 1982, 57(3): 22232234.2 布占宇. 大跨度柔性屋面結(jié)構(gòu)突然開孔時的風振響應(yīng)研究D. 杭州: 浙江大學, 2003.Bu Zhanyu. Study on wind-induced dynamic response of long-span flexible roof structure

32、s in the case of sudden opening D. Hangzhou: Zhejiang University, 2003. (in Chinese)3 Holmes J D. Mean and fluctuating internal pressuresinduced by wind C. Proc. 5th Int. Conf. on Wind Engineering. Fort Collins, USA, 1979: 435450.4 周莉, 席光. 高層建筑群風場的數(shù)值分析J. 西安交通大學學報, 2001, 35(5): 471474.Zhou Li, Xi Gua

33、ng. Numerical analysis of the wind field(上接第67頁)11 Wang Xu, Shen Yapeng. Inclusions of arbitrary shape inmagnetoelectroelastic composite materials J. I. J. Engineering Science, 2003, 41(1): 85102.12 Liu Jinxi, Liu Xianglin, Zhao Yongbin. Green's functionsfor anisotropic magnetoelectroelastic sol

34、ids with an elliptical cavity or a crack J. I. J. Engineering Science, 2001, 39(12): 14051418.13 Tian Wenye, Gabbert U. Multiple crack interactionproblem in magnetoelectroelastic solids J. European Journal of Mechanics A/Solids, 2004, 23: 599614.14 孫建亮, 周振功, 王彪. 功能梯度壓電壓磁材料中斷裂問題分析J. 力學學報, 2005, 37(1)

35、: 914.Sun Jianliang, Zhou Zhengong, Wang Biao. A permeable crack in functionally graded piezoelectric/piezomagnetic materials J. ACTA Mechanica Sinica, 2005, 37(1): 914. (in Chinese) 15 George R. Buchanan. Layered versus multiphasemagneto-electro-elastic composites J. Composites: Part B, 2004, 35: 413420.16 Garcia Lage R, Mota Soares C M. Layerwise partialmixed finite element analysis of magneto-electro-elastic plates J. Computers and Structures, 2004, 82: 12931301.17 Ding Haojiang, Jiang Aimin. A boundar