隨著制造業數字化時代到來，其軟件的重要性日趨顯現，并起著舉足輕重的作用，計算科學發展成為影響國家利益與國家安全的戰略性問題。一些國家將計算機建模與仿真列為優先發展的服務于國家利益的關鍵技術。有國內學者指出應將自主CAE軟件的開發提高到戰略發展高度。國內外開發了大量的通用和專用的CAE軟件，哈爾濱工程大學自2003年起，啟動CFD軟件開發，最初研究目標是針對航空發動機燃燒室內的流動、傳熱與燃燒過程分析，柴油機缸內工作過程模擬，隨著工作的深入和發展需要，發展了船舶水動力學分析模塊。后續在結構聲輻射、水動力噪聲以及結構聲耦合方法和軟件模塊開發方面進行初步工作。

對于數值計算程序開發人員來說，劃分求解域網格，得到求解域網格信息，例如：網格坐標，網格組成網格面系列等。這是一項具有重復性、復雜性的工作。網格生成技術及軟件開發仍然存在挑戰，國內相關單位開展了卓有成效的工作，但尚未形成通用、穩定和可靠的網格生成軟件，也是當前CAE發展的短板之一。有經驗的研究人員能夠自己寫出各種網格信息，但只限于幾何結構簡單，網格數量較少的場合；對于專業人員進行CFD求解器開發的研究人員，自己寫出復雜問題的非結構化網格信息并保證不出錯誤是有一定難度的。因此采用一種折中的方法，團隊最初開發了CGNS（CFD General Notation System）標準格式文件的輸入輸出，自編程序讀取CGNS格式文件的網格信息，處理輸出，應用于本課題組自行開發的求解器，即提供一個標準文件接口。可以采用通用的網格生成器，如ICEM CFD生成CGNS格式網格文件。為了進一步支持主流的CFD軟件的前處理格式文件，提高軟件的通用性，先后開發了前處理軟件GAMBIT生成mesh格式文件和Star-CCM+軟件生成網格文件。

CGNS 系統是一種良好的標準，核心求解器研發與商用前處理和后處理軟件聯系起來。核心求解器開發人員采用CGNS系統作為標準接口，則可以解除網格生成和計算結果后處理問題。作者也采用了CGNS標準作為輸入輸出格式，讀取非結構網格和混合網格，包括二維三角形或四邊形網格，三維四面體、六面體、金字塔形網格。

隨著網格自適應和直角網格技術發展，如商用軟件Star-CCM+和開源軟件Openfoam等廣泛應用，對多邊形和多面體網格前后處理的需求越來越大。作者開發的軟件采用了非結構化網格有限體積法，而且基于面循環的方法形成系數矩陣，該方法適用于多邊形和多面體網格。問題在于目前的多邊形網格和多面體網格生成器難以轉化為CGNS標準格式，作者開發了直接讀取Openfoam網格文件的接口程序，新版tecplot在原有基于單元的結果輸出基礎上，新增了基于面輸出方法，包括多邊形（FEPOLYGON）和多面體單元（FEPOLYHEDRON）。

3.1 基本方程

作者開發GTEA軟件前提是連續介質假設，連續介質力學思想下統一流體力學、固體力學、聲波動控制方程，詳細過程可以參考相關文獻。軟件分為計算流體力學、結構應力波求解、流聲耦合和聲固耦合四大功能模塊，采用非結構化網格離散控制方程，基本思想如圖1所示。

3.2 并行計算流程

GTEA基本流程圖如圖2所示，在串行計算程序的基礎上，進行了并行設計。與串行程序相比，并行計算的增加分為：區域分解與結果重建兩個部分。程序開始運行后，讀取網格，根據計算環境判斷是否存在多個進程，存在多個進程時在主進程中根據進程數目劃分網格并映射到各進程，讀取設置參數后進行離散過程形成線性方程組，利用并行解法進行計算，方程組循環收斂后，對變量通過通信后進行更新，為下一步計算準備。在計算過程收斂后輸出結果，對于存在多進程時，會增加輸出子區計算結果，并在主進程中調用結果重建的環節。程序線性方程并行計算過程中，各相鄰區域之間通過MPI進行通信，線性方程組求解過程需要全局通信。

                                                           圖2 GTEA程序計算基本流程圖

正如前面提及存在多進程計算時，需要將網格根據進程數目劃分成相應的子區。本文調用了Metis軟件系列的partdmesh工具在主進程區域分解部分，Metis軟件包含了不同的工具，如partnmesh和partdmesh，對不同的進程數二者在分區時效率不同。為了支持任意混合網格，作者提出了先將網格文件轉化成圖形文件，然后調用kmetis劃分圖形。兩者的區別是：前者以網格節點為節點轉化成圖形文件，而后者則以單元為節點。選擇partdmesh是因為它可以使各相鄰區域間通信兩最小，而且各進程間負載較partnmesh平衡。采用了基于圖論方法劃分網格的partnmesh工具。

區域分解模塊工作過程中，程序會通過METIS工具根據進程數將計算域分解，采用了主從模式，區域分解模塊在主程序調用，對區域內的每個單元進行標識所在的進程標識，據此形成子區內部單元與ghost單元。因為子區域網格的構建依賴于通信模式的選擇，在介紹子區域網格構造之前，先對通信模式做一下說明。通信模式主要有兩種，以交界面和相鄰單元為基礎的交換模式。

                                                                           a分解前區域

                                                                            b通信模式1

                                                                         c 通信模式2

                                                        圖3 并行計算區域分解和通信模式示意圖

4.1 自然對流與輻射傳熱耦合

在計算流體力學標準算例測試基礎上，對軟件在流動與輻射傳熱方面的應用進行了測試。主要算例為自然對流與輻射耦合分析。重點分析傾斜角a對流體流動和傳熱的影響。參數設置為Ra數為5×106，壁面發射率為e1,2,3,4 = 1，四種情況：t = 0 和t = 1, w = 0, 0.5, 1。腔體的長寬比AR是常數1。傾斜角以每步30度從-60度變化到90度。結果見圖4。

圖4a, b描述了傾斜角為負時，等溫線隨著傾斜角和散射反照率的改變的變化，正角度部分在圖4d、f中描述。當t = 0時，無論傾斜角a如何變化，方腔中心部分的等溫線都垂直于重力方向。對于負角度，等溫壁面附近的溫度梯度相比于正角度要大。從圖4 a, b可看出，w=0和0.5時的等溫線結構非常相似；而且，當w從0.5變化到1時，介質的散射明顯改變了方腔內的熱分布。然而，對于正傾斜角度，方腔內存在熱分層并且中心區域的等溫線是水平的。實際上，這個現象意味著少量的熱傳遞。從圖4 d, f可看出，散射反照率對方腔中等溫線的結構略有影響。

                                                                            (b)

                                                                             (d)

                                                                          (f)

                              圖 4 對t = 0和t = 1, w = 0, 0.5, 1，傾斜角對等溫線的影響：

                             (1) t = 0; (2) t = 1, w = 0; (3) t = 1, w = 0.5; (4) t = 1, w = 1

4.2 發動機缸內燃油噴射過程仿真

       無論是航空發動機還是柴油機，燃油霧化過程模擬是發動機流動、燃燒和傳熱模擬的重要環節之一，下面給出了燃油霧化過程與實驗對比結果，驗證軟件在燃油霧化模型的準確性。

                                                  圖5 不同液滴破碎模型下，燃油噴射模擬結果

4.3柴油機缸內工作過程模擬

為進一步說明，圖6所示為層動區域及滑移邊界的應用示意圖。紅線為靜止的不匹配網格邊界，主要有兩個，一個是氣道區域與氣閥的區域的銜接處，一個是燃燒室區域與氣缸區域的銜接處；藍線為滑移邊界，主要在氣閥內部區域及外部區域的交接處以及氣閥外部區域與氣缸區域的交界處。彩色區域除區域2外均為層動區域，其中區域2為跟隨氣閥運動區域，區域1為適應閥桿運動的層動區域，區域3為適應氣閥上表面運動的層動區域，區域5為適應活塞運動的氣缸層動區域，區域4為適應氣閥底部及活塞運動的層動區域。計算過程不同時刻的網格和流場如圖7和圖8所示。

                                                            圖6層動區域及滑移邊界的應用示意圖

                                                                  圖7 動網格運動示意圖

                                                                    圖8 缸內流動跡線圖

4.4 燃氣輪機燃燒過程仿真

利用GTEA軟件對燃燒室內的流動和燃燒進行了數值模擬，數值模擬斜切徑向旋流器環形燃燒室內冷、熱態流場，并對所得計算結果進行分析。部分的計算結果如圖9，計算所得燃燒室三維整體熱態流場內通過旋流器中心截面的速度矢量圖。由放大圖A圖可知，與雙級軸向旋流器燃燒室一樣，在擴壓器的上下突擴段有明顯的旋渦，這兩個旋渦作為氣動壁面根據進口狀態的變化可自動進行調整，保證火焰筒進口氣流穩定，旋流器出口高速旋轉射流與主燃孔射流相互作用，在火焰筒頭部產生了上下兩個較強的旋渦，形成中心回流區如放大圖B所示。C區放大圖可以看到，主燃孔后，氣流逐漸平穩，冷卻孔形成局部的高速區。

                                                                圖9 中心截面速度分布圖

圖9為計算所得中心截面氣流溫度場，從圖中可以看到，從火焰筒氣膜冷卻孔流出的冷卻空氣對火焰筒壁面起到了很好的冷卻作用，同時在主燃孔后，由于未燃盡燃氣在火焰筒主燃孔后低壓區繼續燃燒，因此主燃孔后火焰筒壁面溫度比周圍的高。

4.5 水下航行體SUBOFF水動力仿真

確立的潛艇模型是以SUBOFF 潛艇模型為基礎進行建立的。由美國國防預研規劃署（DARPA）確立的SUBOFF 項目，SUBOFF 潛艇模型的艇體部分是用回轉體進行模型確立的。關于SUBOFF 項目的提出，為了方便我們對潛艇模型進行模擬計算，分析其水動力性能和潛艇周圍流場。自從1989 年美國國防預研規劃署（DARPA）提出這個模型，已經掌握了大量的實驗數據，為很多國內外的學者提供了研究驗證的對象。

                                                    圖10 SUBOFF模型流向速度速度和壓力分布圖

                                                      圖11  SUBOFF模型中縱剖面上壓力系數分布

4.6結構砰擊水動力仿真

計算采用了文獻中模型，如圖12所示，計算域大小為3.22m×1m×1m，障礙物位于矩形水箱底部，水體前側，大小為0.16m×0.4m×0.16m，距離布置水體側壁面2.48m。初始水體大小為1.228m×1.0m×0.55m。在矩形水箱底部H4(0.582,0.5)設置波高儀監測自由面高度，在障礙物表面P3(2.4,0.475,0.10)設置壓力傳感器監測壓力大小。

                                                                    圖12 計算模型示意圖

計算模型網格取為結構化網格，網格數為161×50×50。圖13為自由面高度監測點H4(0.582,0.5)處自由面高度隨時間變化曲線，圖14壓力監測點P3(2.4,0.475,0.10)處的壓力隨時間變化曲線，可以看到計算結果和實驗結果吻合良好。

                  圖13 H4監測點處自由面高度隨時間變化曲線           圖14 P3處的壓力隨時間變化曲線

4.7復雜消聲器聲學性能仿真

如圖15該消聲器的出口管改到側面，則該消聲器結構的進出口邊界由軸對稱變為只有一個對稱面，這就導致計算網格會成倍的增加，同時為測試FVM方法在大型計算中的能力，劃分網格數為778579，利用8進程進行仿真。對于在本節提到的計算機上，FEM是難以完成的，而FVM方法由于采用顯格式推進，并且不需要計算大型的剛度矩陣，所以可以比較快捷的處理較大型計算問題。為出口管位置不同時消聲器的傳遞損失對比，可以看出，當出口改為側面后在共振頻率處出現多個峰值：

第一個共振峰是由于側面出口管與消聲器端部之間形成一個聲學共振腔引起的；

后面的幾個共振峰是由于側面出口使消聲器變為非對稱結構而激發周向模態，消聲器的高階模態提前激起而形成的，出口改到側面后消聲器的整體聲學性能得到改善，特別是1500Hz以下的聲學性能得到明顯提高。

                                            圖15 復雜消聲器網格圖與消聲器傳遞損失比較

4.8 流場動力噪聲

計算模型和網格劃分見圖16，采用相同的模型和網格，在圓柱壁面附近進行加密。圓柱直徑為D=0.019m，流場馬赫數Ma=0.2，雷諾數Re=200，介質為空氣。時間步長為dt=2×10-5s，總計算時間為t=0.15s，當不可壓流場計算趨于周期變化時，即當t=0.05 s時，開始計算聲場。進行流場計算時，左邊界為速度進口，右邊界為壓力出口，上下為對稱邊界，域內正方形為內部面，在聲場計算中，域內正方形到外邊界為PML區域。監測點設于坐標為A（0 ，10D）、B（0 ，35D）、C（0，44.9D）和D（0 ，45.1D）。

                                                              圖16 計算模型和網格劃分

                                                        圖17 瞬時聲壓云圖和y軸的聲壓分布

圖17（a）為計算區域聲壓分布云圖，從圖中可以看出，圓柱繞流噪聲是個偶極子聲源，PML區域吸收良好，圖17（b）是y軸的聲壓分布，可以看出圓柱上下方聲壓幅值相反。

4.9結構聲耦合算例

考慮如圖18所示的大壩-水庫結構聲耦合系統。大壩的介質屬性為：彈性模量E = 3.437×109 Pa，泊松比μ = 0.25，密度ρs = 2000kg/m3，縱波波速為cp = 1436m/s；水的介質屬性為：聲速ca = 1436m/s，密度ρa = 1000kg/m3。水面高度H = 50m，均勻分布的正弦載荷f(t) = 200×sin(18t) N/m作用在大壩頂部。大壩底部采用固支邊界條件，水庫上側水面為絕對軟邊界，下側固體壁面為絕對硬邊界，右側為吸收邊界。在結構子域、聲學子域均采用顯格式求解。

                                                      圖18 大壩-水庫系統和計算網格

(1)    四邊形單元處理改進前后結果對比分析

網格模型如圖14所示，改進前采用常應變型四邊形單元模型，改進后采用雙線性四邊形單元求解結構聲耦合問題。邊界上網格尺寸為5m，與文獻相同，結構子域中網格的最大邊長Lmax ≈ 9.038m，最小邊長為Lmin ≈ 2.744m，聲學，子域中Lmax = Lmin = 5m。基于顯-顯耦合格式求解，依據穩定性條件，結構子域的時間步長需要滿足?ts ≤ 1.91×10-3s，聲學子域的時間步長需要滿足?ta ≤3.48×10-3s，最終結構子域與聲學子域取相同時間步長為1.5×10-3s。

                      (a)A點位移uy                                  (b)B點聲壓               圖19 四邊形單元計算結果

圖所示為A點位移uy的時間響應與B點壓力的時間響應。將計算結果與文獻中FEM程序的計算結果對比如圖19所示，發現對于該非點源問題，改進處理前的四邊形單元不能得到正確的結果，而改進后的計算結果與文獻結果吻合更好，說明對四邊形單元的雙線性改進處理，提高了傳統常應變單元方法在處理結構聲耦合問題時的正確性。

4.10 結構模態分析

潛艇是一種既能在水上又能在水下航行的艦艇，根據其性能指標需求其外形有不同構形，但主體結構基本相似。對于單層殼潛艇，其耐壓殼體可以看做是由圓柱殼、圓錐殼和半球殼組合組合結構，本節將采用此簡化結構分析艇體主體結構振動特性，計算結果與實驗值及FEM結果進行對比，FVM與FEM分析采用同一套網格。簡化結構實物模型及FVM分析模型如下圖所示：

                                             圖20簡化結構實物模型                        圖21仿真模型

                                                                    圖22 前6階模態結構

4.11 水下結構聲固耦合仿真分析

考慮殼體彈性對聲場的影響，殼體厚度為0.015m，結構為鋼，介質屬性為密度ρs = 7700.0kg/m3，楊氏模量E = 1.44×1011Pa，泊松比v = 0.35，假設為平面應變問題，縱波波速為cp = 5478.5m/s。殼體結構子域908個三角形單元劃分，最小邊長為Lmin ≈ 0.0106m，依據穩定性條件，選擇的時間步長Δtso = 1.36×10-6s；聲學子域采用62084個三角形單元劃分，最小邊長為Lmin ≈ 0.0144m，依據穩定性條件，選擇的時間步長Δtao = 8.22×10-6s；最終采用不同時間步長技術?ta = 5?ts = 5.0×10-6s。

圖23為水下彈性結構時，不同時刻聲壓云圖分布，可以看到點聲源的輻射場到達殼體后，在t = 0.0010s時已經激起結構振動并產生聲輻射，t = 0.0020s后在水下結構頭部可以明顯看到結構振動產生的輻射聲場，這主要是由于聲源在尾部激起結構振動，而結構中縱波的波速遠大于水中的聲速，由此產生波傳遞的速度差，結構中縱波先到達頭部，并且振動產生新的聲輻射。

                                                        圖23 彈性結構不同時刻聲壓(Pa)云圖

5、 結論

      文章介紹了明平劍教授團隊圍繞多物理場耦合問題進行計算研發的基本思想，計算軟件開發的一些實際問題及解決方法，包括介紹了網格格式及并行計算問題等，基于非結構化網格有限體積方法研究可以參見明平劍教授及團隊成員發表的文章。