彈性應(yīng)力波的時(shí)空多尺度計(jì)算方法

彈性應(yīng)力波的時(shí)空多尺度計(jì)算方法

許多復(fù)雜的力學(xué)問(wèn)題，如裂縫邊緣的應(yīng)力場(chǎng)問(wèn)題和熱耦合結(jié)構(gòu)的局部破壞，很難使用單一的數(shù)值方法獲得令人滿(mǎn)意的模擬結(jié)果。發(fā)展多尺度計(jì)算,針對(duì)不同區(qū)域采用不同的數(shù)值計(jì)算方法能夠有效解決這類(lèi)問(wèn)題,現(xiàn)已成為當(dāng)今數(shù)值模擬研究的一個(gè)熱點(diǎn)。在固體結(jié)構(gòu)損傷和破壞的力學(xué)問(wèn)題中,充分考慮離散元法和有限元法的優(yōu)勢(shì)和弊端,即離散元法能夠方便地呈現(xiàn)局部細(xì)觀(guān)非均勻性,但需要大量計(jì)算時(shí)間;有限元法善于快速求解整體結(jié)構(gòu)宏觀(guān)彈性變形,但難以準(zhǔn)確處理局部結(jié)構(gòu)細(xì)微損傷、斷裂等力學(xué)現(xiàn)象,將這兩種方法耦合實(shí)現(xiàn)多尺度計(jì)算,分別應(yīng)用于宏觀(guān)彈性變形區(qū)和局部損傷破壞區(qū),能夠準(zhǔn)確、高效率地分析這些復(fù)雜的力學(xué)問(wèn)題。目前的研究成果主要集中在空間多尺度方面,很少進(jìn)行時(shí)間多尺度以及時(shí)間與空間兼顧的多尺度計(jì)算方法的研究。事實(shí)上,時(shí)間多尺度與空間多尺度對(duì)于提高多尺度計(jì)算效率具有同等重要作用。本文采用力學(xué)與動(dòng)力學(xué)參數(shù)傳遞以及權(quán)值分配法,建立了離散元與有限元耦合的時(shí)空多尺度計(jì)算方法,既能全面實(shí)現(xiàn)空間多尺度,同時(shí)又兼顧時(shí)間多尺度,有效地提高了多尺度數(shù)值計(jì)算的效率。同時(shí),作者將這一時(shí)空多尺度計(jì)算方法應(yīng)用于激光輻照下受拉鋁板破壞行為的數(shù)值模擬研究中,取得了較好的效果。1模型和方法的計(jì)算1.1有限元時(shí)間步長(zhǎng)的限制條件離散元法與有限元法具有相似形式的求解方程,都可以歸結(jié)為一般的牛頓運(yùn)動(dòng)定律形式,課題組曾根據(jù)這一理論基礎(chǔ),通過(guò)引入過(guò)渡層概念,對(duì)元與節(jié)點(diǎn)重疊的有限元和離散元耦合空間多尺度方法(以下簡(jiǎn)稱(chēng)空間多尺度法)進(jìn)行了研究。研究發(fā)現(xiàn),這種多尺度法只部分實(shí)現(xiàn)了空間多尺度,同時(shí)尚未實(shí)現(xiàn)時(shí)間多尺度計(jì)算,數(shù)值模擬仍需要較多計(jì)算時(shí)間,多尺度高效率優(yōu)勢(shì)尚未得到較好發(fā)揮。離散元法與有限元法的時(shí)間步長(zhǎng)直接影響時(shí)間多尺度計(jì)算的實(shí)現(xiàn)。離散元法的時(shí)間步長(zhǎng)存在幾何和物理兩方面的限制。幾何方面,在一個(gè)時(shí)間步長(zhǎng)內(nèi)不允許一個(gè)元與其鄰居元沿它們的中心軸線(xiàn)相互貫穿對(duì)方;物理方面,在一個(gè)時(shí)間步長(zhǎng)內(nèi)所產(chǎn)生的動(dòng)量轉(zhuǎn)移不應(yīng)超過(guò)整個(gè)碰撞過(guò)程中全部動(dòng)量的轉(zhuǎn)移。離散元法計(jì)算的臨界時(shí)間步長(zhǎng)應(yīng)當(dāng)取兩者中的最小值。另外,為了提高碰撞過(guò)程的時(shí)間分辨率,計(jì)算中應(yīng)該選取大于1的時(shí)間步數(shù),這樣離散元數(shù)值模擬的時(shí)間步長(zhǎng)限制為Δtd≤Δtdcr=1kmini,j{ΔtGdcr,ΔtPdcr}(1)Δtd≤Δtdcr=1kmini,j{ΔtdcrG,ΔtdcrΡ}(1)式中:Δtd為離散元時(shí)間步長(zhǎng);Δtdcr為離散元臨界時(shí)間步長(zhǎng);k為時(shí)間步數(shù),k>1;i和j為兩個(gè)相鄰元;ΔtGdcr為離散元幾何限制臨界時(shí)間步長(zhǎng);ΔtPdcr為離散元物理限制臨界時(shí)間步長(zhǎng)。由于離散元法采用顯式差分格式求解,為了耦合兩種方法實(shí)現(xiàn)多尺度計(jì)算,有限元法的計(jì)算采用顯式求解法中的中心差分法。為了保證有限元求解方程的穩(wěn)定性,時(shí)間步長(zhǎng)需要滿(mǎn)足條件Δtf≤Δtfcr=LminE(1?ν)ρ(1+ν)(1?2ν)√(2)Δtf≤Δtfcr=LminE(1-ν)ρ(1+ν)(1-2ν)(2)式中:Δtf為有限元時(shí)間步長(zhǎng);Δtfcr為有限元臨界時(shí)間步長(zhǎng);Lmin為模型中網(wǎng)格單元最小尺寸;ρ為材料密度;E為楊氏模量;ν為泊松比。從離散元與有限元時(shí)間步長(zhǎng)的限制條件可知,離散元時(shí)間步長(zhǎng)取決于元尺寸、鄰居元距離以及元性質(zhì),所以離散元時(shí)間步長(zhǎng)是很小的。有限元時(shí)間步長(zhǎng)不僅與結(jié)構(gòu)或材料的性質(zhì)有關(guān),而且與模型的網(wǎng)格單元最小尺寸成正比。因此,為了實(shí)現(xiàn)時(shí)間多尺度計(jì)算,離散元采用固定較小時(shí)間步長(zhǎng)的同時(shí),有限元需要采用較大的時(shí)間步長(zhǎng),這時(shí)根據(jù)式(2)就需要模型具有較大的有限元網(wǎng)格單元最小尺寸。1.2有限元時(shí)間步長(zhǎng)求解在空間多尺度法的基礎(chǔ)上,基于以上思想作者建立了既能有效實(shí)現(xiàn)時(shí)間多尺度,同時(shí)又能完善空間多尺度的離散元與有限元耦合的時(shí)空多尺度計(jì)算方法,這種方法的計(jì)算模型如圖1所示。從圖中可以看出,在耦合過(guò)渡區(qū)中,有限元網(wǎng)格單元的最小尺寸不必受到元的大小及元間距的限制,這樣網(wǎng)格單元的最小尺寸可以取得較大,從而增大有限元時(shí)間步長(zhǎng),而這時(shí)離散元可以采用較小的時(shí)間步長(zhǎng),完成時(shí)間多尺度計(jì)算。另外,該模型中整體有限元網(wǎng)格的數(shù)量受耦合過(guò)渡區(qū)網(wǎng)格尺度的影響較小,具備完善的空間多尺度計(jì)算能力。在計(jì)算模型的耦合過(guò)渡區(qū)中,離散元和有限元相互為對(duì)方提供邊界條件,并采用權(quán)值分配法進(jìn)行力學(xué)和動(dòng)力學(xué)參數(shù)的傳遞。也就是說(shuō),離散元為有限元提供力邊界條件采用式(3)進(jìn)行計(jì)算,有限元為離散元提供的速度邊界條件采用式(4)進(jìn)行計(jì)算。當(dāng)耦合過(guò)渡區(qū)內(nèi)施加外載荷時(shí),需要將外載荷加到式(3)和式(4)的網(wǎng)格節(jié)點(diǎn)或者元的相應(yīng)變量中。F(D)FE=∑i=1nk(m)iFi(m)DE+∑p∑l<m,i=1qk(l)iFi(l)DE(3)FFE(D)=∑i=1nki(m)FDEi(m)+∑p∑l<m,i=1qki(l)FDEi(l)(3)式中:F(D)FEFE(D)為離散元提供力邊界條件后的網(wǎng)格節(jié)點(diǎn)作用力;k(m)ii(m)為第m網(wǎng)格內(nèi)第i元在該網(wǎng)格節(jié)點(diǎn)的權(quán)值;Fi(m)DEDEi(m)為第m網(wǎng)格內(nèi)第i元受到的合力;n表示第m網(wǎng)格包含的元數(shù)目;k(l)ii(l)為第l(l<m)網(wǎng)格內(nèi)第i元在該網(wǎng)格節(jié)點(diǎn)的權(quán)值;Fi(l)DEDEi(l)為第l(l<m)網(wǎng)格內(nèi)第i元受到的合力;q為第l網(wǎng)格包含的元數(shù)目;p為包含該網(wǎng)格節(jié)點(diǎn)的所有l(wèi)(l<m)網(wǎng)格的數(shù)目。V(F)DE=∑i=13kiViFE(4)VDE(F)=∑i=13kiVFEi(4)式中:V(F)DEDE(F)為有限元提供速度邊界條件后的元速度;ki為該元所屬網(wǎng)格的第i節(jié)點(diǎn)(三角形網(wǎng)格,i=1,2,3)在該元位置處的權(quán)值;ViFEFEi為該元所屬網(wǎng)格單元的第i節(jié)點(diǎn)的速度。為了實(shí)現(xiàn)耦合過(guò)渡區(qū)材料密度的匹配性,需要使過(guò)渡區(qū)內(nèi)網(wǎng)格節(jié)點(diǎn)的質(zhì)量與相應(yīng)元過(guò)渡到該節(jié)點(diǎn)的質(zhì)量相等,計(jì)算方法為M(D)FE=∑i=1nk(m)iMi(m)DE+∑p∑l<m,i=1qk(l)iMi(l)DE(5)ΜF(xiàn)E(D)=∑i=1nki(m)ΜDEi(m)+∑p∑l<m,i=1qki(l)ΜDEi(l)(5)式中:M(D)FEFE(D)為元質(zhì)量傳遞后的網(wǎng)格節(jié)點(diǎn)質(zhì)量;Mi(m)DEDEi(m)為第m網(wǎng)格內(nèi)第i元的質(zhì)量;Mi(l)DEDEi(l)為第l(l<m)網(wǎng)格內(nèi)第i元的質(zhì)量。在時(shí)空多尺度數(shù)值計(jì)算中,離散元時(shí)間步長(zhǎng)設(shè)為dt,有限元時(shí)間步長(zhǎng)設(shè)為fdt,兩者關(guān)系為fdt/dt=N。也就是說(shuō),在每個(gè)有限元時(shí)間步長(zhǎng)求解計(jì)算中,執(zhí)行N次離散元計(jì)算過(guò)程來(lái)完成時(shí)間多尺度計(jì)算。整個(gè)系統(tǒng)的計(jì)算流程如圖2所示。2有限元數(shù)值模擬結(jié)果用彈性應(yīng)力波傳播算例來(lái)驗(yàn)證時(shí)空多尺度計(jì)算方法的有效性和優(yōu)越性。選用一根長(zhǎng)度、寬度和厚度分別為3.2×10-1m、2.4×10-3m和3.2×10-4m的細(xì)長(zhǎng)桿。桿的右端固定,左端施加一個(gè)突然加載和卸載的分布載荷,載荷的大小為1.0×104N/m,脈寬為1.5×10-5s。桿的材料為鋁,密度為2.7×103kg/m3,彈性模量為7.0×1015Pa。在時(shí)空多尺度數(shù)值計(jì)算中,將桿中部4.0×10-3m的區(qū)域用離散元模擬,而其余部分用有限元模擬。在細(xì)長(zhǎng)桿上取3個(gè)不同位置(A、B、C)來(lái)觀(guān)察應(yīng)力波的傳播情況,這3個(gè)位置分別距離桿左端7.80×10-2m,1.60×10-1m和2.42×10-1m,如圖3所示。圖4為使用單一有限元法、時(shí)空多尺度法以及解析法時(shí)應(yīng)力波在細(xì)長(zhǎng)桿A、B、C三個(gè)位置的傳播情況。從圖中可以看出,時(shí)空多尺度法與有限元法相比,兩者在三個(gè)位置獲得的應(yīng)力波傳播情況基本一致,應(yīng)力波能夠在有限元區(qū)(A和C位置)和離散元區(qū)(B位置)之間光滑過(guò)渡。同時(shí),時(shí)空多尺度法計(jì)算的應(yīng)力波到達(dá)三個(gè)位置的時(shí)間及峰值與解析解基本吻合。以上結(jié)果驗(yàn)證了離散元與有限元耦合的時(shí)空多尺度計(jì)算方法的有效性。由于引入了阻尼,波在傳播過(guò)程中有一定的衰減。由于桿的橫向慣性效應(yīng),波形將出現(xiàn)劇烈振蕩,圖4中的模擬結(jié)果很好地說(shuō)明了這一點(diǎn)。為了研究時(shí)空多尺度方法的計(jì)算效率,將空間多尺度法與耦合過(guò)渡區(qū)內(nèi)布設(shè)不同層離散元的時(shí)空多尺度法在彈性應(yīng)力波傳播算例中的計(jì)算時(shí)間進(jìn)行了對(duì)比,如表1所示。根據(jù)式(2)可知,有限元臨界時(shí)間步長(zhǎng)與網(wǎng)格單元最小尺寸成正比。為了對(duì)比分析方便,在時(shí)空多尺度法中,耦合過(guò)渡區(qū)網(wǎng)格尺寸(不大于非耦合過(guò)渡區(qū)網(wǎng)格最小尺寸)為離散單元間距的倍數(shù),令耦合過(guò)渡區(qū)內(nèi)網(wǎng)格尺寸與元間距相等時(shí),有限元與離散元的時(shí)間步長(zhǎng)相等,隨著元層數(shù)增加,有限元時(shí)間步長(zhǎng)相應(yīng)增倍。從表1中可以看出,當(dāng)采用時(shí)空多尺度法時(shí),隨著耦合過(guò)渡區(qū)內(nèi)離散元層數(shù)的增多以及有限元時(shí)間步長(zhǎng)的增大,計(jì)算時(shí)間明顯減小。并且,時(shí)空多尺度法所有情況下的計(jì)算時(shí)間均明顯少于空間多尺度法的計(jì)算時(shí)間,當(dāng)耦合過(guò)渡區(qū)內(nèi)布置五層離散元時(shí),采用時(shí)空多尺度法的計(jì)算時(shí)間僅相當(dāng)于空間多尺度法計(jì)算時(shí)間的1/3?？梢?jiàn),時(shí)空多尺度法在節(jié)省計(jì)算時(shí)間方面具有較大的優(yōu)越性,有效地發(fā)揮了多尺度數(shù)值計(jì)算的高效率優(yōu)勢(shì)。3受拉鋁板微裂紋的有限元模擬研究激光輻照下預(yù)載金屬結(jié)構(gòu)的破壞機(jī)理對(duì)于工業(yè)生產(chǎn)和軍事應(yīng)用具有重要的科學(xué)意義。本文使用離散元與有限元耦合的時(shí)空多尺度計(jì)算方法對(duì)激光輻照下受拉鋁板破壞行為進(jìn)行數(shù)值模擬。選用長(zhǎng)度、寬度和厚度分別為1.0×10-1、5.92×10-2和6.5×10-4m的矩形薄板,板的材料為鋁合金,密度為2.7×103kg/m3,彈性模量為7.0×1015Pa,熔點(diǎn)為993K,熱導(dǎo)率為233W/(m·K),熱膨脹系數(shù)為2.9×10-7K-1。薄板沿y軸受到雙向均勻預(yù)拉載荷,載荷力的大小為該板靜態(tài)破壞載荷Fb(大小為1.4×105N/m)的1/4,同時(shí)薄板中部受到高斯分布的激光輻照,光斑半徑為1.5×10-2m,激光束平均功率密度為2.5×106W/m2,輻照表面吸收系數(shù)為0.15。由于薄板中部受到熱力耦合的作用,將發(fā)生損傷、斷裂等非均勻破壞行為,采用離散元進(jìn)行模擬,共有5289個(gè)離散元;而遠(yuǎn)離激光輻照的區(qū)域,其應(yīng)力狀態(tài)仍在彈性范圍,采用有限元進(jìn)行模擬,共有294個(gè)有限元網(wǎng)格,計(jì)算模型如圖5所示。圖6為激光輻照下受拉鋁板破壞過(guò)程的時(shí)空多尺度數(shù)值模擬圖,圖中的黑色三角形標(biāo)出微裂紋的位置。從圖中可以看出,微裂紋首先出現(xiàn)在垂直于拉伸方向的光斑直徑與光斑邊緣的兩個(gè)交點(diǎn)附近,隨后以該微裂紋處作為裂源,在其周?chē)壬罅啃挛⒘鸭y,并沿橫向迅速向鋁板邊緣和光斑中心雙向擴(kuò)展。微裂紋的擴(kuò)大將形成宏觀(guān)裂紋,并在很短的時(shí)間內(nèi)沿水平方向貫穿整個(gè)鋁板,直至完全斷裂,這些模擬結(jié)果與文獻(xiàn)提供的實(shí)驗(yàn)結(jié)果基本一致。4多尺度計(jì)算方法數(shù)值模擬本文建立了離散元與有限元耦合的時(shí)空多尺度計(jì)算方法,通過(guò)彈性應(yīng)力波傳播算