范志瑞， 閻 軍*， 牛 斌， 隋倩倩， 許 琦， 蔣存存， 趙國忠

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，工程力學(xué)系，大連 116024；2.大連理工大學(xué) 機械工程學(xué)院，大連 116024)

1 引 言

裝備結(jié)構(gòu)變形后整體或局部的幾何形狀往往對其結(jié)構(gòu)性能具有重要的影響，尤其對于考慮幾何大變形效應(yīng)的結(jié)構(gòu)設(shè)計問題。在此類問題中，與傳統(tǒng)基于小變形假設(shè)的結(jié)構(gòu)設(shè)計問題不同，變形前后的構(gòu)型存在較大的差別[1,2]，這使得裝備結(jié)構(gòu)面臨設(shè)計的構(gòu)型和實際在位運行的構(gòu)型不一致的問題。需要說明的是，一般對于承載結(jié)構(gòu)，可近似認為初始構(gòu)型面向設(shè)計和制造，而變形后的構(gòu)型則對應(yīng)結(jié)構(gòu)實際在位運行情況，即在考慮實際載荷(如環(huán)境載荷、安裝載荷和工作載荷等) 作用后，結(jié)構(gòu)會變形為何種構(gòu)型。在實際工程中，變形后的構(gòu)型往往融入了設(shè)計者對特殊性能/功能的需求，如翼面負載和流體阻力作用后機翼表面仍要求保持良好的空氣動力學(xué)性能[3-5]，展開后衛(wèi)星太陽能帆板能夠確保充足的太陽能吸收[6,7]，太空天線的展開形狀能夠使得信號接受質(zhì)量最佳[8,9]，以及變形后的軟體機械手對抓取物體表面的貼合能力[10,11]等。因此，變形后構(gòu)型的精確分析和控制對結(jié)構(gòu)性能及功能的實現(xiàn)具有重要意義。

當采用傳統(tǒng)有限元方法處理幾何大變形問題時，通常以變形前的構(gòu)型為參考，求解變形后的構(gòu)型[1,12],此類分析方法即為正向運動分析方法。由于是以變形前的初始構(gòu)型為結(jié)構(gòu)分析的出發(fā)點，該方法很難適用于上述變形后構(gòu)型的設(shè)計問題。而對于變形后構(gòu)型的設(shè)計問題，更合理的方式是以變形后的構(gòu)型為參考，逆向求解變形前的構(gòu)型，即采用逆向運動分析方法。Schield[13]利用坐標變換的方法以及超彈性方程的對偶性建立了最早的逆向運動分析理論。Chadwick[14]將Eshelby能量-動量張量[15]引入了逆向運動分析理論，建立了逆向運動分析的虛功平衡方程。然而，Eshelby能量-動量張量仍然是以初始構(gòu)型為參考而定義的，其求解過程較為復(fù)雜。為此，Govindjee等[16,17]采用了更為便捷的Cauchy應(yīng)力來對逆向運動問題進行求解。Wallin等[18]對逆向運動問題進行了系統(tǒng)推導(dǎo)，并建立了該類分析問題的有限元求解方法。

在考慮幾何大變形/有限應(yīng)變的有限元分析中，材料的剛度張量、形函數(shù)、內(nèi)力載荷向量以及切線剛度矩陣等均與結(jié)構(gòu)位移場相關(guān)[1]，因此往往需要采用非線性分析的迭代法求解，并且上述物理量在各個迭代步中不斷更新。當結(jié)構(gòu)變形較大時，為了保證分析迭代的收斂性，需要大量的載荷增量步，從而使得非線性迭代的計算成本急劇增加。雖然修正的Newton-Raphson迭代法[19,20]可以采用固定的切線剛度矩陣對非線性問題進行求解，但卻使得求解過程喪失了二階收斂速度，有時甚至導(dǎo)致迭代過程無法收斂。

近年來，隨著計算技術(shù)的不斷發(fā)展，并行計算在結(jié)構(gòu)有限元分析中發(fā)揮著越來越重要的作用[21,22]。根據(jù)內(nèi)存架構(gòu)，并行計算可分為共享式內(nèi)存[23]、分布式內(nèi)存[24,25]和混合式內(nèi)存[26]三種。根據(jù)并行粒度，并行計算通常包含線程并行和進程并行兩種。一般而言，共享式內(nèi)存和線程模式適合單機并行，具有較快的內(nèi)存訪問速度，并且計算資源開銷較小。而采用分布式內(nèi)存和進程模式的并行方案，則在并行的可擴展性及大規(guī)模集群計算方面更具優(yōu)勢。MPI(Massage Passing Interface)[24,25]是基于分布式內(nèi)存和進程模式開發(fā)的標準并行協(xié)議，在并行計算中得到了廣泛應(yīng)用。但是，其涉及復(fù)雜的消息通信過程以及并行協(xié)同機制，在一定程度上限制了其在數(shù)值計算方面的應(yīng)用。為此，Argonne國家實驗室開發(fā)了可移植可擴展科學(xué)計算庫PETSc(Portable,Extensible Toolkit for Scientific Computation)[27]，其對消息通信過程進行了封裝，并且提供了高效的線性方程組求解算法。值得一提的是，目前學(xué)者已經(jīng)基于PETSc解決了涉及十億級[28]，甚至百億級[29]變量的結(jié)構(gòu)優(yōu)化設(shè)計問題。

2 逆向運動分析理論概述

2.1 結(jié)構(gòu)變形描述

以下對逆向分析的理論進行概述，關(guān)于正向運動的非線性分析理論可參考文獻[1]。在結(jié)構(gòu)逆向運動分析中，變形描述與正向運動分析相反。如 圖1 所示，在t時刻，結(jié)構(gòu)體變形后的構(gòu)型Ω已知，并且該構(gòu)型由某一未知的初始構(gòu)型Ω0通過外載荷作用得到。變形前后構(gòu)型的邊界分別記作?Ω0和?Ω。為研究基于逆向運動的結(jié)構(gòu)分析方法，在此定義由變形后的構(gòu)型指向初始構(gòu)型的映射函數(shù)φ，即

X=φ(x)

(1)

式中X∈Ω0為初始構(gòu)型內(nèi)部某一物質(zhì)點的位置，x∈Ω為物質(zhì)點在變形后構(gòu)型的位置。

圖1 逆向運動分析的結(jié)構(gòu)變形描述

與正向運動分析類似，逆向運動分析的變形梯度張量可定義為

f=φ

(2)

f=F-1

(3)

式中F為正向運動的變形梯度。為了保證映射φ的唯一性，要求J=det(f)>0。此外，變形前后的應(yīng)變能密度ω0和ω之間滿足

ω=Jω0

(4)

2.2 虛功方程及其有限元離散化

由于假設(shè)了變形前后構(gòu)型之間滿足微分同胚，故正向運動的平衡方程同樣適用于逆向運動分析過程。當忽略結(jié)構(gòu)的體力且以變形后構(gòu)型為參考時，虛功方程可表述為

(5)

σ=ωI-fT?ω/?f

(6)

(7)

式中I為二階單位張量。

為了能夠采用有限元方法對式(5)進行求解，對其進行Taylor展開，可得

(8)

采用Galerkin方法對上式進行離散化，則有

(9)

(10)

(11)

式中H為變形梯度的向量格式f與節(jié)點位移之間的變換矩陣，即df=Hda，C為逆向運動的剛度張量的Voigt格式，且定義為

=?σ/?f

(12)

根據(jù)1s tPiola-Kirchhoff應(yīng)力張量P之間的關(guān)系σ=JPFT，聯(lián)合式(12)可求得的表達式為

(13)

(14)

根據(jù)式(9)，定義Newton-Raphson迭代過程的不平衡載荷向量為

R=Fint-Fext

(15)

式中Fint和Fext分別為內(nèi)外載荷向量，即

(16，17)

將式(10)和由式(9)推導(dǎo)出的式(15)代入式(8)可得，逆向運動的位移迭代可表示為

Kda(i )=-R，a(i + 1)=a(i )+da(i )

(18，19)

式中下標i表示第i次迭代。

2.3 Neo -Hookean材料分析

本文采用Neo -Hookean超彈性材料模型[30,31]，該模型對于大應(yīng)變情況具有很好的適用性。在逆向分析中，Neo -Hookean 材料模型的應(yīng)變能密度ω的數(shù)學(xué)表達式[18]為

(20)

式中K和G分別為材料的體模量和剪切模量，c=fTf為逆向運動的右Cauchy-Green應(yīng)變張量，tr()表示張量的跡。

將式(7，20)代入式(12)，可得Cauchy應(yīng)力的表達式為

(21)

聯(lián)合式(3,4,14,20)可得Updated-Lagrange迭代法的剛度張量為

(22)

3 基于PETSc的逆向運動并行分析框架

為了提高結(jié)構(gòu)非線性分析的效率，本研究基于PETSc(Portable,Extensible Toolkit for Scientific Computation) 并行函數(shù)庫[27]建立了逆向運動的并行分析程序框架。如圖2所示，PETSc是以Blas/Lapack[32,33]為底層內(nèi)核建立的并行計算庫。同時，信息傳遞采用標準MPI協(xié)議[24]，保證了庫函數(shù)良好的跨平臺性以及對大規(guī)模集群計算的支持。在基礎(chǔ)函數(shù)庫的基礎(chǔ)上，PETSc采用面向?qū)ο笏枷虢⒘讼蛄?Vector) 、矩陣(Matrix)、索引集IS(Index Set)、映射(Mapping)、KSP(Krylov Subspace Package)線性求解器及預(yù)處理算子PC(Preconditioner) 等數(shù)據(jù)結(jié)構(gòu)。通過對數(shù)據(jù)結(jié)構(gòu)和方法(數(shù)學(xué)運算及消息傳遞等)的封裝，為用戶提供了簡便易行的并行解決方案。用戶無需關(guān)心具體數(shù)據(jù)結(jié)構(gòu)及進程之間的消息傳遞模式，而將更多的精力放在頂層并行框架設(shè)計層面。因此，在保證程序穩(wěn)定性的同時，為開發(fā)人員提供了靈活的程序設(shè)計空間，提高了開發(fā)的效率。

圖2 PETSc自下而上的數(shù)據(jù)結(jié)構(gòu)[27]

3.1 數(shù)據(jù)并行存儲與交互

在MPI并行計算中，數(shù)據(jù)會分割成許多片段分布式存儲在不同進程的內(nèi)存空間中。而各個進程的內(nèi)存是獨立的，不能相互直接訪問。進程之間的數(shù)據(jù)交換須在遵循標準MPI協(xié)議的情況下，采用進程通信的方式進行。

在數(shù)據(jù)交互前，需要對每個進程上涉及的單元計算任務(wù)進程合理劃分。劃分的規(guī)則為

ni=nb+nz,i，nb=div(ne/m)

(23，24)

(25)

式中ni為當前進程i的單元計算量，nb和nz,i分別為各個進程基本計算量以及當前進程額外所需的計算量,ne為單元總數(shù)，m為MPI進程數(shù),div()和mod()分別為整數(shù)除法的求商和取余運算。位移向量、內(nèi)外力載荷向量以及標識邊界條件的向量在各個進程上的分布也可以類似采用上述劃分規(guī)則。對于切線剛度矩陣，可遵循上述規(guī)則按行劃分為m個子矩陣，分別存儲于相應(yīng)進程的內(nèi)存中。

對于單元節(jié)點編號和節(jié)點坐標，以及應(yīng)力和應(yīng)變場的存儲需匹配單元計算的劃分方式。以節(jié)點坐標為例，最終的單元坐標信息在各個進程上的劃分和存儲方式如下。

p=[p(1),p(2),…,p(i )，…,p(m )]

(26)

(27)

(28)

此外，最終節(jié)點坐標向量的生成需要借助進程通信。如圖3所示，首先從本地并行讀取已生成的節(jié)點坐標二進制信息，形成初始節(jié)點坐標向量。其次，根據(jù)單元的節(jié)點編號信息生成索引集(IS)，用于確定當前單元節(jié)點坐標信息在初始節(jié)點坐標向量的位置。最后，根據(jù)索引集生成分散器(Sca-tter)，完成進程之間坐標信息的通信和傳遞。

圖3 基于PETSc庫的并行框架的數(shù)據(jù)并行存儲與交互

完成節(jié)點信息交互后，不同進程會同時擁有共有節(jié)點信息的副本。因此，雖然基于MPI的并行計算能夠提高計算效率，但卻增加了額外的內(nèi)存開銷。

3.2 通信鎖死與負載平衡

在并行計算中，需要注意進程之間的負載平衡問題，即要保證每個進程的計算任務(wù)量大致相等。如果某些進程過早地進入下一個計算環(huán)節(jié)，很容易造成進程通信阻塞，進而引起程序卡死現(xiàn)象。以 圖4 所示的計算過程為例，進程1和進程2均參與了單元分析及整體剛度矩陣組集過程。由于兩個進程的任務(wù)量分配不均勻，進程1過早地進入了后續(xù)計算處理環(huán)節(jié)(如邊界條件處理)。如果后續(xù)計算環(huán)節(jié)為全局操作，進程1需要等待進程2的通信響應(yīng)。然而，進程2仍然處于之前的計算環(huán)節(jié)，且該環(huán)節(jié)中計算所得的單元剛度矩陣需要通過進程通信將其組裝至全局剛度矩陣。因此，進程2也需要進程1的通信響應(yīng)后才能完成相應(yīng)計算。由于進程1和進程2都處于進程通信的等待狀態(tài)，此時會造成進程通信鎖死，從而使得計算無法繼續(xù)。

圖4 由于負載不平衡而引起的進程通信鎖死

針對上述問題，可選的解決方案之一是在計算任務(wù)量較少的進程上再執(zhí)行一次空的單元分析，即用各項均為零的單元剛度矩陣執(zhí)行一次全局剛度矩陣的組集任務(wù)。具體地，如果mod(ne/np)≠0，且nz=0時，需要在該進程引入空的單元分析。當然，在實際計算中應(yīng)盡量保證每個進程的計算量一致，其有助于提高分析的效率，因為即使空的單元分析仍然需要花費一定的進程通信成本。

3.3 基于PETSc的逆向運動并行分析框架

圖5給出了基于PETSc并行庫所建立的非線性逆向運動并行分析程序框架。為了提高結(jié)構(gòu)分析的效率，主體程序采用以C++為基礎(chǔ)編程語言的PETSc函數(shù)庫編寫。同時，基于Python語言對ABAQUS商業(yè)有限元軟件[34]進行二次開發(fā)，搭建了完善的前后處理模塊。并行程序的分析流程如下。

(1) 基于ABAQUS軟件生成有限元計算所需的模型信息，并以二進制方式寫入本地磁盤，然后采用標準MPI文件接口讀入主分析程序。

(2) 確定各個進程上單元分析的計算量，并進一步生成用于MPI進程通信的索引集和分散器。

(3) 根據(jù)讀入的有限元模型信息，調(diào)用結(jié)構(gòu)分析模塊對當前Newton-Raphson載荷增量步進行迭代求解。此過程需調(diào)用單元分析和材料分析模塊。

(4)判斷是否滿足當前迭代步的收斂條件。如果不滿足，則繼續(xù)執(zhí)行步驟(3)；否則開始新的載荷增量步。

(5) 當所有載荷增量步迭代均收斂后，將計算結(jié)果以二進制方式并行寫入本地磁盤?；贏BAQUS軟件的Python接口，創(chuàng)建ODB(Output DataBase)文件庫，對計算結(jié)果進行可視化顯示。

圖5 基于PETSc的非線性逆向運動并行分析程序框架

4 數(shù)值算例

算例1逆向運動分析的驗證

如圖6所示長方形板結(jié)構(gòu)，其長度L=200 mm，寬度H=50.0 mm，厚度t=1.0 mm。板的左邊界固定，右端施加豎直向上和大小為8 N/mm的分布載荷F。體模量K=1.45 GPa，剪切模量G=0.81 GPa。在有限元分析中，結(jié)構(gòu)離散化為10000個平面單元，載荷增量步的數(shù)目設(shè)置為200，MPI進程數(shù)為4。計算測試平臺為HP-Z840工作站，操作系統(tǒng)為Ubuntu 18.04 LTS，CPU型號為 Intel?Xeon(R) CPU E5-2697 v4@2.30 GHz×72，內(nèi)存大小為256 GB，內(nèi)存條型號為32 GB DDR4-2133 ECC(LR) RAM，硬盤型號為HP Z Turbo Drive 256 GB SSD?；谡蜻\動和逆向運動分析方法分別對上述結(jié)構(gòu)進行分析。圖6給出了構(gòu)型分別作為正向運動的初始構(gòu)型以及逆向運動的變形后構(gòu)型。兩種方法的分析結(jié)果如圖7所示。

圖6 逆向運動分析方法的驗證算例結(jié)構(gòu)

圖7 基于正向運動和逆向運動方法計算所得圖6結(jié)構(gòu)位移場

正向運動分析的結(jié)果如圖7(a)所示。由于正向運動是用于求解外力驅(qū)動下結(jié)構(gòu)變形后的狀態(tài)，因此結(jié)構(gòu)變形與載荷施加的方向一致。但是，逆向運動則與此相反，其用于求解變形前的狀態(tài)。因此，如圖7(b)所示，結(jié)構(gòu)變形的方向與載荷施加的方向相反。如之前所述，由于在逆向運動分析中結(jié)構(gòu)的切線剛度矩陣為非對稱矩陣，并且其本構(gòu)分析也更為復(fù)雜。因此，逆向運動分析需要耗費更多的CPU計算時間。在此算例中，對于單次Newton-Raphson迭代，正向運動分析所花費的平均CPU時間為0.697 s，而逆向運動分析則需要0.954 s。相比于正向運動分析，逆向運動分析的CPU時間增加了36.881%。此外，如果假設(shè)初始構(gòu)型未知，并且需要通過正向運動分析方法使得結(jié)構(gòu)變形后的構(gòu)型恰好為圖6所示構(gòu)型，此分析問題對于正向運動分析方法而言是十分困難的。因為在正向運動分析中，結(jié)構(gòu)信息均是以初始構(gòu)型為參考的。初始構(gòu)型的未知使得結(jié)構(gòu)信息的定義和求解喪失了先決條件。

為了驗證逆向運動分析的精度，以圖7(b)所示逆向運動分析求解所得初始構(gòu)型作為分析的起點，采用正向運動分析方法求解該構(gòu)型所對應(yīng)的變形后構(gòu)型。最終，求得的變形后構(gòu)型與圖6構(gòu)型的偏差如圖8所示?？芍畹淖畲笾禐?3.934×10-2mm。與最大位移47.150 mm相比，該偏差可忽略不計。因此，本研究發(fā)展的逆向運動并行分析框架和程序可以實現(xiàn)對變形前構(gòu)型的精確求解。

圖8 逆向運動分析所得初始構(gòu)型與指定初始構(gòu)型間的偏差

算例2并行計算效率分析

對逆向運動分析的并行效率進行研究。采用的結(jié)構(gòu)如圖9所示，方形板的長度L=400 mm，寬度H=250 mm，厚度t=1 mm。結(jié)構(gòu)左邊界固定，右邊界中點處施加豎直向上和大小為 5 N 的集中載荷。材料屬性取值與算例1相同。采用逆向運動方法對結(jié)構(gòu)進行有限元分析。分析中，結(jié)構(gòu)離散為ne=100000個單元，外載荷均勻地劃分為4個載荷增量步。選用不同的MPI進程數(shù)對并行效率進行測試，同時為了保證每個進程的負載平衡，MPI進程數(shù)分別設(shè)置為m={1,2,4,5,8,10,16,20,32}。與算例1相同，仍然采用單次Newton-Raphson迭代的CPU時間來衡量并行效率。

圖9 非線性并行程序效率測試結(jié)構(gòu)

(29)

式中T0為串行計算所需的平均CPU時間，T(n)為當MPI進程數(shù)為n時計算所需的平均CPU時間。

圖10 單元數(shù)目ne=100000時，并行效率隨MPI進程數(shù)的變化

圖11 單元數(shù)目ne=1600000時，并行效率隨MPI進程數(shù)的變化

5 結(jié) 論

本研究針對考慮大變形/有限應(yīng)變情況下的變形后構(gòu)型控制問題，引入了逆向運動分析方法，同時為了提高非線性有限元分析過程的效率，基于PETSc庫建立了并行計算框架。對并行框架的模塊劃分、數(shù)據(jù)并行存儲與交互以及通信鎖死與負載平衡等關(guān)鍵問題作了詳解闡述。通過算例分析，得到結(jié)論如下。

(1) 逆向運動分析方法可以實現(xiàn)對具有大變形特征的結(jié)構(gòu)變形前構(gòu)型的精確求解，其誤差與結(jié)構(gòu)的最大位移相比很小。

(2) 由于剛度張量的求解復(fù)雜性，以及切線剛度矩陣的非對稱性，逆向運動分析相比于正向運動分析需要花費更多的計算成本。

(3) 采用并行分析框架可大幅度提高逆向運動分析的效率。在算例2中，與單MPI進程相比，采用并行框架可節(jié)約2/3的計算時間，證明了本研究建立并行框架的高效性。

(4) 合理選擇MPI進程數(shù)目會使得并行效率大幅度提升。相反，選擇不恰當?shù)腗PI進程數(shù)目，則會使得進程通信成本增加，計算效率反而下降。