黏彈性問題的插值型無單元Galerkin方法*

張鵬軒 彭妙娟

(上海大學(xué)土木工程系, 上海 200444)

1 引 言

無網(wǎng)格方法作為一種新的數(shù)值方法, 因其只需要節(jié)點(diǎn)信息, 從而可以避免網(wǎng)格的劃分以及重構(gòu)[1],因此具有前處理簡單、計(jì)算精度高等特點(diǎn), 已經(jīng)成為科學(xué)和工程計(jì)算發(fā)展的重要趨勢之一.當(dāng)前已經(jīng)發(fā)展的無網(wǎng)格方法有無單元Galerkin方法(element-free Galerkin method, EFG)[1-3]、重構(gòu)核粒子方法(reproducing kernel particle method,RKPM)[4-6]、有限點(diǎn)方法 (finite point method,FPM)[7]、無網(wǎng)格局部Petrov-Galerkin方法(meshless local Petrov-Galerkin method, MLPG)[8]、單位分解法[9,10]、復(fù)變量無網(wǎng)格方法[11-15]、維數(shù)分裂無網(wǎng)格方法[16-19]和無網(wǎng)格的邊界積分方程方法[20-26]等.

無單元Galerkin方法是目前應(yīng)用最廣泛的無網(wǎng)格方法之一, 該方法基于移動最小二乘法構(gòu)造逼近函數(shù), 所以具有較高的計(jì)算精度.但該方法的不足之處在于形成方程組時(shí), 有時(shí)是病態(tài)的或者是奇異的.因此, 有可能難以求解或者得到不正確的解.為了改進(jìn)該方法的不足, 程玉民和陳美娟[20]以及Cheng和Peng[21]提出了改進(jìn)的移動最小二乘法.該方法選取正交函數(shù)作為基函數(shù), 從而致使法方程既不病態(tài)也不會奇異, 同時(shí)也不需要求解矩陣的逆, 可以直接得到方程組的解, 降低了計(jì)算量.Zhang等[27-29]采用改進(jìn)的移動最小二乘法建立形函數(shù), 提出了瞬態(tài)熱傳導(dǎo)、波動方程和彈性動力學(xué)等問題改進(jìn)的無單元Galerkin方法.Cheng和Liew[30]以及Cheng和Wei[31]利用改進(jìn)的無單元Galerkin方法求解波動方程和廣義Camassa-Holm方程.文獻(xiàn)[32-35]建立了三維黏彈性、三維彈塑性和彈塑性大變形等問題改進(jìn)的無單元Galerkin方法, 并將改進(jìn)的無單元Galerkin方法用于機(jī)場復(fù)合道面的斷裂力學(xué)分析.Wu等[36]建立了彈性力學(xué)拓?fù)鋬?yōu)化問題改進(jìn)的無單元Galerkin方法.Meng等[37-39]將維數(shù)分裂法與改進(jìn)的無單元Galerkin方法相結(jié)合, 提出了三維勢問題、瞬態(tài)熱傳導(dǎo)問題和波動方程的雜交無單元Galerkin方法.

由于改進(jìn)的移動最小二乘法構(gòu)造的逼近函數(shù)不滿足Kroneckerδ函數(shù)性質(zhì), 使得基于其形成的無網(wǎng)格方法不能直接施加本質(zhì)邊界條件.Lancaster和Salkauskas[40]提出了移動最小二乘插值法, 該方法得到的逼近函數(shù)滿足Kroneckerδ函數(shù)性質(zhì), 從而可以直接施加本質(zhì)邊界條件.

Ren和Cheng[41,42]提出了改進(jìn)的移動最小二乘插值法, 該方法證明了Lancaster的形函數(shù)公式中的一些內(nèi)積可以為零, 可得到更為簡單的形函數(shù)計(jì)算公式, 從而減少了方程個(gè)數(shù)以及未知量個(gè)數(shù),提高了計(jì)算效率.Ren和Cheng[41,42]基于改進(jìn)的移動最小二乘插值法, 建立了勢問題和彈性力學(xué)的插值型無單元Galerkin方法(interpolating elementfree Galerkin method, IEFG).Cheng等[43,44]建立了彈塑性和非線性大變形等問題的插值型無單元Galerkin方法.Deng等[45]建立溫度場問題的插值型復(fù)變量無單元Galerkin方法.

為了克服移動最小二乘插值法因權(quán)函數(shù)奇異導(dǎo)致的計(jì)算困難并避免產(chǎn)生截?cái)嗾`差, Wang等[46]和Sun等[47,48]提出了基于非奇異權(quán)函數(shù)改進(jìn)的移動最小二乘插值法, 建立了勢問題、彈性問題和彈塑性問題改進(jìn)的插值型無單元Galerkin方法.該方法形函數(shù)的待定系數(shù)比傳統(tǒng)的移動最小二乘法少一個(gè), 求逆矩陣的階數(shù)比移動最小二乘法少一階, 提高了計(jì)算效率和計(jì)算精度.基于非奇異權(quán)函數(shù)改進(jìn)的移動最小二乘插值法, Wang等[49,50]研究了兩點(diǎn)邊值問題改進(jìn)的插值型無單元Galerkin方法的誤差估計(jì)和插值型移動最小二乘法的誤差估計(jì).Sun等[51]分析了在n維空間的插值型移動最小二乘法的誤差估計(jì).Liu等[52-54]建立了彈性大變形、彈塑性大變形和凝膠非均勻溶脹大變形等問題改進(jìn)的插值型無單元Galerkin方法.

黏彈性問題在力學(xué)和工程中具有十分重要的應(yīng)用.由于其具有非線性的特點(diǎn), 目前國內(nèi)外主要是運(yùn)用有限元法或者邊界元法等方法來求解黏彈性問題.但是, 近年來無網(wǎng)格方法發(fā)展快速, 越來越多的學(xué)者運(yùn)用此方法來解決黏彈性問題.Yang和Liu[55]把無單元Galerkin方法和時(shí)域精細(xì)算法結(jié)合, 用于求解黏彈性問題.Canelas和Sensale[56]運(yùn)用無網(wǎng)格邊界節(jié)點(diǎn)法分析了彈性和黏彈性材料中的簡諧波問題.彭妙娟和程玉民等建立了改進(jìn)的無單元Galerkin方法[32]、黏彈性問題復(fù)變量無單元Galerkin方法[57]和改進(jìn)的復(fù)變量無單元Galerkin方法[58].

本文基于改進(jìn)的移動最小二乘插值法構(gòu)建形函數(shù), 根據(jù)Galerkin積分弱形式建立方程, 提出了黏彈性問題的插值型無單元Galerkin方法.通過數(shù)值算例討論了影響域、節(jié)點(diǎn)數(shù)對計(jì)算精確性的影響, 說明了該方法具有較好的收斂性; 與無單元Galerkin方法以及有限元的解或解析解比較, 說明了該方法具有提高計(jì)算效率的優(yōu)點(diǎn).

2 改進(jìn)的移動最小二乘插值法

對于改進(jìn)的移動最小二乘插值法[41,42], 定義奇異權(quán)函數(shù)

其中參數(shù)α是正整數(shù),ρI是影響域xI的影響半徑.

定義任意函數(shù)f(x) 和g(x) 的內(nèi)積為

首先對基函數(shù)pi(x) 做標(biāo)準(zhǔn)化處理.在空間span(p1,p2,...,pm)中, 在x點(diǎn)將基函數(shù)p1(x)≡1單位化,

再將p2(x),p3(x),...,pm(x)與q1(x) 正交化,可得

將新的基函數(shù)qi(x) 應(yīng)用于移動最小二乘法可得

將(10)式代入(6)式可得

由于改進(jìn)的移動最小二乘插值法建立的形函數(shù)滿足Kroneckerδ函數(shù)性質(zhì), 因此在建立離散方程時(shí), 可以直接施加本質(zhì)邊界條件, 從而減少了計(jì)算時(shí)間, 提高了計(jì)算效率.

3 黏彈性問題的基本方程

黏彈性問題的應(yīng)力-應(yīng)變關(guān)系與時(shí)間相關(guān), 其材料參數(shù)是時(shí)間的函數(shù).假定σ(t) 在整個(gè)加載過程中恒定, 然而ε(t) 和位移u(t) 隨著時(shí)間變化.如果給定求解域?內(nèi)的體力b, 應(yīng)力邊界Γt上的面力及位移邊界Γu上的位移(邊界條件與時(shí)間無關(guān)),下面給出二維黏彈性問題的基本方程和邊界條件.

平衡方程

其中L是微分算子矩陣

σ是域內(nèi)任意一點(diǎn)的應(yīng)力

b是域內(nèi)任意一點(diǎn)的體力

幾何方程

其中ε和u分別為域內(nèi)任意點(diǎn)的應(yīng)變和位移,

J(t)為不同流變模型的蠕變?nèi)崃? 對Maxwell模型,

對Kelvin 模型,

對三參數(shù)模型,

其中G,G1和G2分別是各模型中彈簧的剪切彈性模量,η是各模型中黏壺的黏性系數(shù),K為體積彈性模量.

邊界條件

n1和n2是邊界點(diǎn)的外法線方向余弦.

由于可以直接施加本質(zhì)邊界條件, 因此黏彈性問題的Galerkin積分弱形式為

4 黏彈性問題的插值型無單元Galerkin方法

設(shè)問題求解域?內(nèi)配有M個(gè)節(jié)點(diǎn), 節(jié)點(diǎn)的影響域?I(I=1,2,···,M)的并集覆蓋了整個(gè)域?.由改進(jìn)的移動最小二乘插值法, 位移向量

黏彈性平面問題的應(yīng)變張量可表示為

黏彈性平面問題的應(yīng)力張量可表示為

是正應(yīng)變矩陣;

是偏應(yīng)變矩陣.

對平面應(yīng)力情況, 即σ33=0 , 可得

其中β1和β2是常數(shù),

則(53)式可寫為

對平面應(yīng)變情況, 即ε33=0 , 可得

將(49)式和(50)式代入(48)式得到

將(44)式和(61)式代入(35)式得到

由δUT的任意性, 得到最后的離散系統(tǒng)方程為

(63)式中, 對于平面應(yīng)力問題,

對于平面應(yīng)變問題,

以上即為黏彈性問題的插值型無單元Galerkin方法.

5 數(shù)值算例

為了驗(yàn)證上述理論的準(zhǔn)確性, 以下通過該方法對4個(gè)算例進(jìn)行了計(jì)算, 并將計(jì)算結(jié)果與無單元Galerkin方法以及有限元方法的計(jì)算結(jié)果或解析解進(jìn)行對比.在算例中, 構(gòu)造移動最小二乘插值法的形函數(shù)采用線性基函數(shù)和奇異權(quán)函數(shù).在每個(gè)積分單元, 采用 4×4 的Gauss積分.

定義相對誤差公式為

其中L2范數(shù)定義為

定義方差公式為

算例1受均布荷載的懸臂梁

圖1所示的受均布荷載作用的懸臂梁, 梁長L=8m , 梁高D=2m , 取單位厚度, 按平面應(yīng)力計(jì)算.材料體積變化是彈性的,E=1.0×108Pa ,v=0.25, 剪切變形流變性質(zhì)滿足Kelvin黏彈性模型, 其參數(shù)為G=2.0×108Pa ,η=6.0×108Pa·s ,均布荷載為p=3.0×104Pa , 不計(jì)自重.

圖1 受均布荷載的懸臂梁Fig.1.A cantilever beam subjected to a distributed loading.

為保證有限元法(FEM)數(shù)值解的精確性, 進(jìn)行不同節(jié)點(diǎn)布置以得到較為精確的數(shù)值解.采用四邊形4節(jié)點(diǎn)單元, 具體節(jié)點(diǎn)布置為: 1 7×9 ,21×11, 2 9×11 , 3 3×11 , 3 3×13 , 3 7×17 , 圖2給出了不同節(jié)點(diǎn)分布下有限元法解的方差.從圖2可以看出, 節(jié)點(diǎn)分布越密方差越小, 說明有限元法得到的解是收斂的.在此選用 3 3×13 的節(jié)點(diǎn)布置.

圖2 不同節(jié)點(diǎn)分布下有限元法解的方差Fig.2.The variances of the solutions of FEM under different node distributions.

在不同節(jié)點(diǎn)分布 9×3 , 1 1×4 , 1 3×5 , 1 5×6 ,17×7和 1 9×8 的情況下, 利用本文提出的插值型無單元Galerkin方法計(jì)算, 圖3給出了不同節(jié)點(diǎn)分布下數(shù)值解的相對誤差.從圖3可以看出, 節(jié)點(diǎn)分布越密相對誤差越小, 即插值型無單元Galerkin方法具有較好的收斂性.在考慮計(jì)算效率的情況下, 采用圖4所示的 1 7×7 節(jié)點(diǎn)布置, 背景積分網(wǎng)格選取 1 6×6.

圖3 不同節(jié)點(diǎn)分布下的相對誤差Fig.3.The relative error under different node distributions.

圖4 節(jié)點(diǎn)布置Fig.4.Node distribution.

不同影響域比例參數(shù)對本文方法數(shù)值解的精度具有一定影響, 如圖5所示, 當(dāng)dmax=1.7-1.9時(shí), 該方法數(shù)值解的精度較高.本算例選用dmax=1.7.

圖5 不同影響域比例參數(shù)下的相對誤差Fig.5.The relative error for different scale parameters of influence domains.

采用插值型無單元Galerkin方法求解時(shí), 節(jié)點(diǎn)布置為 1 7×7 , 背景積分網(wǎng)格選取 1 6×6 ,dmax=1.7, 此時(shí)中軸線各點(diǎn)撓度和右端中點(diǎn)撓度隨時(shí)間變化的總相對誤差為 1.29% , 計(jì)算時(shí)間為31.19 s.

采用無單元Galerkin方法求解時(shí), 節(jié)點(diǎn)布置為 1 7×7 , 背景積分網(wǎng)格選取 1 6×6 ,dmax=3.0 ,α=2.3×1010, 從而可以得到與插值型無單元Galerkin方法相近的計(jì)算精度, 計(jì)算時(shí)間為34.67 s.

上述兩種方法所得的數(shù)值解與有限元法的數(shù)值解的對比如圖6和圖7所示, 可以看出, 兩種方法得到的數(shù)值解與有限元法的數(shù)值解較為符合, 但插值型無單元Galerkin方法能夠提高計(jì)算效率.

圖6 t=20s 時(shí)懸臂梁中軸線上各點(diǎn)的撓度Fig.6.Vertical displacements of nodes on the neutral axis of the beam when t=20s.

圖7 梁右端中點(diǎn)的撓度隨時(shí)間的變化Fig.7.Time history of vertical displacement of midpoint in the right end of the beam.

算例2受純彎曲的梁

圖8為受純彎曲的梁, 幾何參數(shù)為L=5m ,H=2m, 單位厚度.材料體積變化是彈性的,E=1.0×106Pa ,ν=0.3 , 剪切變形的流變性質(zhì)滿足三參數(shù)模型, 其參數(shù)為G1=5.0×105Pa ,G2=1.0×106Pa 和η=2.0×106Pa·s.其所受三角形分布荷載大小如圖8所示, 不計(jì)體力, 按平面應(yīng)力問題計(jì)算.

圖8 純彎曲的梁Fig.8.A beam subjected to simple bending.

由于對稱性和反對稱性, 只取1/4梁進(jìn)行計(jì)算.在x1軸和x2軸上, 位移u1均為 0, 為了消除結(jié)構(gòu)在x2方向上剛體位移, 可令原點(diǎn)處u2為0.有限元法采用四邊形4節(jié)點(diǎn)單元, 節(jié)點(diǎn)布置為 2 1×13.無單元Galerkin方法和插值型無單元Galerkin方法均采用圖9所示的 1 5×7 的節(jié)點(diǎn)布置, 背景積分網(wǎng)格選取 1 4×6.

圖9 節(jié)點(diǎn)分布Fig.9.Node distribution.

采用插值型無單元Galerkin方法求解時(shí),dmax=2.0, 此時(shí)中軸線各點(diǎn)撓度和右端中點(diǎn)撓度隨時(shí)間變化的總相對誤差為 0.30% , 計(jì)算時(shí)間為21.88 s.

時(shí)采用無單元Galerkin方法求解時(shí),dmax=3.0 ,α=1.5×108, 從而可以得到與插值型無單元Galerkin方法相近的計(jì)算精度, 計(jì)算時(shí)間為28.76 s.

上述兩種方法所得數(shù)值解與有限元法數(shù)值解的對比如圖10和圖11所示.可以看出, 插值型無單元Galerkin方法和無單元Galerkin方法都可以達(dá)到較高的精度, 但插值型無單元Galerkin方法能夠提高計(jì)算效率.

算例3受均布內(nèi)壓的圓環(huán)

受均布內(nèi)壓的圓環(huán)如圖12所示, 其內(nèi)表面承受均勻分布的壓力p=30kPa.幾何參數(shù)為a=1m ,b=5m.材料體積變化是彈性的,E=1.0×107Pa ,ν=0.25, 剪切變形的流變性質(zhì)采用Kelvin模型.其參數(shù)為G=5.0×105Pa ,η=2.0×106Pa·s.不計(jì)體力, 按平面應(yīng)變問題計(jì)算.

圖10 t=30s 時(shí)梁中軸線上的節(jié)點(diǎn)撓度Fig.10.Vertical displacements of nodes on the neutral axis of the beam when t=30s.

圖11 梁右端中點(diǎn)的撓度隨時(shí)間 t 的變化Fig.11.Time history of vertical displacement of midpoint in the right end of the beam.

圖12 受均布內(nèi)壓的厚壁圓筒Fig.12.Circular ring under a distributed inner pressure.

在極坐標(biāo)系下, 采用Kelvin模型, 圓筒徑向位移隨時(shí)間變化的解析解為

其中

(70)式表明Kelvin模型加載瞬時(shí)位移為零, 當(dāng)時(shí)間足夠長, 位移將趨于穩(wěn)定.

如圖13, 由于對稱性, 只取1/4區(qū)域?yàn)檠芯繉ο?無單元Galerkin方法和插值型無單元Galerkin方法均選用如圖14所示的的節(jié)點(diǎn)布置, 背景積分網(wǎng)格選取.

圖13 受均布內(nèi)壓1/4圓筒Fig.13.A quarter of the circular ring under a distributed inner pressure.

圖14 1/4圓筒的節(jié)點(diǎn)分布Fig.14.Node distribution of a quarter of the circular ring.

采用插值型無單元Galerkin方法求解時(shí),dmax=1.01, 此 時(shí)t=30s 時(shí) 沿x2=0 線上節(jié)點(diǎn)的位移和點(diǎn) (2,0) 的徑向位移隨時(shí)間變化的總相對誤差為 1.10% , 計(jì)算時(shí)間為23.38 s.

采用無單元Galerkin方法求解時(shí),dmax=1.5 ,α=8.0×107, 此時(shí)可得到與上述插值型無單元Galerkin方法相近的相對誤差, 計(jì)算時(shí)間為28.66 s.

兩種方法得到的數(shù)值解與解析解的對比如圖15和圖16所示.從以上分析和圖15和圖16的對比可看出, 插值型無單元Galerkin方法和無單元Galerkin方法的數(shù)值解均與解析解符合得較好,但插值型無單元Galerkin方法能夠提高計(jì)算效率.

圖15 t=30s 時(shí)沿 x2=0 線上節(jié)點(diǎn)的位移Fig.15.Radial displacements at x2=0 when t=30s.

圖16 點(diǎn) (2,0) 的徑 向位移隨時(shí)間 t 的變化Fig.16.Time history of radial displacement at point (2,0).

算例4工程算例: 受靜水壓力的混凝土水壩

為了證明本方法在較復(fù)雜幾何求解域的可行性, 本算例計(jì)算一個(gè)受到靜水壓力的混凝土水壩,水壩的幾何參數(shù)如圖17所示.壩高100 m, 上游水位90 m, 下游無水.壩底部完全固定, 僅考慮壩體受到的上游水荷載.壩體材料為混凝土, 所以需要考慮建造完工后期混凝土的流變效應(yīng).假定混凝土材料體積變化是彈性的,E=2.0×1010Pa ,ν=0.2 ,混凝土的早期流變性質(zhì)滿足三參數(shù)模型, 其參數(shù)為G1=8.33×109Pa ,G2=4.0×1010Pa 和η=8.0×1016Pa·s.不計(jì)體力, 按平面應(yīng)變問題計(jì)算.

圖17 受靜水壓力的混凝土水壩Fig.17.A concrete dam under hydrostatic pressure.

有限元法采用四邊形4節(jié)點(diǎn)單元, 共布置234個(gè)節(jié)點(diǎn).無單元Galerkin方法和插值型無單元Galerkin方法均在求解域內(nèi)布置如圖18所示的255個(gè)節(jié)點(diǎn), 背景積分網(wǎng)格選取231個(gè).

圖18 混凝土水壩的節(jié)點(diǎn)分布Fig.18.Node distribution of a concrete dam.

采用插值型無單元Galerkin方法求解時(shí),dmax=1.2,t=500 d時(shí)沿x1=15 方向上節(jié)點(diǎn)的水平位移和點(diǎn) (15,50) 的水平位移與時(shí)間關(guān)系的總相對誤差為 0.27% , 計(jì)算時(shí)間為50.42 s.

采用無單元Galerkin方法求解時(shí),dmax=3.0 ,α=2.0×1012, 其總相對誤差為 1.36% , 計(jì)算時(shí)間為94.44 s.

兩種方法的數(shù)值解與有限元法的數(shù)值解的對比如圖19和圖20所示.從以上分析和圖19和圖20的對比可以看出, 與無單元Galerkin方法相比, 插值型無單元Galerkin方法具有更高的精度和效率, 說明該方法在解決復(fù)雜工程問題時(shí)能夠提高計(jì)算精度和計(jì)算效率.

圖19 t=500 d時(shí)沿 x1=15 方向上節(jié)點(diǎn)的水平位移Fig.19.Horizontal displacements at x1=15 when t=500d.

圖20 混凝土壩上點(diǎn) (15,50) 的水平位移與時(shí)間的關(guān)系Fig.20.Time history of horizontal displacement of the point (15,50).

6 結(jié) 論

本文基于改進(jìn)的移動最小二乘插值法構(gòu)造形函數(shù), 建立了黏彈性問題的插值型無單元Galerkin方法.由于改進(jìn)的移動最小二乘插值法的形函數(shù)滿足Kroneckerδ函數(shù)性質(zhì), 避免了傳統(tǒng)無單元Galerkin方法要利用Lagrange乘子法或者罰函數(shù)法來施加本質(zhì)邊界條件, 大大減少了計(jì)算量, 從而有效地提高了在求解黏彈性問題時(shí)的計(jì)算效率.通過數(shù)值算例討論了影響域大小、節(jié)點(diǎn)數(shù)對計(jì)算精確性的影響, 說明了該方法具有較好的收斂性; 將計(jì)算結(jié)果和有限元解或解析解進(jìn)行對比, 說明了插值型無單元Galerkin方法較無單元Galerkin方法在求解黏彈性問題上具有提高計(jì)算效率的優(yōu)點(diǎn).