移動最小二乘法在高階連續(xù)問題數(shù)值模擬中的應(yīng)用

陳根生， 楊子勝， 孫玉周

(中原工學(xué)院, 鄭州 450007)

摘要：以彎曲梁為例，應(yīng)用移動最小二乘法構(gòu)造具有高階連續(xù)特征的形函數(shù)，建立無網(wǎng)格數(shù)值計算框架。研究移動最小二乘法在傳統(tǒng)彎曲梁和高階連續(xù)彎曲梁數(shù)值模擬中的應(yīng)用，并對施加邊界條件等相關(guān)關(guān)鍵問題進(jìn)行討論。

關(guān)鍵詞：彎曲梁；高階連續(xù)；移動最小二乘法；無網(wǎng)格法；罰數(shù)法

中圖分類號：TU53

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

DOI:10.3969/j.issn.1671-6906.2015.03.019

Abstract:The shape function with the high-order continuum is constructed with the moving least-square method, and a mesh-free computational scheme is then established for the bending beams. The proposed method is used to investigate the application of the moving least-square approximation in the numerical simulation of the classical bending beam and high-order continuum bending beam. Some key issues such as the imposition of the boundary conditions are studied and discussed.

收稿日期：2014-03-10

基金項目：河南省科技攻關(guān)計劃項目(122102210492)

作者簡介：朱付保(1974-)，男，河南柘城人，副教授，博士，主要研究方向為智能信息處理，空間數(shù)據(jù)庫、地理信息系統(tǒng)、數(shù)據(jù)挖掘。

文章編號：1671-6906(2015)03-0085-05

對于傳統(tǒng)的連續(xù)梁來說，曲率是一個基本變量，它被近似為橫向位移的二階導(dǎo)數(shù)。在應(yīng)變梯度連續(xù)理論[1-2]、偶應(yīng)力理論[3]等高階連續(xù)理論中，需要考慮位移的二階導(dǎo)數(shù)。在微/納機(jī)電系統(tǒng)中，經(jīng)常把元器件處理為高階連續(xù)梁[4-6]，這時需要考慮撓度的三階甚至更高階的導(dǎo)數(shù)。高階連續(xù)問題給數(shù)值模擬帶來新的挑戰(zhàn)，利用有限元法構(gòu)造高階連續(xù)形函數(shù)需要做大量工作[7]。通過應(yīng)用適當(dāng)?shù)幕瘮?shù)和權(quán)函數(shù)，利用移動最小二乘法可以很容易地構(gòu)造出具有高階連續(xù)特征的形函數(shù)[8-9]，位移的高階導(dǎo)數(shù)可以直接用節(jié)點位移來近似，給數(shù)值離散帶來很大的方便。但是，它也帶來一些新的問題，比如，在數(shù)值計算中，高階應(yīng)變不再是顯式的變量。如何基于位移高階導(dǎo)數(shù)施加轉(zhuǎn)角約束，如何施加高階面力邊界條件，都是人們面臨的新的問題。本文以彎曲梁為例，研究移動最小二乘法在高階連續(xù)問題數(shù)值模擬中的應(yīng)用，并討論如何施加邊界條件等。

1移動最小二乘法

(1)

對于本文研究的彎曲梁問題，選用的一維三次基函數(shù)為：

(2)

(3)

式中：wI(x-xI)是權(quán)函數(shù)，當(dāng)x在xI影響域內(nèi)部時，wI(x)>0；當(dāng)x在xI影響域的邊界和外部時，wI(x)=0。

對(3)式求導(dǎo)：

=0,j=1,2，…,m

(4)

令

B(x)=[w1(x)b(x1)w2(x)b(x2)…wN(x)b(xN)]

(5)

移動最小二乘法的形函數(shù)[8-9]可以計算為：

(6)

令

(7)

形函數(shù)φ(x)的1～3階導(dǎo)數(shù)可計算為：

(8)

權(quán)函數(shù)w(x)在移動最小二乘近似中具有很重要的作用，它在節(jié)點xI處的值最大，并且具有緊支撐性，即當(dāng)r=‖x-xI‖/dmI>1時,w(r)=0,dmI為每個節(jié)點控制的支撐域直徑。本文采用三次樣條權(quán)函數(shù)[8-9]，即

(9)

圖1為移動最小二乘法形函數(shù)和它的1、2階導(dǎo)數(shù)的圖像。

(a)形函數(shù)φ

(b)一階導(dǎo)數(shù)φ ,x

(c)二階導(dǎo)數(shù) φ ,xx 圖1　移動最小二乘法形函數(shù)及其導(dǎo)數(shù)

2傳統(tǒng)彎曲梁的無網(wǎng)格法

圖2為一個簡支梁模型，在有限元方法中，經(jīng)常把轉(zhuǎn)角θ處理為節(jié)點變量，曲率由節(jié)點的撓度和轉(zhuǎn)角插值得出[7]。

圖2　簡支梁模型

考慮到移動最小二乘法形函數(shù)具有高階連續(xù)性，本文把梁的應(yīng)變能表示為：

(10)

撓度的二階導(dǎo)數(shù)直接用節(jié)點的撓度插值，即

(11)

其中：w,xx為撓度的二階導(dǎo)數(shù)；φI,xx為移動最小二乘法形函數(shù)的二階導(dǎo)數(shù)；wI(x)為節(jié)點撓度；n為計算點影響域內(nèi)的節(jié)點數(shù)。

由于轉(zhuǎn)動慣量I=∫Ay2dA，式(10)可以變?yōu)椋?/p>

(12)

彎曲梁剛度矩陣可計算為：

(13)

在計算剛度矩陣時，積分區(qū)間的布置與節(jié)點的布置分別獨立進(jìn)行，為了方便計算，本文中積分區(qū)間的布置與節(jié)點的布置一致(相鄰兩節(jié)點確定一個積分區(qū)間)，然后應(yīng)用高斯積分法對每個積分區(qū)間進(jìn)行積分。通過求解下列方程組可以得到節(jié)點的撓度:

KU=F

(14)

其中:U為節(jié)點撓度向量(不包含節(jié)點轉(zhuǎn)角)；F為節(jié)點力矢量。

(K+Kp)U=F+Fp

(15)

(16)

罰數(shù)α可被取為比剛度矩陣元素大5～6個量級的正數(shù)。

對于轉(zhuǎn)角邊界條件，將轉(zhuǎn)角用節(jié)點撓度和形函數(shù)的一階導(dǎo)數(shù)插值后，可由罰數(shù)法來施加。

用一個簡支梁的例子驗證方法的收斂性和精度。假設(shè)圖2中簡支梁的跨中受集中力F作用，梁截面高度h=1m，厚度b=0.01m，梁中點撓度δ可以采用經(jīng)典力學(xué)方法計算為FL3/48EI。為了驗證節(jié)點數(shù)量對結(jié)果的影響，每個積分單元布置3個高斯點，節(jié)點的影響域半徑定為3.0。表1所示為不同節(jié)點數(shù)時梁中點撓度值。

表1　不同節(jié)點數(shù)時梁中點撓度值　10 -4 m

由表1可以看出，所得結(jié)果收斂性很好。為了分析影響域大小對計算結(jié)果的影響，固定梁上節(jié)點數(shù)為61，節(jié)點的影響域半徑(R=Dmax(xj+1-xj))從2.0(xj+1-xj)變化到5.0(xj+1-xj)。

梁中點撓度值如表2所示。

表2　不同D max數(shù)值時的梁中點撓度值　10 -4 m

由表2可以看出,Dmax取3.0時可以得到最好的結(jié)果。

對圖3所示的兩個懸臂梁進(jìn)行計算。

(a)受集中力的懸臂梁 　　(b) 受角位移的懸臂梁 圖3　懸臂梁計算簡圖

梁在自由端分別受集中力和角位移作用，角位移用罰數(shù)法施加，梁的橫截面尺寸均為0.6m×0.6m，集中力P=100kN，角位移θ=0.1rad，梁的長度l=5m，彈性模量E=2.06×1011Pa，梁上布置的節(jié)點數(shù)固定為101。圖4為梁的撓度和轉(zhuǎn)角隨橫截面位置的變化曲線。

(a)受集中力作用的懸臂梁撓度曲線

(b) 受集中力作用的懸臂梁轉(zhuǎn)角曲線

(c)受角位移作用的懸臂梁撓度曲線

(d)受角位移作用的懸臂梁轉(zhuǎn)角曲線 圖4　彎曲梁的撓度曲線和轉(zhuǎn)角曲線

由圖4可以看出，本文得到的撓度和轉(zhuǎn)角與經(jīng)典力學(xué)的計算得到的結(jié)果相符較好，說明用罰數(shù)法施加位移邊界條件和轉(zhuǎn)角邊界條件均可取得很好的效果。

3高階連續(xù)梁的無網(wǎng)格法

在高階連續(xù)梁模型中，梁的應(yīng)變和高階應(yīng)變[4-5]分別定義為：

(17)

應(yīng)力和高階應(yīng)力定義為：

σxx=Eεxx+lxεxxx,τxxx=g2Eεxxx+lxEεxx,

τyxx=g2Eεyxx

(18)

式中:E為彈性模量;g和lx分別為高階連續(xù)梁的體本征尺寸影響因子和面本征尺寸影響因子。

應(yīng)變能可以表示為：

(19)

式中，A表示梁的橫截面。

對于高階連續(xù)梁，出現(xiàn)了撓度的三階導(dǎo)數(shù)，本文直接將其近似為：

(20)

高階連續(xù)梁的剛度矩陣計算為：

(21)

可以采用類似于本文上述方法建立無網(wǎng)格計算框架，數(shù)值計算直接給出節(jié)點的撓度值。

下面通過數(shù)值計算分析尺寸因子產(chǎn)生的影響。構(gòu)建一個簡支梁，其高度等于厚度(h=b)，梁上布置的節(jié)點總數(shù)固定為61，決定節(jié)點的影響域半徑固定為3.0，邊界條件的施加與本文上述方法相同。分析尺寸因子g的影響時，固定尺寸因子lx為0.001 25b，令g從0變化到5.0。從圖5可以看出，隨著g值的不斷增大。其對計算結(jié)果的影響不斷增大。圖5中縱坐標(biāo)為梁中點撓度相對于傳統(tǒng)力學(xué)結(jié)果(兩個尺寸因子均為0)的相對變化量。分析尺寸因子lx的影響時，將尺寸因子g固定為0.187 5b，令lx從0變化到0.8。從圖6可以看出，當(dāng)lx值在某一范圍內(nèi)時，對計算結(jié)果有較大的影響。圖6中縱坐標(biāo)同樣為梁中點撓度相對于傳統(tǒng)力學(xué)計算結(jié)果的相對變化量。

圖5　尺度影響因子g對結(jié)果的影響

圖6　尺度影響因子lx對結(jié)果的影響

4結(jié)語

高階連續(xù)問題中出現(xiàn)了位移的高階導(dǎo)數(shù)，給有限元法等傳統(tǒng)數(shù)值方法帶來了很大麻煩。本文應(yīng)用移動最小二乘法構(gòu)造具有高階連續(xù)特征的形函數(shù)，位移的高階導(dǎo)數(shù)直接由節(jié)點位移插值得到，可以方便地建立無網(wǎng)格計算框架，數(shù)值計算結(jié)果說明本文的方法有效。

在數(shù)值執(zhí)行中，節(jié)點位移直接通過求解方程組得到，任一位置的轉(zhuǎn)角可基于節(jié)點位移插值得到。由于轉(zhuǎn)角為位移的一階導(dǎo)數(shù)，轉(zhuǎn)角邊界條件可以較方便地通過罰數(shù)法來施加。對于高階連續(xù)問題中的高階面力等復(fù)雜邊界條件，由于可以表示為位移導(dǎo)數(shù)的函數(shù)，所以也可以用罰數(shù)法來施加，這將在后續(xù)的工作中進(jìn)行研究。

參考文獻(xiàn):

[1]XiaZC,HutchinsonJW.CrackTipFieldsinStrainGradientPlasticity[J].JournalofMechanicsPhysicsofSolids, 1996(44): 1621-1648.

[2]FleckRD,HutchinsonJW.StrainGradientPlasticity[J].AdvancedAppliedMechanics, 1997(33): 295-361.

[3]ToupinRA.ElasticMaterialswithCoupleStresses[J].Arch.Ration.Mech.Anal., 1962(11): 385-414.

[4]CivalekO,DemirC.StrainGradientElasticityandModifiedCoupleStressModelsforBucklingAnalysisofAxiallyLoadedMicro-scaledBeams[J].InternationalJournalofEngineeringScience, 2011(49): 1268-1280.

[5]CivalekO,DemirC.BucklingandBendingAnalysesofCantileverCarbonNanotubesUsingtheEuler-BernoulliBeamTheorybasedonNon-localContinuumModel[J].AsianJournalofCivilEngineering, 2011(12): 651-661.

[6]李志宏. 微納機(jī)電系統(tǒng) (MEMS/NEMS) 前沿[J].中國科學(xué), 2012, 42(12): 1599-1615.

[7]梁醒培, 王輝. 應(yīng)用有限元分析[M]. 北京：清華大學(xué)出版社，2003.

[8]BelytschkoT,KrongauzY.MeshlessMethod:anOverviewandRecentDevelopments[J].Comput.MethodsAppl.Mech.Engrg., 1996(139): 3-47.

[9]BelytschkoT,LuYY,GuL.Element-freeGalerkinMethods[J].InternationalJournalforNumericalMethodsinEngineering, 1994(37):229-256.

[10]SunYZ,LiewKM.BendingBucklingofSingle-walledCarbonNaotubes:HigherOrderGradientContinuumandMesh-freeMethod[J].ComputerMethodinAppliedMechanicsandEngineering, 2008(197): 3001-3013.

[11]SunYZ,LiewKM.ApplicationoftheHigh-orderCauchy-BornRuleinMesh-freeContinuumandMultiscaleSimulationofCarbonNanotubes[J].InternationalJournalforNumercalMethodsinEngineering, 2008, 75(10): 1238-1258.

(責(zé)任編輯：姜海芹)

The Application of the Moving Least-square Approximation in the

Numerical Simulation of the Higher-order Continuum Problems

CHEN Gen-sheng， YANG Zi-sheng, SUN Yu-zhou

(Zhongyuan University of Technology，Zhengzhou 450007, China)

Key words:bending beam; higher-order continuum; moving least-square approximation; mesh-free method; penalty method