張玉輝，姜清輝

(武漢大學(xué)土木建筑工程學(xué)院，武漢 430072)

傳統(tǒng)彈性力學(xué)法求解邊值問(wèn)題建立在聯(lián)立平衡微分方程、幾何方程以及物理方程的基礎(chǔ)上。其中平衡微分方程和幾何方程由靜力學(xué)平衡和運(yùn)動(dòng)學(xué)導(dǎo)出，無(wú)不確定性。而物理方程是在對(duì)觀測(cè)數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)總結(jié)和擬合的基礎(chǔ)上建立的，包括根據(jù)實(shí)驗(yàn)數(shù)據(jù)識(shí)別現(xiàn)有模型參數(shù)，或者擬合提出一種新的材料物理模型，至今依舊是力學(xué)研究的熱點(diǎn)[1-3]，具有經(jīng)驗(yàn)性和不確定性。

距離最小化數(shù)據(jù)驅(qū)動(dòng)計(jì)算方法是在大數(shù)據(jù)時(shí)代背景下由Kirchdoerfer 和Ortiz[4]首先提出，通過(guò)在應(yīng)力-應(yīng)變相空間內(nèi)仿照應(yīng)變能密度和應(yīng)變余能密度的形式定義空間距離，在材料的宏觀力學(xué)響應(yīng)數(shù)據(jù)庫(kù)內(nèi)為每個(gè)積分點(diǎn)搜索并匹配最可能的應(yīng)力-應(yīng)變數(shù)據(jù)點(diǎn)。該方法在計(jì)算過(guò)程中不再依賴顯式表達(dá)式形式的材料本構(gòu)關(guān)系，將聯(lián)立彈性力學(xué)三大方程的顯式求解問(wèn)題轉(zhuǎn)變?yōu)樗阉髌ヅ渥羁赡艿撵o力許可應(yīng)力場(chǎng)、變形許可應(yīng)變場(chǎng)和數(shù)據(jù)分配方案的最優(yōu)化問(wèn)題，最大程度上保留了觀測(cè)數(shù)據(jù)的原始性。各個(gè)搜索對(duì)象分別對(duì)應(yīng)平衡微分方程約束下的可行解集、幾何方程約束下的可行解集和物理本構(gòu)關(guān)系約束下的可行解子集。

針對(duì)數(shù)據(jù)點(diǎn)的質(zhì)量，Kirchdoerfer 和Ortiz[5]在原有框架的基礎(chǔ)上，提出了熵最小化數(shù)據(jù)驅(qū)動(dòng)計(jì)算方法，降低了數(shù)據(jù)噪音對(duì)計(jì)算精度的影響；針對(duì)邊值問(wèn)題的類型，該方法的適用范圍被陸續(xù)推廣到動(dòng)力學(xué)問(wèn)題[6]、非線性彈性問(wèn)題[7]以及非彈性問(wèn)題[8]。在理論上，Conti 等[9]論證了該方法在數(shù)學(xué)層面上的一系列問(wèn)題，Nguyen 等[10]給出了該方法在變分框架下的一般形式；在應(yīng)用上，Yang 等[11]和Xu 等[12]將該方法與多尺度有限元法相結(jié)合，拓展了求解復(fù)合材料等效連續(xù)介質(zhì)問(wèn)題的新思路，Stainier 等[13]基于數(shù)字圖像關(guān)聯(lián)技術(shù)提出的數(shù)據(jù)識(shí)別方法能夠生成更有效的材料數(shù)據(jù)庫(kù)；在運(yùn)行效率上，分層搜索方法[14]，k-d 樹(shù)搜索方法[15]的提出有效降低了計(jì)算開(kāi)銷。

對(duì)比傳統(tǒng)彈性力學(xué)法，該計(jì)算方法除了以數(shù)據(jù)庫(kù)替換本構(gòu)方程外，同時(shí)還要求輸入定義空間距離時(shí)使用的平衡參數(shù)。在一維桁架問(wèn)題中，表現(xiàn)為一常數(shù)因子，記為C，用于平衡應(yīng)力和應(yīng)變單位上巨大的量級(jí)差；在平面和空間問(wèn)題中表現(xiàn)為正定的常數(shù)矩陣，以保證距離度量始終為正，記為C，在平衡應(yīng)力應(yīng)變量綱的基礎(chǔ)上為各個(gè)方向上的分量也定義了權(quán)重。已發(fā)表的文獻(xiàn)多是針對(duì)數(shù)據(jù)庫(kù)容量[4-7,11-12]及數(shù)據(jù)點(diǎn)質(zhì)量[5]對(duì)計(jì)算精度和運(yùn)行效率的影響展開(kāi)討論，對(duì)于平衡因子取值影響的討論卻很少，或是針對(duì)特定的模型和邊界條件，僅以計(jì)算結(jié)果的好壞給出取值建議[7,11-12]。

本文分別對(duì)線彈性桁架問(wèn)題和線彈性平面問(wèn)題中，平衡因子取值對(duì)算式的收斂速率和性能的影響進(jìn)行了研究。結(jié)果表明：盡管平衡因子只是數(shù)值而沒(méi)有物理上的意義，不代表任何的材料力學(xué)響應(yīng)，為保證計(jì)算精度，一定程度上仍需要反映材料自身的物理本構(gòu)關(guān)系，這種依賴性將隨著數(shù)據(jù)點(diǎn)密度的上升而降低。

1 距離最小化數(shù)據(jù)驅(qū)動(dòng)計(jì)算方法

上述整形約束二次規(guī)劃問(wèn)題采用的求解策略是交替極小化算法。第一個(gè)極小化問(wèn)題的描述為：從滿足物理本構(gòu)關(guān)系的應(yīng)力-應(yīng)變場(chǎng)出發(fā)，找到相距最近的應(yīng)力可行解和應(yīng)變可行解；第二個(gè)極小化問(wèn)題的描述為，從滿足力系平衡和幾何條件的應(yīng)力-應(yīng)變可行解出發(fā)，找到相距最近的應(yīng)力-應(yīng)變數(shù)據(jù)點(diǎn)。由此形成一個(gè)定點(diǎn)迭代過(guò)程：

式中：Zj表示第j步迭代中的數(shù)據(jù)點(diǎn)分配；PM表示子問(wèn)題一向所有應(yīng)力-應(yīng)變可行解的集合(EL/KL)映射；PC表示子問(wèn)題二向所有數(shù)據(jù)分配方案的集合(CL)映射。如圖1 所示。

圖1 迭代流程Fig. 1 Iterative procedure

算法的目標(biāo)在于，找到應(yīng)力-應(yīng)變相空間內(nèi)，圖1 中兩個(gè)集合的交集或使二者距離取得最小值時(shí)的數(shù)據(jù)分配Zult。當(dāng)Zj不再發(fā)生變化時(shí)，認(rèn)為迭代收斂，滿足平衡和幾何約束條件的可行解輸出為力學(xué)邊值問(wèn)題的近似解。輸入不同的迭代初值Z0、距離度量方式以及數(shù)據(jù)庫(kù)時(shí)，輸出的解不一定相同。

解線性方程式(8)和式(10)可得到節(jié)點(diǎn)位移列陣和拉格朗日乘子列陣，代入式(2)和式(9)即可得到應(yīng)變可行解和應(yīng)力可行解。

第二個(gè)極小化問(wèn)題中約束條件式(2)、式(3)已滿足，只剩約束條件式(4)。有限元格式下，可以將問(wèn)題分治到各個(gè)積分點(diǎn)，每個(gè)積分點(diǎn)均匹配最近的材料數(shù)據(jù)點(diǎn)時(shí)，全局罰函數(shù)最小，對(duì)應(yīng)的數(shù)學(xué)模型為：

算法的本質(zhì)是將碎片化后材料力學(xué)響應(yīng)作為整形約束引入，從而找到滿足平衡方程和幾何方程，且與整形約束沖突程度最小的近似解。

2 桁架問(wèn)題中常數(shù)C 的影響

文中使用的桁架和平面數(shù)據(jù)驅(qū)動(dòng)計(jì)算程序均基于Visual Studio 19.0 平臺(tái)C++語(yǔ)言開(kāi)發(fā)。

2.1 松弛問(wèn)題的收斂速率

對(duì)收斂速率展開(kāi)的分析建立在數(shù)據(jù)庫(kù)直接由材料宏觀力學(xué)特性表征的基礎(chǔ)上，即約束條件式(4)中去掉整形約束，數(shù)據(jù)庫(kù)替換為材料本構(gòu)方程曲線：

式中，m=C/E為平衡常數(shù)與彈性模量的比例系數(shù)。

式(34)和式(35)表明，應(yīng)力-應(yīng)變可行解分別以恒定的收斂比線性收斂。以應(yīng)力作軸，并將應(yīng)變乘以彈性模量E投影到軸上，如圖2 所示。

圖2 桁架問(wèn)題迭代過(guò)程Fig. 2 Iterative process of truss problem

當(dāng)?shù)絢趨于無(wú)窮時(shí)，應(yīng)力-應(yīng)變可行解都將收斂于精確解，而參數(shù)m只影響二者收斂的快慢。式(36)表明，分配的數(shù)據(jù)點(diǎn)始終位于該軸上應(yīng)力和應(yīng)變可行解之間，并最終恢復(fù)為有限元解。邊值條件對(duì)迭代路徑的影響主要體現(xiàn)在對(duì)可行解集合的限制，如靜定桁架結(jié)構(gòu)的應(yīng)力可行解始終保持為參考解(存在且只有一個(gè)靜力許可應(yīng)力場(chǎng))，在迭代過(guò)程中僅應(yīng)變可行解發(fā)生改變。顯然，該類問(wèn)題中平衡常數(shù)C的最優(yōu)取值趨近于0，以加速應(yīng)變解收斂。

2.2 原始問(wèn)題的收斂性能

2.1 節(jié)中的討論基于數(shù)據(jù)點(diǎn)連續(xù)無(wú)間斷的假設(shè)，但是在實(shí)際情況中，數(shù)據(jù)點(diǎn)是離散且有限的，迭代過(guò)程不可能一直持續(xù)下去。該節(jié)主要討論數(shù)據(jù)庫(kù)為均勻分布的離散數(shù)據(jù)點(diǎn)時(shí)，平衡常數(shù)C對(duì)算式收斂性能的影響?；謴?fù)整形約束后的數(shù)據(jù)庫(kù)集合為：

模型選用圖3(a)所示的代表性桁架[16]，位移約束和邊界荷載如圖所示，彈性模量E取10 MPa，兩桿橫截面面積取0.01 m2，荷載大小取3 kN，生成數(shù)據(jù)點(diǎn)的平均應(yīng)力間距取10 kPa。平衡常數(shù)C取2×107時(shí)，對(duì)相同的迭代初值，松弛問(wèn)題和原始問(wèn)題對(duì)應(yīng)的迭代路徑如圖3 所示。

對(duì)比圖3(b)及圖3(c)，當(dāng)?shù)綌?shù)較低時(shí)，二者路徑基本相同，但是在迭代末期，圖3(c)所示迭代過(guò)程在遠(yuǎn)離參考解的位置收斂。這是由于原始問(wèn)題中數(shù)據(jù)點(diǎn)作為整形約束，其間距限制了相鄰迭代步間數(shù)據(jù)點(diǎn)所能產(chǎn)生的最小改變量，形成了天然的收斂準(zhǔn)則[17]。當(dāng)自變量的變化不足以克服數(shù)據(jù)點(diǎn)間距時(shí)，算法將過(guò)早收斂。

以2.1 節(jié)松弛問(wèn)題中應(yīng)力-應(yīng)變可行解的收斂速率聯(lián)立式(15)，應(yīng)力數(shù)據(jù)在相鄰兩次迭代中的改變量計(jì)算如下：

2)數(shù)據(jù)點(diǎn)間距 σdis；給定迭代初值以及平衡常數(shù)后，應(yīng)力-應(yīng)變可行解將以恒定的收斂比線性收斂。迭代步數(shù)越多，收斂解的精度越高，這就要求數(shù)據(jù)點(diǎn)間距盡可能小。

3)系數(shù)m；式(40)表明，m=1 時(shí)，數(shù)據(jù)點(diǎn)的線性收斂比為50%；當(dāng)m≠1 時(shí)，數(shù)據(jù)點(diǎn)的漸進(jìn)線性收斂比將與應(yīng)力-應(yīng)變可行解中較大的一方相等(如圖3(b)中，m=2 時(shí)，應(yīng)力可行解在第4 迭代步時(shí)已逼近參考解，此時(shí)數(shù)據(jù)點(diǎn)的收斂比與應(yīng)變可行解一致)。當(dāng)收斂階均為線性時(shí)，收斂比越大，收斂速率越慢，相鄰迭代步間自變量的變化量越小，越容易滿足停步準(zhǔn)則，早熟收斂現(xiàn)象愈發(fā)顯著。

圖3 松弛問(wèn)題和原始問(wèn)題迭代路徑對(duì)比Fig. 3 Comparison of iterative paths between relaxation problem and original problem

從計(jì)算精度的角度看，當(dāng)m＞1 時(shí)，經(jīng)過(guò)充分迭代后，可以忽略應(yīng)力項(xiàng)將式(42)改寫為：

m越大，收斂時(shí)應(yīng)變誤差越大；同理m＜1 時(shí)有：

m越小，收斂時(shí)應(yīng)力誤差越大；

本節(jié)分析表明，數(shù)據(jù)點(diǎn)的離散性將迭代步數(shù)k限制為有限的整數(shù)，這就要求C的取值能夠在有限的迭代步內(nèi)，使應(yīng)力-應(yīng)變可行解均充分收斂至參考解鄰域。當(dāng)m等于1 時(shí)，應(yīng)力-應(yīng)變解均以50%的收斂比線性收斂，二者相對(duì)精度不易產(chǎn)生較大差異。同時(shí)，數(shù)據(jù)點(diǎn)的漸進(jìn)線性收斂比最小，迭代早熟收斂的可能性低，收斂解的可靠性高。當(dāng)然，也可以根據(jù)實(shí)際工程中位移場(chǎng)和應(yīng)力場(chǎng)計(jì)算精度需求和主要的邊界約束類型等靈活取值。

3 平面問(wèn)題中矩陣C 的影響

對(duì)于平面問(wèn)題，平衡系數(shù)從一個(gè)增加到九個(gè)，若按照對(duì)稱正定的要求給出系數(shù)矩陣C，剩余六個(gè)獨(dú)立參數(shù)。

3.1 松弛問(wèn)題的收斂速率

從3.2 節(jié)生成的數(shù)據(jù)庫(kù)(數(shù)據(jù)庫(kù)A)中隨機(jī)選取初始數(shù)據(jù)分配方案Zini，然后，將其固定為每次試驗(yàn)的迭代初值。將D中的每個(gè)元素分別乘上0 到2 的隨機(jī)系數(shù)，從而生成50 組滿足正定對(duì)稱約束的系數(shù)矩陣C用于試驗(yàn)。為減輕量級(jí)變化帶來(lái)的影響，各矩陣均放縮至與D的二范數(shù)相等。圖5給出了第5 和第10 迭代步中，勢(shì)能和余能隨r變化的散點(diǎn)圖和擬合趨勢(shì)。

圖5 第5 和第10 迭代步的勢(shì)能和余能Fig. 5 Potential energy and complementary energy for thefifth and tenth iteration steps

從擬合趨勢(shì)上不難看出，隨r上升，勢(shì)能和余能均呈上升趨勢(shì)，這意味著算式收斂到精確解將需要更長(zhǎng)的迭代步。至此，發(fā)現(xiàn)與m對(duì)應(yīng)力-應(yīng)變可行解收斂速率影響相反所不同的是，隨r的增大，以余能和勢(shì)能度量的應(yīng)力-應(yīng)變解整體收斂速率呈下降趨勢(shì)。

3.2 原始問(wèn)題的收斂性能

為討論當(dāng)數(shù)據(jù)庫(kù)為有限的離散數(shù)據(jù)點(diǎn)時(shí)，r對(duì)收斂性能的影響，構(gòu)建了如下數(shù)據(jù)庫(kù)：離散數(shù)據(jù)點(diǎn)在域內(nèi)均勻分布，應(yīng)力數(shù)據(jù)取值的上下限為有限元參考解最大應(yīng)力水平的1.2 倍。每個(gè)獨(dú)立的應(yīng)力分量上均勻布置有20 個(gè)采樣點(diǎn)，正交組合后共生成8000 個(gè)數(shù)據(jù)點(diǎn)。在生成數(shù)據(jù)庫(kù)的過(guò)程中，應(yīng)變數(shù)據(jù)直接由應(yīng)力數(shù)據(jù)和給定的材料物理方程獲取。生成的數(shù)據(jù)庫(kù)記為數(shù)據(jù)庫(kù)A，其在應(yīng)力空間內(nèi)的投影如圖6 所示。

圖6 數(shù)據(jù)庫(kù)平面投影Fig. 6 Plane projection of database

第一節(jié)中指出算法的本質(zhì)是找到滿足平衡方程和幾何方程，且與整形約束沖突程度最小的近似解。引入罰函數(shù)的數(shù)學(xué)意義是對(duì)數(shù)據(jù)點(diǎn)整形約束的放松，松弛后將形成以各個(gè)數(shù)據(jù)點(diǎn)為中心的凸包，但整個(gè)數(shù)據(jù)庫(kù)松弛后仍然是非凸集合。凸包的形態(tài)決定了跳出局部最優(yōu)的難度，尖銳部分越有可能成為局部陷阱，“捕獲”迭代過(guò)程從而導(dǎo)致算法在局部極小值處收斂，如圖7 所示。

圖7 局部陷阱和自定義歐氏空間圖示Fig. 7 Visual illustration of local trap and custom Euclidean space

表1 基矢量方向上的方差Table 1 Variance in direction of basis vector

圖8 為Zini作為迭代初值時(shí)，應(yīng)力數(shù)據(jù)的各個(gè)獨(dú)立分量與參考解的偏差絕對(duì)值統(tǒng)計(jì)直方圖。結(jié)果表明，方差越大的矢量方向上，收斂點(diǎn)的偏差幅度也越大。原因是方差一定程度上反映了給定數(shù)據(jù)庫(kù)在定義的歐氏空間內(nèi)平均間距的大小，而在平均間距較大的矢量方向上，圖7 所展示的局部陷阱問(wèn)題越嚴(yán)重。A1～A5 中的六個(gè)正交基矢量與應(yīng)力-應(yīng)變空間內(nèi)相應(yīng)分量之間的最大夾角不超過(guò)6°，在后續(xù)分析中近似為平行處理。

在試驗(yàn)A2 中，數(shù)據(jù)庫(kù)A 在Ec2 和Ec3 基矢量方向上的方差最大，反映在計(jì)算結(jié)果(圖8(b)和圖8(c))中就是σy和τxy分量的偏差幅度最大；A3 的試驗(yàn)現(xiàn)象與A2 類似；在試驗(yàn)A4 和A5 中，數(shù)據(jù)庫(kù)A 分別在Ec3 和Ec6 基矢量方向上的方差最大，但是由于應(yīng)力-應(yīng)變數(shù)據(jù)是共同參與迭代的，因此，τxy分量的偏差幅度均增大。

圖8 分配的數(shù)據(jù)點(diǎn)各應(yīng)力分量與參考解之差Fig. 8 Difference of each stress component between assigned data point and reference solution

圖9 為某代表性積分點(diǎn)的應(yīng)力數(shù)據(jù)在松弛問(wèn)題和原始問(wèn)題內(nèi)的迭代路徑對(duì)比。其中標(biāo)號(hào)C11=4C22對(duì)應(yīng)試驗(yàn)A2，標(biāo)號(hào)C11=2C22則是在A2 的基礎(chǔ)上降低了r。在松弛問(wèn)題中，C11=C22對(duì)照組仍沿直線路徑逼近精確解，各個(gè)分量的線性收斂比均為50%。隨C11/C22上升，σy和τxy分量收斂到參考解鄰域所需迭代步數(shù)顯著上升。從原始問(wèn)題的迭代路徑觀察到，數(shù)據(jù)點(diǎn)在σx分量上均收斂至參考解鄰域，但在σy和τxy分量上，隨著C11/C22上升，與參考解的偏差幅度也相應(yīng)上升。

圖9 某代表性積分點(diǎn)的應(yīng)力數(shù)據(jù)迭代路徑Fig. 9 Iterative path of stress data at a representative integral point

為進(jìn)一步說(shuō)明早熟收斂性質(zhì)的影響，由圖4 中相同邊界條件但網(wǎng)格剖分更精細(xì)的有限元模型計(jì)算生成了第二個(gè)數(shù)據(jù)庫(kù)。共包含8044 個(gè)數(shù)據(jù)點(diǎn)，記為數(shù)據(jù)庫(kù)B，其在應(yīng)力空間內(nèi)的投影見(jiàn)圖6。

圖4 平面問(wèn)題模型Fig. 4 Model of plane problem

圖10 為使用不同數(shù)據(jù)庫(kù)驅(qū)動(dòng)仿真計(jì)算時(shí)，各積分點(diǎn)的應(yīng)力-應(yīng)變相對(duì)誤差統(tǒng)計(jì)直方圖，計(jì)算過(guò)程中平衡因子取值均為D。分析表明：精確解附近的數(shù)據(jù)點(diǎn)密度比數(shù)據(jù)庫(kù)容量對(duì)算式精度的影響更大。從圖10 中可以看出，雖然數(shù)據(jù)庫(kù)容量相近，但使用數(shù)據(jù)庫(kù)B 的計(jì)算誤差明顯低于數(shù)據(jù)庫(kù)A。一方面，由于精確解附近的可行解集合被拓展，可能引入了新的全局最優(yōu)解；另一方面，停步準(zhǔn)則中允許出現(xiàn)更小的改變量，減小了在迭代過(guò)程中跳出局部最優(yōu)的難度。合并兩個(gè)數(shù)據(jù)庫(kù)后，計(jì)算誤差仍與單獨(dú)使用數(shù)據(jù)庫(kù)B 相近，說(shuō)明結(jié)合松弛問(wèn)題在迭代末期、參考解鄰域內(nèi)的性質(zhì)來(lái)評(píng)價(jià)原始問(wèn)題的收斂性能是有效可行的。

圖10 積分點(diǎn)誤差統(tǒng)計(jì)Fig. 10 Error statistics of integral points

4 結(jié)論

在線彈性距離最小化數(shù)據(jù)驅(qū)動(dòng)算法的框架下，分別以桁架和平面問(wèn)題為例，討論了平衡因子取值關(guān)于彈性常數(shù)的量級(jí)m和平行關(guān)系偏離程度r對(duì)算法的收斂速率及收斂性能的影響：

(1)在去掉整形約束的松弛問(wèn)題中，對(duì)于給定的m，應(yīng)力可行解和應(yīng)變可行解分別以1/(1+m2)和m2/(1+m2)的收斂比線性收斂，并在迭代終止時(shí)作為近似解輸出。m越大于1，應(yīng)力解的相對(duì)精度更高，m越小于1，應(yīng)變解的相對(duì)精度更高。隨r的增大，以余能和勢(shì)能度量的應(yīng)力-應(yīng)變可行解整體收斂速率均呈下降趨勢(shì)。

(2)在包含整型約束的原始問(wèn)題中，迭代路徑與松弛問(wèn)題的主要區(qū)別在于：應(yīng)力-應(yīng)變數(shù)據(jù)點(diǎn)在相鄰迭代步間的最小改變量受整型約束限制，算法早熟收斂。r為0 時(shí)，應(yīng)力-應(yīng)變數(shù)據(jù)點(diǎn)沿直線逼近參考解，m取1 時(shí)的漸進(jìn)線性收斂比最小，為50%；r不為0 時(shí)，數(shù)據(jù)點(diǎn)在各個(gè)分量方向上收斂到參考解鄰域所需的迭代步長(zhǎng)出現(xiàn)較大差異。對(duì)于收斂較慢的分量，迭代終止時(shí)偏差幅度上升，早熟收斂問(wèn)題愈發(fā)突出。

(3)除了m和r外，邊值約束條件的主要類型、隨機(jī)初值的選取、數(shù)據(jù)點(diǎn)間距也是影響算法收斂性能的因素。其中，精確解附近的數(shù)據(jù)點(diǎn)密度能夠很好解決數(shù)據(jù)點(diǎn)間距和平衡因子取值帶來(lái)的早熟收斂問(wèn)題，從而保障算法的計(jì)算精度。