常亮　段世慧　王立凱　羅利龍

摘要：

針對飛機設計精細化數(shù)值分析模型自由度已經(jīng)達到億級，對高性能計算的要求也越來越高的問題，圍繞大規(guī)模并行計算環(huán)境下結構分析和優(yōu)化的若干關鍵問題，研究滿足高性能計算體系特點的區(qū)域分解并行算法、超大規(guī)模結構變量敏度高效求解和結構非線性振動特性求解等關鍵技術.對國產(chǎn)CAE軟件HAJIF進行并行化改造，初步實現(xiàn)基于最大航程的氣動結構綜合優(yōu)化設計和基于精細化模型的復合材料機翼綜合優(yōu)化設計.HAJIF的計算效率和精度得到明顯提高.

關鍵詞：

區(qū)域分解算法； 超大規(guī)模結構變量； 解析敏度； HAJIF

中圖分類號： V221

文獻標志碼： B

0 引 言

進入21世紀以來，以高性能計算機為基礎的計算科學得到長足的發(fā)展，在航空航天、地震預測、核爆模擬等許多工業(yè)領域得到應用，高性能計算水平的高低已經(jīng)成為衡量企業(yè)、研究機構科技水平的重要標志，更是檢驗一個國家高科技能力的試金石.高性能計算是未來高精度數(shù)值仿真和數(shù)值試驗的基本要求，其基礎是并行計算機硬件與其上運行的高性能計算軟件，兩者缺一不可.隨著天河系列、曙光系列等高性能計算機的發(fā)展，我國在高性能計算的硬件研究和制造方面已經(jīng)達到國際先進水平，但是在軟件方面，仍缺乏自主研發(fā)的高性能計算系統(tǒng).從飛行器結構強度分析領域來看，目前國內(nèi)可解決飛行器綜合優(yōu)化設計問題的系統(tǒng)性軟件平臺還處于空白.另一方面，由于CAE軟件市場長期由國外壟斷，授權費用昂貴，戰(zhàn)略性、高水準的工業(yè)工程數(shù)值模擬軟件發(fā)展受到限制.

1 高性能并行結構分析與優(yōu)化的若干關鍵技術

現(xiàn)代工業(yè)設計的產(chǎn)品復雜性越來越高，對高性能計算的要求也越來越高，例如我國自主研發(fā)的某型飛機，其精細化模型自由度已經(jīng)達到數(shù)億量級.這對高性能計算能力，包括并行計算、超大規(guī)模高效數(shù)值算法等提出更高的要求，涉及到以下幾類關鍵技術.

1.1 滿足高性能計算體系特點的區(qū)域分解并行算法

并行計算方法的發(fā)展一般遵循2個途徑：重構和分而治之.按照算法設計特點的不同，現(xiàn)有的結構并行分析技術也歸為2類：一類稱為算法級的并行有限分析技術，根據(jù)有限元分析過程中不同環(huán)節(jié)的計算特點，通過重構計算序列來實施并行求解；另一類稱為子結構級的并行有限元分析技術，是指將一個大結構分解成多個子結構，然后根據(jù)分而治之的策略進行并行求解.區(qū)域分解法是后者的典型代表，其算法是高度并行的，即計算的主要步驟在各子域內(nèi)獨立進行，同時允許在不同子域中選用不同的數(shù)學模型，以便整體模型更適合于工程物理實際，具有很好的工程適用性.[13]

采用區(qū)域分解算法進行結構計算時，依照大規(guī)模并行計算特點，先按全機大部件（機翼、機身、中央翼、垂尾、平尾等）分別建立結構的大部件分析模型.每個大部件分別進行自動區(qū)域劃分，按單元的不同計算量分別考慮不同的權數(shù)，形成靜態(tài)負載平衡的若干子結構（與高性能計算機的計算節(jié)點數(shù)匹配）.在每個子結構內(nèi)部將子結構全部自由度縮聚到子結構的邊界點上，形成以子結構邊界點為未知量的子結構剛度矩陣和載荷向量矩陣，按自由度形成全局界面方程.用分布式預處理共軛梯度（Distributed Preconditioned Conjugate Gradient， DPCG）法求解界面方程，為提高并行效率，按結構的受載工況數(shù)量將界面方程劃分為若干個單工況的界面方程進行并行求解.最后，將邊界點的位移回代到每個子結構計算子結構內(nèi)部自由度.[46]求解流程見圖1.

1.2 超大規(guī)模結構變量敏度高效求解技術

經(jīng)過幾十年的發(fā)展，參數(shù)優(yōu)化依然是航空結構優(yōu)化設計的最重要手段.隨著飛機綜合性能和設計生產(chǎn)自動化水平越來越高，超大規(guī)模變量優(yōu)化成為飛機設計的關鍵技術.[3]采用自然網(wǎng)格進行全機有限元模擬時，模型節(jié)點數(shù)一般為數(shù)千至上萬，而采用精細網(wǎng)格進行全機有限元模擬時，節(jié)點數(shù)可達數(shù)十萬至上千萬，精細化模型自然帶來結構設計變量的大幅增長，當前國際先進優(yōu)化設計水平的設計變量為10萬個.[78]

敏度求解的常用方法有3種：有限差分法、半解析法和解析法.按照約束性質(zhì)的不同，超大規(guī)模（10萬以上）結構變量的敏度求解方法分為2類：對數(shù)量較少的全局約束（位移、模態(tài)）采用直接解析法進行求解，對數(shù)量眾多的局部約束（應力、應變）采用伴隨矩陣法等半解析法進行求解.敏度求解過程中需要對剛度矩陣求逆，流程見圖2.

但是，實際工程結構剛度矩陣往往規(guī)模太大或奇異性較強，不能直接求逆，虛載荷法可有效解決該問題.在虛載荷計算中提前建立約束及實位移、虛位移關系表，將實、虛位移挑選過程從敏度計算循環(huán)過程中分離，按邊界條件計算虛載荷，可有效提高敏度計算效率.同時，在結構微分剛度矩陣組裝中建立設計變量與關聯(lián)的屬性卡及單元編號關系表，避免微分剛度計算時反復查找單元和屬性.

1.3 結構非線性振動特性求解技術

在變形較大的情況下，大展弦比機翼位移與應變關系的線性假設不再適用，此時小應變假設不能正確反映其物理本質(zhì)，必須采用大應變理論.在求解結構發(fā)生大變形時的切線剛度矩陣時，基于更新的拉格朗日列式的非線性方程的求解方法NewtonRaphson法，由Gaussia消去法求解增量位移的線性方程組，最終誤差向量由Jacobian法形成結構的切線剛度矩陣（幾何大變形的貢獻）.在得到結構大變形下的切線剛度矩陣后，根據(jù)求解結構的不同特點，分別采用逆冪法、子空間迭代和FEAST法進行非線性振動特性的計算，求解流程見圖3.分別計算載荷包線內(nèi)所有工況下的變形及其非線性振動特性，形成非線性振動特性庫，供耦合分析調(diào)用.[910]

2 高性能并行結構分析與優(yōu)化應用

2.1 結構分析優(yōu)化平臺簡介

航空結構強度分析與優(yōu)化系統(tǒng)HAJIF是中航工業(yè)強度所研制推出的目前國內(nèi)航空領域功能最為全面的大型CAE軟件系統(tǒng).該系統(tǒng)以強度試驗數(shù)據(jù)庫為支撐，提供飛行器結構基礎分析、優(yōu)化設計、氣動彈性分析、熱分析等功能[11]，見圖4.系統(tǒng)提供圖形前后置功能，集成自動快速建模、切面剛度計算、細節(jié)強度分析等航空特色模塊.系統(tǒng)采用先進的開放式、可擴充的軟件架構，可方便地為用戶提供個性化定制開發(fā)和服務.

基礎分析模塊包括線性靜力分析、模態(tài)分析和線性屈曲分析，支持33種結構單元和非結構單元，可滿足航天航空結構常用分析的需求，可解決1 000萬自由度結構的靜力分析、動力固有特性分析、瞬態(tài)動力分析、屈曲模態(tài)分析及屈曲臨界載荷分析等問題.優(yōu)化設計模塊包括準則法對結構進行“滿應力/滿應變”優(yōu)化設計，針對廣義位移、自振頻率、顫振速度、翼面效率、翼面發(fā)散速度等變量進行計算的大規(guī)模變量敏度求解，對結構進行考慮“靜、動、顫、彈”的多約束優(yōu)化設計的規(guī)劃法，對機身、機翼、尾翼等典型結構的結構布局優(yōu)化設計等.系統(tǒng)支持大規(guī)模并行結構分析與優(yōu)化設計，已在中國飛機強度研究所10萬億次計算中心和中國航空研究院100 TFlops高性能計算中心得到驗證.

2.2 基于最大航程的氣動結構綜合優(yōu)化設計

突破滿足真實飛機設計需求、以飛機總體性能為目標的優(yōu)化設計，基于氣動結構多學科優(yōu)化的綜合優(yōu)化技術，建立相應的工作流程，采用基于工程綜合優(yōu)化的方法以及系統(tǒng)級和學科級的兩級優(yōu)化技術實現(xiàn)綜合優(yōu)化.為實現(xiàn)10萬個設計變量的綜合優(yōu)化問題，解決復雜流程中各學科模型間的數(shù)據(jù)交換和管理，突破彈性飛機的氣動力與結構耦合的型架外形計算方法，采用并行度較高的人工智能算法.針對復雜系統(tǒng)的多學科優(yōu)化流程快速建模和工作流引擎的生成問題，采用基于多模型分層次驅(qū)動的體系架構實現(xiàn)計算.

對于亞聲速、高亞聲速氣動力設計，為獲得最大的升阻比，機翼的氣動力分布應當滿足橢圓分布，見圖5中的虛線.但是，對于飛機設計而言，優(yōu)化的終極目標并不是最大升阻比，而是能夠反映飛機整體性能的目標函數(shù).

式（1）表示飛機設計中阻力與飛機結構質(zhì)量之間的取舍關系.如果考慮氣動、結構的綜合設計，在飛機總質(zhì)量保持一定的前提下，增加巡航航程，氣動力分布應該變成圖5中的實線所示.該分布能夠降低機翼根部彎矩，從而減輕結構質(zhì)量.在巡航狀態(tài)下，增加升阻比、減輕結構質(zhì)量可以最大化巡航航程.基于最大航程的氣動結構綜合優(yōu)化設計流程見圖6.

2.3 復合材料機翼精細化模型綜合優(yōu)化設計

全復合材料常規(guī)布局的翼盒，沿展向布置若干根肋，前后布置2根大梁，有限元分析模型共有8萬個單元.優(yōu)化模型中，設計變量約有11萬個（復合材料分層厚度、等效金屬桿元面積等），約束約有150萬個（強度、剛度、穩(wěn)定性、長桁與蒙皮面積比、蒙皮厚度連續(xù)性等），見圖7.在此模型上建立敏度計算與優(yōu)化設計模型，敏度計算時設計變量為復合材料蒙皮的分層厚度，按非均衡設計，每個屬性區(qū)4個變量，從機翼中部到翼稍取8個節(jié)點的垂直弦平面方向位移作為敏度計算約束.從預處理到敏度計算結束共用時74 545 s，其中敏度計算部分74 350 s，占總用時的99.7%，采用并行求解與單核相比計算效率大幅提升.

3 結束語

引入?yún)^(qū)域分解法進行結構并行求解和軟件的并行化改造，可大大提高結構求解規(guī)模和求解效率，為精細化結構分析奠定很好的基礎.通過引入直接解析法敏度求解、虛位移求解、數(shù)據(jù)結構組織等技術，實現(xiàn)超過10萬設計變量敏度的高效求解和超大規(guī)模變量結構優(yōu)化設計.

參考文獻：

[1] HODGSON D C， JIMACK P K. A domain decomposition preconditioner for a parallel finite element solver on distributed unstructured grids[J]. Parallel Computing， 1997， 23（8）： 11571181. DOI： 10.1016/S01678191（97）000471.

[2] GIBBS W R. A parallel/recursive algorithm[J]. Journal of Computational Physics， 2004， 201（2）： 573585. DOI： 10.1016/j.jcp.2004.06.008.

[3] ZUCCHINI A. A parallel preconditioned conjugate gradient solution method for finite element problems with coarsefine mesh formulation[J]. Computers and Structures， 2000， 78（6）： 781787. DOI： 10.1016/S00457949（00）000535.

[4] BASERMANN A， REICHEL B. Preconditioned CG methods for sparse matrices on massively parallel machines[J]. Parallel Computing， 1997， 23（3）： 381398. DOI： 10.1016/S01678191（97）000057.

[5] 苗新強， 金先龍. 結構動力數(shù)值仿真兩級并行計算系統(tǒng)開發(fā)及應用[J]. 計算機輔助設計與圖形學學報， 2015， 27（6）： 1118.

MIAO X Q， JIN X L. Development and application of twolevel parallel computing system for structural dynamic numerical simulation[J]. Journal of ComputerAided Design & Computer Graphics， 2015， 27（6）： 1118.

[6] BASERMANNA A， STEFFENB B. New preconditioned solvers for large sparse eigenvalue problems on massively parallel computers[J]. Advances in Parallel Computing， 1998（12）： 565572. DOI： 10.1016/S09275452（98）800725.

[7] QIU Q Y， LI B. Decomposition method of complex optimization model based on global sensitivity analysis[J]. Chinese Journal of Mechanical Engineering， 2014， 27（4）： 722728. DOI： 10.3901/CJME.2014.0516.096.

[8] 菅帥. 求解對稱特征值問題的修正Jacobi共軛預處理梯度法[J]. 應用數(shù)學與計算數(shù)學學報， 2013， 27（2）： 261287. DOI： 10.3969/j.issn.10066330.2013.02.005.

JIAN S. Modified Jacobiconjugate preconditioned gradient method for symmetric eigenvalue problems[J]. Communication on Applied Mathematics and Computation， 2013， 27（2）： 261287. DOI： 10.3969/j.issn.10066330.2013.02.005.

[9] BERTOLAZZI E， FREGO M. Preconditioning complex symmetric linear systems[J]. Mathematical Problems in Engineering， 2014， 5（1）： 1120.

[10] 姚成寶， 郭永輝， 張柏華， 等. 基于區(qū)域分解的并行有限元程序及其在水下爆炸中的應用[J]. 計算力學學報， 2013， 30（S1）： 150154. DOI： 10.7511/jslx2013z032.

YAO C B， GUO Y H， ZHANG B H， et al. Parallel finite element code based on domain decomposition and its application in underwater explosion[J]. Chinese Journal of Computational Mechanics， 2013， 30（S1）： 150154. DOI： 10.7511/jslx2013z032.

[11] 孫俠生， 段世慧， 陳煥星. 堅持自主創(chuàng)新， 實現(xiàn)航空CAE軟件的產(chǎn)業(yè)化發(fā)展[J]. 計算機輔助工程， 2010， 19（1）： 2126. DOI： 10.3969/j.issn.10060871.2010.01.003.

SUN X S， DUAN S H， CHEN H X. Keeping independent innovation， implementing industrialization development of aviation CAE software[J]. Computer Aided Engineering， 2010， 19（1）： 2126. DOI： 10.3969/j.issn.10060871.2010.01.003.

（編輯 武曉英）