莊 茁, 嚴子銘, 姚凱麗, 崔一南, 柳占立

(清華大學 航天航空學院, 北京 100084)

1 引 言

具備跨尺度Spanning the scale力學分析和計算能力,提升從微納米到宏細觀尺度的計算能力,建立從離散到連續(xù)的材料本構模型,解決力學行為分析和結構設計等科學和工程問題,是力學學者長期為之奮斗的目標。盡管針對研究對象和具體問題,研究者提出了許多可行的研究方法和技術路線,但是仍缺少強勁的普適性解決方案。在跨尺度計算中,目前存在的困難是前提條件是否滿足物理機制,計算方法是否可行可靠,變量之間轉換是否合適合理,只有突破上述瓶頸,才能實施跨尺度計算分析。本文討論了固體力學領域的若干理論分析和計算問題,提出了幾種不同類型的跨尺度計算方法。同時,也考慮到在機器學習技術發(fā)展的新態(tài)勢下,以此拋磚引玉,探討解決更多的固體力學科學和結構工程設計難題。

2 建立位錯動力學與有限元耦合的計算方法,提出跨納米到微米尺度有限變形晶體塑性流動應力

晶體材料在現(xiàn)代工業(yè)中扮演了重要角色,隨著微納米力學和材料科學的不斷發(fā)展,對晶體材料的研究重點已經(jīng)從宏觀轉向細觀及納觀尺度。1994年,針對幾十至幾百微米的細銅絲扭轉和薄梁彎曲產(chǎn)生的越細越硬的實驗結果,研究者提出了應變梯度塑性理論[1,2],基于幾何必需位錯機制建立流動應力理論公式,揭示了強度的尺寸效應,但該理論只適用于10微米至毫米尺度的晶體塑性問題。

通常將幾百納米到10微米的尺度范圍稱作亞微米尺度。實驗研究發(fā)現(xiàn),亞微米尺度晶體塑性流動具有許多新的特性[3]。以微納米晶體柱單軸壓縮實驗為例,如圖1所示,在沒有應變梯度的前提下,單晶柱的屈服應力隨試樣尺寸的減小和溫度的升高而顯著增加,不但產(chǎn)生強度和變形的尺寸效應,而且還產(chǎn)生反常規(guī)的本構行為。一些研究者通過納米尺度的分子動力學MD模擬得到了一些有意義的結果,由于MD計算能力的限制難以揭示強度和變形的力學機制,因此,需要建立適用于亞微米尺度的理論和計算方法,作為向上尺度聯(lián)系應變梯度塑性理論和宏觀尺度經(jīng)典塑性力學,向下尺度聯(lián)系分子動力學模擬,跨越多個尺度揭示材料的強度和變形的尺寸效應,填補晶體塑性理論的空白。

經(jīng)典位錯動力學DDD計算模型是建立在晶體內部位錯運動的基礎上,其計算結果得到的晶體材料不是有限變形體,而真實的實驗材料是有限變形體,如圖1所示。晶體材料中含有高達1013/m2的位錯,位錯的伯格斯矢量只有約0.25 nm長度。位錯靜止是彈性問題,只有位錯運動才產(chǎn)生塑性應

圖1 微納米晶體柱單軸壓縮實驗,在無應變梯度條件下發(fā)生越細越硬的強度尺寸效應[3]

變εp。于是提出如何描述位錯的運動,建立位錯機制的晶體塑性有限變形理論的問題。我們提出了三維的位錯動力學與有限元耦合的計算模型DDD-FEM,首先將每一增量步的塑性應變代入有限元的單元積分點處,乘以切線剛度C得到塑性應力σp,然后將塑性應力返回到DDD模型,形成驅動位錯運動的力。重復這個耦合計算,直至完成計算過程,如圖2所示[4]。該計算模型突破了基于位錯機制計算有限變形晶體塑性的技術瓶頸,建立了微納尺度(200 nm～10 μm)晶體塑性流動應力理論,結合試驗數(shù)據(jù)驗證了無應變梯度下強度和變形的尺寸效應[5]。

圖2 位錯動力學與有限元耦合計算模型(DDD-FEM)

應用該計算模型,根據(jù)位錯源機制預測流動應力的演化。研究表明,激活單臂源的臨界分切應力由三部分組成,即晶格摩擦力以及與位錯密度ρ有關的彈性相互作用力和線張力,其中,線張力與平均有效位錯源長度λ有關,即位錯釘扎點至自由邊界或晶體界面的距離。基于位錯機制,我們提出了臨界分切應力τ的表達式為[6]

(1)

式中b為伯格斯矢量,μ為切線模量,k為系數(shù)。式(1)等號右端第一項為初始流動應力,第二項為基于應變梯度理論的流動應力,第三項為亞微米尺度的單晶體塑性流動應力。如果考慮多晶體材料,可以加上一項考慮晶界對流動應力影響的背應力[7]。圖3給出實驗中的分切應力數(shù)據(jù)是[111]取向單晶鎳柱在單軸壓縮下應變值達到3%時的結果。為方便對比,式(1)理論預測和離散位錯動力學模擬得到的穩(wěn)定流動應力均轉化成分切應力τ。可以看出,實驗、計算和理論得到的結果吻合良好,單晶柱直徑d的尺度跨越200 nm～1000 nm,反映了越細越硬的強度尺寸效應。Ryu等[8]肯定了該理論在預測流動應力尺寸效應的貢獻,實驗驗證了我們提出的穩(wěn)定單臂位錯源長度正比于試樣尺寸以及導致應變硬化消失的結論。

圖3 微納米單晶鎳柱單軸壓縮試驗數(shù)據(jù)與理論和位錯動力學計算結果對比[6]

塑性流動應力的理論還成功預測了單晶微柱變形模式的尺寸依賴性,揭示了微柱變形模式轉變的力學機理,得到了剪切帶形成的臨界轉變尺寸。隨著直徑尺寸減小,由30 μm至1 μm,變形由鼓形至S形,如圖4所示[5]。即在微米到納米轉變的亞微米尺度下,FCC晶體從雙剪切帶對稱滑移變形(鼓形)轉變?yōu)閱渭羟袔Х菍ΨQ滑移變形(S形),說明在10 μm以下的小直徑尺度下,微柱已經(jīng)沒有足夠的能量啟動第二個剪切帶滑移系,反映了亞微米尺度變形的尺寸效應,這是在經(jīng)典塑性理論和應變梯度塑性理論中無法揭示的變形機制。

圖4 單晶柱直徑由30 μm至1 μm,變形由桶形至S形,展示了變形模式的尺度依賴性[5]

3 建立粗粒化分子動力學模型,獲得準連續(xù)介質下聚脲材料的存儲模量和損耗模量

在高分子材料中,由微相結構可控的兩種材料組成的聚脲屬于嵌段共聚物,其軟相為氨基聚醚,硬相為異氰酸酯。由于軟相和硬相的團簇,產(chǎn)生各向異性的力學性質以及儲能和耗能機制。利用聚脲優(yōu)異的力學性質,可以設計防護爆炸沖擊波的輕質材料結構。通過全原子模型All-atom MD和粗?；Ｐ虲G-MD[9-11],計算和設計軟相/硬相微結構,如圖5所示。圖中的全原子模型邊長為4 nm,粗?；Ｐ偷倪呴L為40 nm,在相同勢函數(shù)和分子鍵力的前提下,粗?；Ｐ吞岣吡?個時間尺度量級的計算效率;相對于代表性體積單元,粗粒化模型提高了3個空間尺度量級,因此大大提高了計算分析的效率。由此研制出模量、黏彈性和硬度在較寬頻段和應變率范圍可控的聚脲材料,實現(xiàn)調控衰減沖擊波能量密度的設計目標[12,13]。

通過圖5所示的分子動力學粗?；Ｐ虲G-MD計算,得到維里應力(能量密度)為

(2)

式中m為原子質量,v為運動速度,F為原子之間的鍵力。采用Green-Kubo方法,得到時域下的剪切松弛函數(shù)為

G(t)=V/(kBT)〈ααβ(t)ααβ(0)〉

(3)

由此定義了剪切復模量,通過對式(3)進行傅里葉變換,轉換為頻域下的存儲模量,式中的ω為圓頻率,

(4)

和損耗模量

(5)

由此建立了損耗因子

tanδ=G″(ω)/G′(ω)

(6)

式(6)反映了聚合物材料的力學行為,若分母遠大于分子,反映聚合物趨近于固體性質;反之,說明聚合物趨近于流體性質。動態(tài)模量隨頻率變化,計算結果與超聲試驗和DMA試驗數(shù)據(jù)對比,如圖6所示[12,13]。

圖5 通過全原子模型和粗粒化模型計算和設計聚脲材料[12]

圖6 頻域下的存儲模量和頻域下的損耗模量

該結果形成細觀尺度材料模型的初始模量值[14],可以代入連續(xù)介質力學作為材料的初始彈性模量部分,由此跨越了納米尺度到微米尺度,應用于計算一類聚合物材料的能量衰減與耗散性能。這種粗粒化模型能夠跨越尺度的原則是放大了空間尺度和放慢了時間尺度,采用的方法是通過CG-MD計算的數(shù)值獲得了維里應力,通過維里應力計算時域下剪切松弛函數(shù),其也是跨尺度的變量。通過剪切松弛函數(shù)的傅里葉變換獲得頻域下的存儲和損耗模量。在保持原問題物理狀態(tài)不變的前提條件下,提高了計算效率。

4 建立數(shù)據(jù)驅動微觀CT與普通CT影像方法,獲得松質骨組織的等效模量和設計人工假體骨點陣結構

在骨科骨缺損的重建治療中,通常采用均勻設計或簡單梯度設計的人工假體骨,其力學性質與患者骨缺損部位的實際力學性質差異較大,無法實現(xiàn)假體骨匹配缺損部位力學性質的設計目標。因此,骨科骨缺損領域一方面需要解決臨床CT影像下對骨骼力學性質精準預測的問題,另一方面需要考慮如何設計具有良好力學性質匹配度的假體,形成正向的骨傳導信號,促進新骨組織的生長與發(fā)育。

近年來,圍關節(jié)骨缺損重建成為力學和醫(yī)學領域的研究熱點,同時也面臨諸多挑戰(zhàn)性難題,一是松質骨非均質性和各向異性規(guī)律復雜,生長方向性明顯,難以建立合適的力學理論模型、計算方法和實驗技術,實現(xiàn)定量描述缺損區(qū)的骨結構及其功能;二是骨缺損患者的病歷信息主要來自普通CT(Computed Tomography)影像,分辨率較低,難以辨識骨微觀組織結構的拓撲分布和力學性能;三是未進行個性化設計的骨假體與缺損部位組織之間難以形成應力重分配、變形協(xié)調和阻抗匹配的良性機制,使骨應力傳導達到刺激骨生成的水平。這些難題構成了骨缺損臨床治療的瓶頸。

隨著精準醫(yī)療概念在骨缺損治療領域的不斷發(fā)展,以及力學與醫(yī)學領域的交叉融合,本文通過數(shù)據(jù)驅動辨識Micro-CT高分辨率(0.06 mm)和臨床-CT低分辨率(0.6 mm)影像的特征值,建立了跨尺度的數(shù)據(jù)驅動方法,如圖7所示。應用數(shù)據(jù)驅動方法,一方面通過醫(yī)學影像獲取松質骨的非均勻密度和微結構取向分布等與松質骨材料力學性能強相關的形態(tài)學特征;另一方面可建立無監(jiān)督聚類模型,篩選與松質骨各向異性行為強相關的力學特征,發(fā)展骨組織各向異性本構模型。通過基于高分辨率Micro-CT影像建立高精度有限元模型,計算基于代表性體元RVE的各向異性彈性矩陣標簽數(shù)據(jù)庫。將數(shù)據(jù)驅動方法提供關鍵松質骨微結構形態(tài)學特征,與建立的骨組織非均質性、各向異性以及生長方向性的等效本構模型相融合,由此形成數(shù)據(jù)驅動與力學模型融合的多尺度建模方法,如圖7所示。

本文從獲取體積胞元的方式展開,首先,基于動物骨組織的高精度Micro-CT影像進行數(shù)字重建,獲取基于高分辨率影像的松質骨體積胞元模型,分別通過圖像分析方法和有限元計算獲取不同胞元對應的形態(tài)學參數(shù)(結構張量M)和各向異性力學性質(彈性矩陣C),并建立三維自編碼器獲取胞元對應的微結構特征向量fHR(Step1),以及fHR與對應胞元的結構張量M和彈性矩陣C的映射關系。進而,基于動物骨組織的低分辨率臨床CT影像建立以宏觀骨密度表示的松質骨體積胞元模型,并通過上采樣與神經(jīng)網(wǎng)絡方法,將對應的fHR作為標簽對模型進行訓練,實現(xiàn)僅通過臨床CT影像便可獲取微結構特征向量fHR(Step2),進而對結構張量M和彈性矩陣C的映射關系進行預測(Step3)。最終,基于獲得的圍關節(jié)松質骨小梁的等效模量和結構張量,可以進行松質骨組織增材制造的點陣設計[15],實現(xiàn)了個性化骨缺損的重建[16-18]。

圖7 從Micro-CT影像中獲取松質骨單元模型的形態(tài)學參數(shù)和微結構特征規(guī)律,并通過數(shù)據(jù)驅動普通CT影像的微結構特征規(guī)律,獲取等效各向異性力學性質

松質骨的力學性質基本滿足正交各向異性力學行為,具有9個獨立的彈性張量參數(shù)。為了從物理上證明本文方法的預測精度,選擇了松質骨立方體模型,比較了基于Micro-CT圖像的有限元方法和基于低分辨率臨床CT的預測模型給出的彈性矩陣C和結構張量M的結果。本文彈性矩陣C通過Voigt標記簡化為

(7)

結構張量M表示為

(8)

從圖8可以看出,預測的三組立方體體積分數(shù)為0.29,0.35和0.47。棕色立方體模型使用0.06 mm/pixel空間分辨率的Micro-CT影像,以骨小梁微結構形態(tài)表征。相應的深灰色立方體提取0.6 mm/pixel低分辨率的臨床CT影像?？梢?兩者數(shù)值基本吻合,說明基于影像的數(shù)據(jù)驅動方法轉換材料力學特征的正確性。

圖8 Micro-CT有限元計算結果與臨床低分辨率CT預測結果的對比

上述研究建立了臨床CT影像下精準預測松質骨等效力學參數(shù)的數(shù)據(jù)驅動模型,將結構張量作為特征量,依次關聯(lián)微米級的高分辨率micro-CT影像下的松質骨微結構形態(tài)特征規(guī)律和低分辨率臨床CT下的宏觀骨密度分布,實現(xiàn)了松質骨結構分析從微米尺度到毫米尺度的跨越,最終僅依靠臨床低分辨率CT影像便可對結構張量精準預測,進而對松質骨非均質各向異性力學性能進行表征,解決了骨科臨床中難以表征骨組織結構非均質性與各向異性的難題。

5 結 論

本文展示了在固體力學領域跨尺度計算的若干問題和研究進展,提出了三種不同的跨尺度計算方法。一是提出了位錯動力學與有限元耦合模型和計算方法,突破了基于位錯機制計算分析有限變形晶體塑性理論的瓶頸,建立了微納尺度(200 nm～10 μm)晶體塑性流動應力公式,結合試驗數(shù)據(jù)驗證在無應變梯度下強度和變形的尺寸效應,填補了小尺度晶體塑性流體應力理論的空白。二是應用納米尺度分子動力學模型和粗?；椒?計算獲得聚脲材料在時域和頻域下的存儲模量和損耗模量,通過動態(tài)加載分析試驗和超聲波試驗的數(shù)據(jù)驗證,解決了連續(xù)介質尺度下微相分離高分子共聚物的設計難題。三是建立了數(shù)據(jù)驅動變焦距方法,辨識高分辨率的微米尺度CT影像和臨床低分辨率的毫米尺度CT影像的特征值,建立了圍關節(jié)松質骨小梁的等效模量和結構張量,為非均質和各向異性骨組織結構增材制造點陣設計和個性化骨缺損重建,奠定了力學理論和計算的基礎。

以上針對三類不同科學和工程問題,采用了三種跨尺度計算方法,其共同特點是分別通過關聯(lián)兩個尺度下具有相同物理描述且不受實際計算所選取的時間尺度和空間尺度影響的物理量,依次為塑性應變、存儲與損耗模量、結構張量,實現(xiàn)了跨尺度計算。因此,根據(jù)問題導向,如何尋找并確定不同尺度計算方法下與時間尺度和空間尺度無關的關鍵特征物理量,將其作為橋梁,是解決跨尺度計算問題的一個關鍵因素。

后 記

鐘萬勰先生是中國計算力學科學領域的開拓者和領導者,是教書育人的楷模。第一次見到鐘先生是1990年在北京友誼賓館的國際會議上,鐘先生擔任大會主席,神采奕奕,思維敏捷,令我十分崇拜。在國內外的力學學術會議上,多次聆聽鐘先生充滿激情的學術報告,學習先生獨立創(chuàng)新的研究精神。也多次邀請鐘先生來清華講學,每次都受益匪淺。謹以此文紀念鐘萬勰院士九十誕辰。

2007年,在北京的全國力學大會與鐘先生合照