耦合場分析方法簡介、理論、分析方法和數(shù)值仿真

1、耦合場分析方法簡介、理論、分析方法和數(shù)值仿真課程主要內(nèi)容 引言 耦合場分析方法簡介 間接耦合分析方法 直接耦合分析方法 一些耦合場分析專題 （包括基本模型、解析和數(shù)值分析方法等） 力-熱耦合問題 流(氣)固耦合問題 (電)磁力耦合問題 風(fēng)沙電等耦合問題課程主要目的 了解和掌握多場問題的基本特征、耦合的性質(zhì)與意義； 接介紹幾類典型的多場耦合問題的基本模型與特征； 以幾類典型的多場耦合問題為例，介紹其分析的方法 和思路，從中體會和學(xué)習(xí)基本的方法； 用于自己的研究工作或者今后可能遇到的多場耦合 問題中 課程學(xué)習(xí)方式 課堂講授（主）+ 學(xué)生個人課后閱讀（輔） 講授內(nèi)容是一些具有典型的文章資料 （包括我

2、自己的相關(guān)專題的研究經(jīng)歷） 學(xué)生個人閱讀主要是提供的資料 + 個人從事的耦合問題的內(nèi)容 + 個人感興趣的耦合場方向； 課程成績（2部分組成平時考勤 + 期末 以提交課程報告的方式）多場耦合問題：(Muti-fields coupling problem) 研究兩個或者兩個以上的場通過相互作用而形成的物理（或力學(xué)）現(xiàn)象的問題。引 言 普遍存在于客觀世界 普遍存在于工程應(yīng)用領(lǐng)域 常見的耦合問題：結(jié)構(gòu)-熱耦合、流-固耦合、 結(jié)構(gòu)-電、結(jié)構(gòu)-磁耦合等越來越多的耦合問題：與智能材料關(guān)聯(lián)智能材料(Intelligent material、Smart material、daptive material an

3、d structure)是二十世紀(jì)90年代迅速發(fā)展起來的一類新型復(fù)合材料。智能材料就是指具有感知環(huán)境（包括內(nèi)環(huán)境和外環(huán)境）刺激，對之進行分析、處理、判斷，并采取一定的措施進行適度響應(yīng)的智能特征的材料。智能材料需具備以下內(nèi)涵：具有感知功能，能夠檢測并且可以識別外界(或者內(nèi)部)的刺激強度，如電、光、熱、應(yīng)力、應(yīng)變、化學(xué)、核輻射等； 具有驅(qū)動功能，能夠響應(yīng)外界變化； 能夠按照設(shè)定的方式選擇和控制響應(yīng)； 反應(yīng)比較靈敏、及時和恰當(dāng)； 當(dāng)外部刺激消除后，能夠迅速恢復(fù)到原始狀態(tài)。智能材料主要種類 形狀記憶合金（復(fù)合材料）； 電流變體和磁流變體(液體和彈性體、膠體)； 磁致伸縮材料（復(fù)合材料） ； 鐵電、壓電

4、陶瓷； 電致伸縮陶瓷； 智能材料系統(tǒng)（電、磁、溫度等敏感）； 光、電致變色材料等； 多場耦合作用下的材料功能研究是科學(xué)技術(shù)發(fā)展前沿！ 耦合的分類： 從耦合的空間屬性上分類 Felippa et al. Comput. meth. Appl. Mech. Engrg.，2001 區(qū)域耦合 整個區(qū)域或部分區(qū)域內(nèi)多場共存， 各場間無邊界。 如：結(jié)構(gòu)-熱、結(jié)構(gòu)-電(磁) 耦合 邊界耦合 各場間有明顯的邊界，場之間通過 邊界作用實現(xiàn)相互作用。 如：流-固耦合、空氣-彈性、壓電-結(jié)構(gòu) 多物理場：(Multiphysics)場尺度域“多場(Multi-field)”“多區(qū)域(Multi-domain)”“多

5、尺度(Multi-scale)”系統(tǒng)同時存在多個物理場的激勵和響應(yīng)系統(tǒng)的各個具有不同特征的連續(xù)體通過邊(交)界之間的相互作用系統(tǒng)中不同尺度下從微觀到宏觀行為的連續(xù)一致跨越“一個”、或“多個”“納尺度”、“微尺度”、“宏觀尺度”耦合問題領(lǐng)域幾個發(fā)展方向（“十一五”學(xué)科發(fā)展規(guī)劃） 力-電-磁-熱耦合場的分析理論； 智能材料的本構(gòu)關(guān)系； 智能結(jié)構(gòu)動力學(xué)與主被動控制； 耦合場的破壞力學(xué)、失效機理 與智能器件的可靠性； 風(fēng)-沙耦合、風(fēng)沙電耦合問題 凍土、巖石，應(yīng)力場-溫度場-流場 -空氣泡耦合等 蘭州大學(xué)電磁固體力學(xué)研究組、風(fēng)沙物理研究組間接耦合 & 直接耦合方法 間接耦合方法 ，或稱順序耦合、序貫耦合

6、方法 按照順序進行兩次或更多次的相關(guān)場分析。 物理或變形場 1物理場 2物理場 3耦合場分析方法例：結(jié)構(gòu)-熱分析溫度場分析熱載荷結(jié)構(gòu)變形場分析數(shù)學(xué)描述（以兩個場耦合為例） 物理或力學(xué)變形場1： 場變量u 物理或力學(xué)變形場2： 場變量輸入初始值： 01 0獲得場1的解： ui作為初始值輸入獲得場2的解： 1 直接耦合方法 兩場或更多場的同時求解，以獲得耦合場的解。 物理或變形場 1物理場 2物理場 3+耦合場的解數(shù)學(xué)描述（以兩個場耦合為例） 物理或力學(xué)變形場1： 場變量u 物理或力學(xué)變形場2： 場變量同時獲得場變量：u、？兩種耦合分析方法的比較 直接耦合方法 間接耦合方法 迭代思想 分場求解、方

7、程階數(shù)低 適合非線性程度不高的問題 每個場分析中均采用收斂條件 可能出現(xiàn)結(jié)果發(fā)散現(xiàn)象 逆算子思想 “合場”求解、方程階數(shù)高 適合高度非線性問題 “合場”方程建立困難 高維非線性問題帶來的困難理論上講，不受問題限制，適合任何耦合場分析如：壓電-結(jié)構(gòu)耦合、流動-熱傳導(dǎo)耦合、電路-電磁場耦合等耦合場的分析方法： 解耦方法 順序求解各個物理場或者力學(xué)變形場，將獲得了上一個場的相關(guān)信息后代入下一個場進行分析，最后獲得多場作用下的總效果。 單向非雙向、考慮作用但非相互作用與影響 并非真正的耦合，意義？ 實際上我們熟悉了太多這樣的問題：溫度應(yīng)力問題、早期的電磁結(jié)構(gòu)變形分析、小變形、低溫、低頻、低電磁場下結(jié)構(gòu)

8、分析等 可以給出一些解析解，可作為考慮耦合效應(yīng)的考據(jù) 解耦單向分析思路考慮了雙向的作用與影響就是順序耦合思想耦合問題的求解（間接、直接耦合分析） 主要適合解耦場分析、低維、低非線性 可在某些條件下的線性化問題分析中 解析、半解析求解耦合問題 數(shù)值求解耦合問題 目前的主要手段，適合多個場分析，稍高維、非線性 分為網(wǎng)格方法（有限元法、邊界元法、有限差分法、有限體積法等）和無網(wǎng)格方法（再生核質(zhì)子方法、有限點方法、MPLG法等） 數(shù)值仿真軟件 商業(yè)軟件：具有一定的耦合場分析功能 FEMLAB：基于偏微分方程基礎(chǔ)的軟件，最新V3.2,可求解聲場、擴散、電磁場、流體力學(xué)、結(jié)構(gòu)力學(xué)問題或耦合問題； ANSY

9、S:最初為解決固體力學(xué)和結(jié)構(gòu)力學(xué)問題，最新V10.0,陸續(xù)加入了對流場、聲場、熱場、電磁場的仿真功能，以及多場耦合的仿真算法； MSC.DYTRAN:高度非線性、流體-結(jié)構(gòu)耦合、瞬態(tài)動力響應(yīng)問題仿真； ALGOR: 功能包括結(jié)構(gòu)，流體，熱，電磁分析以及目前主流有限元分析軟件中最為便捷的多物理場耦合分析:流固耦合分析和熱結(jié)構(gòu)耦合分析,最新V14； ABQAS: 結(jié)構(gòu)（應(yīng)力/位移）問題，以及工程領(lǐng)域的熱傳導(dǎo)、質(zhì)量擴散熱電耦合分析、聲學(xué)分析、巖土力學(xué)分析（流體滲透/應(yīng)力耦合分析）及壓電介質(zhì)分析等。 開發(fā)耦合分析模塊，或者商業(yè)軟件的二次開發(fā)一些耦合場分析實例1. 力 - 磁、力 - 磁 - 熱 耦合問

10、題 區(qū)域耦合問題 解析解法 數(shù)值解法（有限元）間接耦合分析方法 多重非線性迭代技術(shù)Background & ObjectiveApplications: Magnetic fusion, Energy storage device, Magnetohydrodynamic system (MHD), Magnetic forming, Magnetically levitated vehicles (MLV) Magnetic guns or cannons in military field, Nuclear-magnetic-resonance measurement (NMR) for

11、medical use Problems Induced & Objectives:Stress in electromagnetic structures induced by electromagnetic forcesMagneto-elastic stabilityMechanics behaviour of electromagnetic structures under coupled multi-fields, such as magnetic, thermal, fluid fields and so on (1) Magnetic energy of ME system (M

12、agnetization nonlinearity)Mathematic Modeling (板殼的力-磁變分理論,多非線性) Magnetoelastic generalized variational principle(2) Strain energy of plate (Geometrical nonlinearity)(3) Total generalized energy of ME system (4) Magnetoelastic generalized variational principle Magnetic FieldGoverning equations& Bound

13、ary conditionsMechanics Deformation FieldGoverning equations &Boundary conditions for plates(5) Equivalent magnetic forces exerted on SFM platesExplanation: Transformation from the magnetic energyto the mechanical energy of the system. Magnetic FieldMechanics Deformation FieldNumerical Method - Coup

14、led FEM for Multi-fields FEM for Magnetic FieldN-R Method for nonlinearity of MFFEM for Deformation FieldN-R Method for nonlinearity of MDFIteration Method for nonlinearity of coupling fields: SolutionsNumerical Simulation Results (1) Linear magnetization and linear deformation for ME systemNote: i)

15、 For cantilevered SFM plate in transverse magnetic field. ii) Edge effect of magnetic field included.Zheng XJ,Zhou YH,Wang X, et al ASCE J. Eng. Mech. 1999(2) Nonlinear magnetization and linear deformation for ME systemSFM simply supported plate in oblique magnetic field: Wmax vs. B0.SFM simply supp

16、orted plate in oblique magnetic field: Wmax vs. .Zheng XJ, Wang X, INT. J. Solids Struct. 2001(3) Linear magnetization and nonlinear deformation for ME systemZheng XJ, Wang X, ASCE J. Eng. Mech. 2003Effect of incident angle on B0cr for SFM cantilevered plate(b) Post-buckling(a) Pre-buckling(4) Mathe

17、matic Modeling and simulation For SFM ShellsStrain energy of shell對已有實驗的模擬數(shù)值模型Zheng XJ, Wang X, INT. J. Solids Struct. 2003(1) Magnetic energy of MTE system (Magnetization nonlinearity)(2) Total mechanical energy of thermoelasticity for MTE system Mathematic Modeling (廣義磁熱彈性變分理論) Magneto-thermo-elas

18、tic generalized variational principle (3) Heat potential energy of thermal flux of MTE system (4) Functional of total generalized energy of MTE system (5) Magneto-thermo-elastic generalized variational principle Magnetic FieldGoverning equations& Boundary conditionsDeformation FieldGoverning equatio

19、ns& Boundary conditionsThermal FieldGoverning equations& Boundary conditionsAnalysis (Magneto-thermo-elastic buckling of SFM plate ) 解析分析磁熱彈性板的屈曲問題 Equations : Equations for magnetic field; T: Equations for thermal field; w: Bending equation of plate: Solutions For simply supported rectangular SFM p

20、lates(without edge effect)磁場溫度場結(jié)構(gòu)變形場+線性化、攝動理論 BucklingCase (i). Magneto-elasticity:Case (ii). Thermo-elasticity:Case (iii). Magneto-thermo-elasticity:Wang X, Zhou YH, Zheng XJ, Int. J. Eng. Sci. 2002For rectangular SFM plates (with edge effect)+磁 場Simulations (5) Repeat (1)(4)，until the following co

21、nditions (1) (2) (3) (4) Wang X, Zheng XJ, Lee J.S., INT. J. Solids Struct. 2003一些耦合場分析實例2. 力 - 磁耦合動力學(xué)問題 力-磁耦合動力學(xué)模型與數(shù)值分析 復(fù)雜動力學(xué)行為:非線性、磁阻尼、混沌 Xingzhe Wang, et al. ASCE Journal of Engineering Mechanics, 2006, 132(4):422-428 Xingzhe Wang, et al. Int J of Mechanical Sciences, 2006,48(8):889-898 Xingzhe W

22、ang, Int Conference on enhancement and promotion of Computational Methods in Engineering Science and Mechanics, 2006, Aug, Changchun,China 一些耦合場分析實例3. 空氣-彈性、空氣-彈性-控制耦合問題 邊界耦合問題 半解析半數(shù)值解法Background & Objective High density and high speed HDD, VCD/DVD, Floppy Disk. 應(yīng)用背景: Flexible rotating blades, gas t

23、urbines, circular saws 容量(磁盤密度): 增長了25M倍， 100% /每年驅(qū)動電機： 幾百轉(zhuǎn)/分鐘 7500轉(zhuǎn)/分鐘 上萬尺寸： 24英寸 1.0 英寸盤片： 24片 2-3片、單片 現(xiàn)代高密磁盤發(fā)展 磁盤工業(yè)與設(shè)計中(HDD)的力學(xué)問題磁頭懸臂的振動、動力穩(wěn)定性；讀寫磁頭的懸浮、定位與控制；磁記錄介質(zhì)表面摩擦學(xué)；磁盤的噪聲與控制；高速旋轉(zhuǎn)磁盤空氣彈性失穩(wěn)顫振，及其控制.Rotating Disk Flutter（顫振）? Hydrodynamic instability caused by aeroelastic coupling between rotating

24、disk and surrounding airflow. Critical speed for disk flutter Flutter Speed (臨界旋轉(zhuǎn)速度)旋轉(zhuǎn)振動圓盤具有穩(wěn)定性:1) 屈曲失穩(wěn)(行波頻率之一 等于零)2) 顫振失穩(wěn)(負阻尼)目前盤片顫振抑制研究 By enhancing the disk stiffness Flutter speedBy designing the base casting Heo et al. (2000) 實驗Reduce disk-rim to shroud gapSmooth the shroudBy employing air sque

25、eze film Bittner and Shen (1999), Ono and Maeda (2000)， Deeyiengyang and Ono (2001).實驗 本文研究 磁盤顫振穩(wěn)定性分析、臨界轉(zhuǎn)速的預(yù)測； 提出一種主動控制方式，抑制磁盤顫振失穩(wěn)； 進行相關(guān)實驗并給出實驗結(jié)果； 與理論預(yù)測和數(shù)值模擬結(jié)果進行對比。THEORETICAL MODELING （理論模型）Disk Actuator(Speaker)Disk SensorEnclosureFeature: Indirect, non-contact method 非接觸控制 封閉或周邊開口問題描述 Air-couplin

26、g force, Acoustic pressure loading, Control acoustic force.Description of rotating disk vibrationBoundary Conditiona) At the clamped edge: b) At the free edge: Mathematical Modeling 數(shù)學(xué)模型 空氣壓力粘性旋轉(zhuǎn)流體壓力控制力旋轉(zhuǎn)磁盤振動(1). Aerodynamic Force Induced By Rotating Disk-Airflow Coupling -Rotating Damping Model C:

27、Damping coefficientd/ : Rotation speed ratioDescription of loadings of system(2). Acoustic Force Induced By Acoustic-Structure Coupling The boundary conditions, match conditions on the disk surface and at the clearance between the disk rim and enclosure: Where a is the acoustic velocity potential, a

28、nd is governed by聲場(3). Acoustic Control Force Induced By ActuatorWhere c is the acoustic velocity potential, and is governed by控制聲場Description of loadings of system+(1)(2)(3)兩個問題：（1）旋轉(zhuǎn)盤片的空氣彈性動力學(xué)行為？（2）顫振反饋控制的實施？Yasuda et al. (JSME Int J. 1992): Aerodynamic force, “damping” and “l(fā)ift” ，the ratio prop

29、ortional to . Kim et al. (J Sound Vib. 2000),Hansen et al.(J Fluids Struct. 2001) in an aeroelastic model to include a series of parameters determined for each mode at each rotation speed. Wang X,Huang,X, Acta Mechanica Sinica, 2006Renshaw et al. (J Sound Vib. 1994) taking into account the air coupl

30、ing and examining the eigenvalues of the whole disk systems. ? Solutions 半解析半數(shù)值求解 Characteristics Difficulties arising from couplings All equations should be solved synchronouslyTransverse DisplacementDisturbed Acoustic FieldsCouplingsMatching conditions of velocitiesOn the surface of disk:and the s

31、urfaces of actuators:Equation of disk vibration Solutions 半解析半數(shù)值求解 假設(shè)含有參數(shù)的變形場、聲場的解 satisfy partial boundary conditions or governing equations (m,n): Disk vibration mode m =: Nodal circle number n =: Nodal diameter number =: EigenvalueGalerkins Method - free vibration of rotating disk - aerodynamic,

32、acoustic, control forces - unknown coefficient matrix Natural FrequencyRe() DampingIm() 0 0Stable Unstable (Flutter)Nontrivial Solution Eigenvalues (FTW, BTW)Material & geometric properties DiskDensity, d (Kg/m3)7.8X103Outer radius, ro (m)0.178Clamping ratio, 0.3Thickness, h (m)0.775Youngs modulus,

33、E (GPa)200Possions ratio, 0.3EnclosureRadius, re/ro1.2Height, ze/ro0.5AirflowDensity, a (Kg/m3)1.21Speed of sound, a (m/s)340 Simulation results The properties of disk same as the ones used in DAngelo et als experiment.(b) Imaginary part of eigenvalue or dampingVerifications & Observations of Disk F

34、lutter(a) Real part of eigenvalue or mode frequencyCritical speed(C=0.02; d/=2/3)The feedback control method proposed can suppress disk flutter effectively The actuating system has a large operation region and therefore it is robust Huang,X, Wang X，J. Fluid & Structures, 2004Control Performance & Op

35、timizationCase (a).One Piezo-patchCase (b).Two Piezo-patchesCase (c).Three Piezo-patchesSchematic diagram of Piezo-patch(es) arrangement on the upper cover plate surface Case (a) - One piezo-patch actuator Control performance for one actuator with r1=0.7,r2=1.0, and =100(a) Damping vs. Rotation spee

36、d(b) Stability Map in G- planeEffect of actuator size on control performance for one actuator(a) Gmin vs. Radial width (Fixed =100)(b) Gmin vs. Sector angle (Fixed r1=0.7,r2=1.0)Case (a) Disable for control: =1200 for mode (0,3); =900, 1800 for mode (0,4); =720, 1440, 2160 for mode (0,5). Actuatorr2

37、 -r1Control performance for two actuators(a) Gmin vs. Sector angle (Fixed =900)_ Fixed relative angleChange sector angle Case (b) - Two piezo-patches (G1= G2, 1=2)Control performance for two actuators(b) Gmin vs. Relative angle ( =400/2)_ By dividing one patch into two patches and arranging with a p

38、roper relative angle, the control gains are smaller.Fixed sector angle Relative angleChange relative angleCase (b)Case (c) - Three piezo-patches (G1 = G2 = G3)Control performance for three actuators(a) Gmin vs. Sector angle (1=2=3) Flutter control for modes (0,4) and (0,5) is disable.Change sector a

39、ngleCase (c)Control performance for three actuators(b) Gmin vs. Relative phase shift 123 Fixed locations. Fixed size =400/3. Different phase shift: 1, 2=1+123, 3=2+123. Wang X, Huang,X, AIAA Journal, 2006程序設(shè)計相關(guān)問題1. 程序設(shè)計的一般原則 針對具體問題的特點選擇合適的數(shù)值方法 問題的復(fù)雜程度，數(shù)學(xué)描述的可選途徑 理論知識 各種方法的優(yōu)缺點要有認(rèn)識 閱讀與學(xué)習(xí) 多種方法的聯(lián)合使用 比較與思

40、考 可優(yōu)先考慮自己所熟悉掌握的方法 自身優(yōu)勢、事半功倍 使用商業(yè)軟件或二次開發(fā)、或完全編寫代碼 直接使用了解商業(yè)軟件的功能，特別是耦合場分析能力 二次開發(fā)商業(yè)軟件的開放性、可開發(fā)性 編寫代碼各個物理或者力學(xué)場分析方法的實現(xiàn)、零散代碼 的利用與集成 相關(guān)數(shù)值算法的掌握 基本的：矩陣運算、微積分運算、排序、特殊函數(shù)等 特別的： 迭代算法（收斂性、穩(wěn)定性） 特征值問題（線性、非線性） 動力學(xué)問題 （Willson-、Newmark方法等） 非線性方程算法（ Newton 法、Newton-Raphson法、 最速下降法、共軛梯度法等） 圖形圖像處理等 程序設(shè)計思想 分場、依次對涉及的多長問題進行程序設(shè)計，確定輸入量、 輸出量 采用合適的算法處理場-場耦合，迭代算法的收斂條件等 聯(lián)系各個場之間的數(shù)據(jù)傳輸與處理 最終輸出結(jié)果的表征（圖像圖像等） 程序調(diào)試與結(jié)果的可靠性判斷 最耗時、長期的、也是計算能力的體現(xiàn) 合適的考題驗證數(shù)值編碼的正確性 局部驗證每一個子程序、每一個場的驗證：解析解、退化 結(jié)果、對稱性、不同的軟件、算法求解統(tǒng)一問題比較 等、結(jié)果穩(wěn)定性的驗證（不同的單元剖分、節(jié)點數(shù)目 的變化、邊界條件的變化、空間和時間步長