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

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

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