航空結(jié)構(gòu)大規(guī)模并行分析與優(yōu)化應(yīng)用

常亮　段世慧　王立凱　羅利龍

摘要：

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

關(guān)鍵詞：

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

中圖分類號(hào)： V221

文獻(xiàn)標(biāo)志碼： B

0 引 言

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

1 高性能并行結(jié)構(gòu)分析與優(yōu)化的若干關(guān)鍵技術(shù)

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

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

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

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

1.2 超大規(guī)模結(jié)構(gòu)變量敏度高效求解技術(shù)

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

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

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

1.3 結(jié)構(gòu)非線性振動(dòng)特性求解技術(shù)

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

2 高性能并行結(jié)構(gòu)分析與優(yōu)化應(yīng)用

2.1 結(jié)構(gòu)分析優(yōu)化平臺(tái)簡(jiǎn)介

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

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

2.2 基于最大航程的氣動(dòng)結(jié)構(gòu)綜合優(yōu)化設(shè)計(jì)

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

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

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

2.3 復(fù)合材料機(jī)翼精細(xì)化模型綜合優(yōu)化設(shè)計(jì)

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

3 結(jié)束語(yǔ)

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

參考文獻(xiàn)：

[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] 苗新強(qiáng)， 金先龍. 結(jié)構(gòu)動(dòng)力數(shù)值仿真兩級(jí)并行計(jì)算系統(tǒng)開發(fā)及應(yīng)用[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 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] 菅帥. 求解對(duì)稱特征值問題的修正Jacobi共軛預(yù)處理梯度法[J]. 應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué)學(xué)報(bào)， 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ū)域分解的并行有限元程序及其在水下爆炸中的應(yīng)用[J]. 計(jì)算力學(xué)學(xué)報(bào)， 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] 孫俠生， 段世慧， 陳煥星. 堅(jiān)持自主創(chuàng)新， 實(shí)現(xiàn)航空CAE軟件的產(chǎn)業(yè)化發(fā)展[J]. 計(jì)算機(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.

（編輯 武曉英）