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

陳根生， 楊子勝， 孫玉周

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

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

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

中圖分類(lèi)號(hào)：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

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

作者簡(jiǎn)介：朱付保(1974-)，男，河南柘城人，副教授，博士，主要研究方向?yàn)橹悄苄畔⑻幚?，空間數(shù)據(jù)庫(kù)、地理信息系統(tǒng)、數(shù)據(jù)挖掘。

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

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

1移動(dòng)最小二乘法

(1)

對(duì)于本文研究的彎曲梁?jiǎn)栴}，選用的一維三次基函數(shù)為：

(2)

(3)

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

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

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

(4)

令

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

(5)

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

(6)

令

(7)

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

(8)

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

(9)

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

(a)形函數(shù)φ

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

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

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

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

圖2　簡(jiǎn)支梁模型

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

(10)

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

(11)

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

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

(12)

彎曲梁剛度矩陣可計(jì)算為：

(13)

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

KU=F

(14)

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

(K+Kp)U=F+Fp

(15)

(16)

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

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

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

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

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

梁中點(diǎn)撓度值如表2所示。

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

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

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

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

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

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

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

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

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

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

3高階連續(xù)梁的無(wú)網(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表示梁的橫截面。

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

(20)

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

(21)

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

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

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

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

4結(jié)語(yǔ)

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

在數(shù)值執(zhí)行中，節(jié)點(diǎn)位移直接通過(guò)求解方程組得到，任一位置的轉(zhuǎn)角可基于節(jié)點(diǎn)位移插值得到。由于轉(zhuǎn)角為位移的一階導(dǎo)數(shù)，轉(zhuǎn)角邊界條件可以較方便地通過(guò)罰數(shù)法來(lái)施加。對(duì)于高階連續(xù)問(wèn)題中的高階面力等復(fù)雜邊界條件，由于可以表示為位移導(dǎo)數(shù)的函數(shù)，所以也可以用罰數(shù)法來(lái)施加，這將在后續(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].中國(guó)科學(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