黏彈性模型微分型本構(gòu)方程的矩陣形式1)

彭培火 黃朝軍

*(北京建筑大學(xué)理學(xué)院,北京 102612)

?(中國(guó)建筑第二工程局有限公司,北京 100160)

黏彈性材料廣泛應(yīng)用于各類工程中,如塑料、橡膠等工業(yè)材料,土壤、瀝青等地質(zhì)材料,肌肉、血液等生物材料,不勝枚舉。由于黏彈性材料兼具彈性與黏性兩種特性,其力學(xué)性能依賴于環(huán)境溫度、加載時(shí)間、加載速率和幅值等條件,因此建立合適的模型,準(zhǔn)確地描述材料的黏彈性行為,對(duì)工程應(yīng)用及理論研究都具有十分重要的意義。材料黏彈性特性的描述和本構(gòu)方程的表達(dá)是黏彈性力學(xué)的重要研究?jī)?nèi)容之一。為了描述黏彈性材料的應(yīng)力?應(yīng)變規(guī)律,常利用彈性元件和黏性元件以及它們的各種組合構(gòu)成的模型來進(jìn)行模擬。如Maxwell 模型(圖1(a))、Kelvin 模型(圖1(b))、Poynting–Thomson 模型(圖1(c))、Burgers 模型(1(d))、廣義Maxwell 模型、廣義Kelvin 模型等等。多年來,人們開發(fā)了各種各樣的模型來表征材料的黏彈性特性,對(duì)黏彈性材料的本構(gòu)方程的數(shù)學(xué)描述進(jìn)行了大量的研究[1]。利用這些模型,諸多學(xué)者在各個(gè)工程領(lǐng)域開展了大量的研究。例如,將土體視為黏彈性材料,鄭灶鋒等[2]利用Merchant 模型研究了穩(wěn)態(tài)載荷頻率和黏滯性對(duì)地基振動(dòng)的影響,研究電站地下廠房洞室圍巖的變形規(guī)律時(shí),楊林德等[3]基于Kelvin–Voigt 模型來模擬和計(jì)算圍巖的黏彈性變形特征。范家參[4]采用Poynting–Thomson 模型對(duì)斷裂動(dòng)力學(xué)問題進(jìn)行了研究。Ren 等[5]利用改進(jìn)的Burgers 模型,探究了黏彈性復(fù)合材料結(jié)構(gòu)的振動(dòng)特征。Huang 等[6]利用改進(jìn)的Burgers 模型來描述瀝青混凝土的黏彈性本構(gòu)行為。為了更準(zhǔn)確地表征血液的流變特性,更好地與高切變率下的實(shí)驗(yàn)數(shù)據(jù)相吻合,人們也常常采用由更多元件組成的黏彈性模型來進(jìn)行相關(guān)研究[7]。候之超等[8]利用分?jǐn)?shù)微分方程描述泡沫材料的黏彈性特征,研究了單向恒速率加載下模型的力學(xué)響應(yīng)。李兆生等[9]利用廣義Maxwell 模型構(gòu)建損傷力學(xué)本構(gòu)模型,研究了瀝青材料的凍融損傷過程。詹小麗等[10]采用分?jǐn)?shù)階導(dǎo)數(shù)Maxwell 模型對(duì)瀝青的動(dòng)態(tài)黏彈性力學(xué)性能進(jìn)行了研究。潘文瀟等[11]研究了兩平板之間的廣義Maxwell 黏彈性流體非定常流動(dòng)的解析解。不同的模型具有不同的優(yōu)缺點(diǎn),如Maxwell 模型可以模擬應(yīng)力松弛現(xiàn)象,但不能表示蠕變過程；Kelvin模型能夠模擬蠕變,卻不能表示應(yīng)力松弛[12-13]。一般來說,隨著彈簧和阻尼器數(shù)量的增多,模型能夠模擬的黏彈性特征也會(huì)越來越貼近實(shí)際材料[14]。但是,對(duì)于多個(gè)彈簧和阻尼器以任意連接方式組成的黏彈性模型,其微分型本構(gòu)方程如何快速計(jì)算出,是各種基于黏彈性材料結(jié)構(gòu)計(jì)算的基礎(chǔ)。圖論算法可方便快捷地計(jì)算路徑等問題,而有向圖的路徑和閉包圍又與黏彈性模型的應(yīng)力方程和應(yīng)變方程存在聯(lián)系,以計(jì)算由10 個(gè)彈性元件和黏性元件構(gòu)成的較為復(fù)雜的模型的微分本構(gòu)方程為例,推導(dǎo)得到了任意線性黏彈性模型的微分型本構(gòu)方程的一般矩陣形式,總結(jié)了一套適合計(jì)算機(jī)編程的固定范式。并利用Python語(yǔ)言編寫了圖形界面用戶程序,可快速地構(gòu)建任意復(fù)雜的線性黏彈性模型和計(jì)算其應(yīng)力?應(yīng)變微分關(guān)系式。該計(jì)算方法對(duì)準(zhǔn)確地構(gòu)建實(shí)際工程材料的黏彈性計(jì)算模型、精確地?cái)M合實(shí)際工程材料的黏彈性特性及應(yīng)力?應(yīng)變規(guī)律具有一定的參考意義。

1 本構(gòu)方程的矩陣形式表述

一般的黏彈性力學(xué)教材中都給出了各種基本模型的微分型本構(gòu)方程[14],如：Maxwell 模型、Kelvin模型、Poynting–Thomson 模型、Burgers 模型等等。由于彈性元件和黏性元件數(shù)量不多,它們的本構(gòu)方程推導(dǎo)過程并不難。為了實(shí)現(xiàn)計(jì)算機(jī)編程計(jì)算本構(gòu)方程,以下利用圖論算法和矩陣運(yùn)算,通過推導(dǎo)由10 個(gè)彈性元件和黏性元件構(gòu)成的較為復(fù)雜模型的應(yīng)力?應(yīng)變微分關(guān)系,總結(jié)和建立任意線性黏彈性模型的微分型本構(gòu)方程的一般形式。

如圖2(a) 所示的10 元件黏彈性模型,可以看成由Maxwell 單元連接組成,如圖2(b) 中所示。單個(gè)Maxwell 單元的應(yīng)力?應(yīng)變微分關(guān)系為

圖2 Fig.2

基于公式(1),只要寫出10 元件黏彈性模型的應(yīng)變?chǔ)拧⒛Ｐ偷膽?yīng)力σ與各個(gè)Maxwell 單元的應(yīng)變?chǔ)舏(i= 1,2,··· ,5),應(yīng)力σi(i= 1,2,··· ,5) 之間的關(guān)系式,將它們適當(dāng)求導(dǎo),并寫成矩陣形式,模型的微分型本構(gòu)方程即可以通過各系數(shù)矩陣的運(yùn)算得到。為了計(jì)算便捷和程式化,借助圖論中有向圖的相關(guān)工具,可迅速將求解過程編寫成程序?，F(xiàn)將整個(gè)求解過程描述如下。

(1)根據(jù)黏彈性模型中的彈簧和阻尼器的連接方式,寫出模型的總應(yīng)變?chǔ)排c各個(gè)Maxwell 單元的應(yīng)變?chǔ)舏之間的關(guān)系式,并求j階導(dǎo)數(shù)。將模型的總應(yīng)變?chǔ)偶捌鋵?dǎo)數(shù)與各Maxwell 單元的應(yīng)變?chǔ)舏及其導(dǎo)數(shù)的關(guān)系式以矩陣形式記為

式中,記Djε表示由模型的總應(yīng)變?chǔ)偶捌?～j階導(dǎo)數(shù)組成的列向量,有j+ 1 個(gè)元素,形如式(3)所示

記Djεi表示各個(gè)Maxwell 單元的應(yīng)變?chǔ)舏及其1～j階導(dǎo)數(shù)組成的列向量,有i×(j+1) 個(gè)元素,i為Maxwell 單元的數(shù)量,形如式(4) 所示

系數(shù)矩陣A有(j+1)×r行,i×(j+1) 列,為行滿秩矩陣,即R(A) = (j+1)r。Ir,j+1是形如式 (5) 所示的矩陣。因模型的總應(yīng)變?chǔ)排c各個(gè)Maxwell 單元的應(yīng)變?chǔ)舏之間的關(guān)系式共有r個(gè)獨(dú)立方程,故乘以該Ir,j+1矩陣來將其合寫成一個(gè)矩陣方程。

式中,dr為共有r個(gè)元素且值全為1 的向量。將模型轉(zhuǎn)化為有向圖,如圖3 所示。模型的總應(yīng)變?chǔ)排c各個(gè)基礎(chǔ)元件的應(yīng)變?chǔ)舏之間的r個(gè)獨(dú)立應(yīng)變關(guān)系式可以根據(jù)由Maxwell 單元組成的模型中由起點(diǎn)到終點(diǎn)的獨(dú)立路徑的數(shù)量進(jìn)行確定。有向圖中獨(dú)立路徑的數(shù)量即為r,如圖3(a) 中所示,r= 3?？梢詫懗? 個(gè)總應(yīng)變?chǔ)排c各個(gè)Maxwell 單元的應(yīng)變?chǔ)舏之間的關(guān)系式,即

圖3 Fig.3

根據(jù)圖論中求解有向圖路徑的方法,先寫出圖的關(guān)聯(lián)矩陣為

列出方程Mx=b,然后求其解(x的元素只能是1 或者0,向量b表示從起點(diǎn)v1到終點(diǎn)v4)。

方程(10) 的滿足條件的解中,有3 個(gè)線性無關(guān)的解向量為

以式(11) 中的每一個(gè)解表示一條從起點(diǎn)v1到終點(diǎn)v4的路徑,解的第幾個(gè)元素為1 表示該路徑經(jīng)過第幾條有向邊。如：向量x2的第1,4,5 行元素為1,其余元素為0,表示圖3(a)中路徑 2○是一條經(jīng)過e1,e4,e5從起點(diǎn)v1到終點(diǎn)v4的路徑。由x1,x2,x3構(gòu)成的矩陣即為獨(dú)立路徑矩陣,記為R。

對(duì)比式(2)～式(8),式(12),可以得出式(2)中的矩陣A可由獨(dú)立路徑矩陣R構(gòu)造而成。

(2)根據(jù)黏彈性模型中的彈簧和阻尼器的連接方式,寫出模型的總應(yīng)力σ與各個(gè)Maxwell 單元的應(yīng)力σi之間的關(guān)系式,并求m階導(dǎo)數(shù),將模型的總應(yīng)力σ及其導(dǎo)數(shù)與各Maxwell 單元的應(yīng)力σi及其導(dǎo)數(shù)的關(guān)系式以矩陣形式記為

式中,模型的總應(yīng)力σ與各個(gè)基礎(chǔ)元件的應(yīng)力σi之間的關(guān)系式共有n個(gè)獨(dú)立方程,利用乘以In,m+1來將其合寫成一個(gè)矩陣方程。記Dmσ表示由模型的總應(yīng)力σ及其1～m階導(dǎo)數(shù)組成的列向量,有m+1個(gè)元素,形如式(15) 所示。

記Dmσi表示各個(gè)Maxwell 單元的應(yīng)力σi及其1～m階導(dǎo)數(shù)組成的列向量,有i×(m+1)個(gè)元素,i為Maxwell 單元的數(shù)量,形如式(16) 所示。

系數(shù)矩陣B有(m+1)×n行,i×(m+1) 列,為行滿秩矩陣,即R(B)=(m+1)n。模型的總應(yīng)力σ與各個(gè)基礎(chǔ)元件的應(yīng)力σi之間的n個(gè)獨(dú)立關(guān)系式可以根據(jù)有向圖中包含起點(diǎn)v1的閉包圍的數(shù)量進(jìn)行確定(其中閉包圍須滿足條件：與閉包圍相交的有向邊的方向應(yīng)指向外)。如圖3(b)中所示,n=3?？梢詫懗? 個(gè)總應(yīng)力σ與各個(gè)Maxwell 單元的應(yīng)力σi之間的應(yīng)力關(guān)系式,即

因?yàn)橛邢驁D的任一閉包圍必與每一條路徑相交,故閉包圍矩陣的任一行與獨(dú)立路徑矩陣R的每一行求點(diǎn)積必為1,于是閉包圍矩陣可通過求方程Rx=c的解獲得(非齊次項(xiàng)c為全1 向量,其行數(shù)與R的行數(shù)相同,解x的元素只能是1 或者0)。

方程(20) 的滿足條件的解中,有3 個(gè)線性無關(guān)的解向量為

以上式(21)中的每一個(gè)解表示有向圖中一條包含起點(diǎn)v1的閉包圍,解的第幾個(gè)元素為1 表示該閉包圍與第幾條有向邊相交。如：向量x2的第2,3,4 行元素為1,表示圖3(b)中閉包圍 2○與有向邊e2,e3,e4相交。由以上解x1,x2,x3構(gòu)成的矩陣即為閉包圍矩陣,記為E。

對(duì)比式(14)～式(19),式(22),可以得出式(14)中的矩陣B可由閉包圍矩陣E構(gòu)造而成。

(3) 將Maxwell 單元的應(yīng)力?應(yīng)變微分關(guān)系式(1) 求適當(dāng)階數(shù)的導(dǎo)數(shù)后代入式(2),可得到模型的總應(yīng)變?chǔ)排c各個(gè)Maxwell 單元的應(yīng)力σi之間的關(guān)系式,以矩陣形式記為

式中,記Djε′表示對(duì)向量Djε再求一階導(dǎo)數(shù),故Djε′是由模型的總應(yīng)變?chǔ)诺?～j+1 階導(dǎo)數(shù)組成的列向量,有j+1 個(gè)元素,形如式(25) 所示。

由黏彈性模型中各彈簧的彈性模量Ei及各阻尼器的黏度系數(shù)ηi構(gòu)造Y矩陣和V矩陣分別為

不難得出,式(24) 中的系數(shù)矩陣H為

(4) 聯(lián)立求解方程(14) 與(24),消去中間變量Dmσi,即可求得模型的應(yīng)力σ與應(yīng)變?chǔ)胖g的微分關(guān)系式(需滿足m=j+1,保證Dmσi=Dj+1σi)。具體求解過程如下。

將式(14) 看成是非齊次線性方程組,將向量Dmσi看成待求量,In,m+1Dmσ為非齊次量,則Dmσi可由Dmσ表示為

式中,null(B) 表示矩陣B的零空間矩陣,P為式(14) 的任意一特解,可按式(30) 方法構(gòu)造

其中矩陣F為

式中,列向量f可取獨(dú)立路徑矩陣R的任一行通過轉(zhuǎn)置得到,F矩陣有i×(m+1) 行,m+1 列。將式(29) 代入式(24),可得

記

求出CT的零空間矩陣 null(CT)。 式 (33)中,null(B) 為i ×(m+ 1) 行,(i ?n)×(m+ 1)列矩陣,H為r ×(j+ 1) 行,i ×(j+ 2) 列矩陣(因?yàn)閙=j+ 1,H ×null(B) 符合矩陣乘法規(guī)則)。C為r ×(j+1) 行,(i ?n)×(m+1) 列矩陣,故CT的零空間矩陣null(CT) 有r ×(j+ 1)行,r×(j+1)?(i ?n)×(m+1) 列。記null(CT)的某一個(gè)列向量為b,將其轉(zhuǎn)置并左乘式(32),則可以得到黏彈性模型的應(yīng)力σ與應(yīng)變?chǔ)胖g的微分關(guān)系式的矩陣形式表述為

(5)上述計(jì)算過程中,關(guān)于r個(gè)應(yīng)變方程及n個(gè)應(yīng)力方程求多少階導(dǎo)數(shù)最合適的問題,即m和j的取值為多少合適。從上述的應(yīng)變方程和應(yīng)力方程以及推導(dǎo)過程可知,應(yīng)力的求導(dǎo)階數(shù)比應(yīng)變應(yīng)該高一階,故應(yīng)有m=j+1。為了使CT的零空間矩陣null(CT) 至少含有一個(gè)向量b,則應(yīng)有

式中,i為黏彈性模型中Maxwell 單元的數(shù)量。

(6) 為了將上述計(jì)算方法程式化, 需要將黏彈性模型中的所有彈性元件和黏性元件都配對(duì)成Maxwell 基本單元。當(dāng)由彈簧和阻尼器構(gòu)成的連接圖中存在有非Maxwell 基本單元時(shí),先把非Maxwell單元的部分?jǐn)U展成Maxwell 單元。例如,圖1(c) 所示的Poynting–Thomson 模型,可以通過添加相應(yīng)的彈簧和阻尼器進(jìn)行擴(kuò)展成如圖4(a) 所示,則該模型中的非Maxwell 單元的部分都擴(kuò)展成了Maxwell 單元,如圖4(b) 所示。先計(jì)算出圖4 所示的黏彈性模型的微分本構(gòu)方程,然后令η2為∞,即與原模型等價(jià)。

圖4 Fig.4

2 Python 計(jì)算程序示例

基于以上推導(dǎo)的計(jì)算方法,利用Python 編寫了圖形用戶界面程序(如圖5 所示),可以實(shí)現(xiàn)自由搭建模型和計(jì)算本構(gòu)方程等功能。如將軟件界面右側(cè)的彈簧或阻尼器拖入左側(cè)圖框中,即可構(gòu)造以任意的連接方式組成的線性黏彈性模型。程序在鼠標(biāo)拖拽的同時(shí)會(huì)自動(dòng)識(shí)別各個(gè)彈性元件或黏性元件的連接結(jié)構(gòu),若自動(dòng)識(shí)別的連接結(jié)構(gòu)不符合要求,還可以通過<編輯>菜單進(jìn)行調(diào)整。黏彈性模型構(gòu)建完畢后即可轉(zhuǎn)換得到與之相對(duì)應(yīng)的有向圖,并得到關(guān)聯(lián)矩陣M。在<計(jì)算>菜單中可設(shè)置相關(guān)參數(shù)并進(jìn)行計(jì)算,如采用符號(hào)計(jì)算或是數(shù)值計(jì)算,零空間矩陣的求解方法等等。程序按照?qǐng)D6 所示的流程圖可快速計(jì)算出黏彈性模型的微分形式的本構(gòu)方程。并在此基礎(chǔ)上,由本構(gòu)方程計(jì)算松弛模量、蠕變?nèi)岫?、?fù)模量等描述黏彈性材料的基本物理量。計(jì)算結(jié)果輸出在程序界面下方的圖框內(nèi),也可輸出中間計(jì)算結(jié)果(如獨(dú)立路徑矩陣R、閉包圍矩陣E等等)。通過多個(gè)模型的測(cè)試,利用該程序均能快速正確地計(jì)算出黏彈性模型的本構(gòu)方程。其中表1 為2 個(gè)簡(jiǎn)單黏彈性模型的算例,依據(jù)圖6 所示的計(jì)算流程,給出了中間過程及最終結(jié)果。程序還可繼續(xù)擴(kuò)展,開發(fā)具有其他功能的代碼。

圖5 程序界面及10 元件模型的計(jì)算結(jié)果Fig.5 Program interface and results of 10-component model

圖6 線性黏彈性模型本構(gòu)方程計(jì)算流程圖Fig.6 Calculation flow chart of constitutive equation of linear viscoelastic model

表1 中右側(cè)所示模型,若令η2和E3為無窮大,表示阻尼η2處和彈簧E3處為剛性連接,即與圖1(d)所示的Burgers 模型等價(jià)。按照式(34) 計(jì)算得到的本構(gòu)方程中,η2和E3都出現(xiàn)在分母的位置,令η2和E3為無窮大,則倒數(shù)為零,可以很方便計(jì)算得到Burgers 模型的本構(gòu)方程為

表1 兩個(gè)簡(jiǎn)單模型的計(jì)算示例(續(xù))Table 1 Calculation examples of two simple models (continued)

3 結(jié)語(yǔ)

諸多工程材料如塑料、橡膠、瀝青等具有隨時(shí)間變化的力學(xué)響應(yīng),如恒載荷下的蠕變、恒定變形下的應(yīng)力松弛和卸載時(shí)的延遲應(yīng)變恢復(fù)等現(xiàn)象。材料的這些黏彈性特征對(duì)它們承載時(shí)的力學(xué)性能有著重要的影響。為了從數(shù)學(xué)上對(duì)黏彈性材料的應(yīng)力?應(yīng)變規(guī)律進(jìn)行描述,常利用彈性元件和黏性元件以及它們的組合來建立模型,模擬材料的黏彈性行為。一般來說,彈簧和阻尼器數(shù)量越多,模型能夠模擬的黏彈性特征也會(huì)越來越貼近實(shí)際材料。但是,元件越多,模型的本構(gòu)方程的計(jì)算也越復(fù)雜。對(duì)于多個(gè)彈簧和阻尼器以任意連接方式組成的黏彈性模型,本文推導(dǎo)出了其微分型本構(gòu)方程的一般矩陣形式的計(jì)算方法。首先通過添加彈性參數(shù)或黏性參數(shù)為無窮大的虛加元件,將存在非Maxwell 基本單元的模型擴(kuò)展為由Maxwell 基本單元構(gòu)成的標(biāo)準(zhǔn)模型。然后將標(biāo)準(zhǔn)模型轉(zhuǎn)化為與之對(duì)應(yīng)的有向圖,根據(jù)圖論的有關(guān)知識(shí),基于有向圖的特征,由獨(dú)立路徑的形式可以寫出基本應(yīng)變方程,由閉包圍的形式可以寫出基本應(yīng)力方程,并根據(jù)式(35) 對(duì)其求適當(dāng)階數(shù)的導(dǎo)數(shù),按照?qǐng)D6 所示的計(jì)算流程可以得到獨(dú)立路徑矩陣、閉包圍矩陣及其他中間結(jié)果,通過各系數(shù)矩陣的運(yùn)算,即可根據(jù)式(34)直接得到模型的微分型本構(gòu)方程的矩陣形式表述。本文總結(jié)了一套適合計(jì)算機(jī)編程的固定范式,編寫了Python 程序,可方便搭建黏彈性模型及快速計(jì)算其本構(gòu)關(guān)系,并在此基礎(chǔ)上可繼續(xù)開發(fā)具有其他功能的代碼。