杜成斌 李欣竹 田新冉 黃文倉

(河海大學(xué) 力學(xué)與材料學(xué)院，南京 211100)

在數(shù)值計算中，網(wǎng)格劃分的質(zhì)量是影響計算精度和計算效率的重要因素.常規(guī)有限元網(wǎng)格劃分無論是利用軟件還是人工剖分，網(wǎng)格均不允許存在懸掛節(jié)點(diǎn)，因此在材料分區(qū)和復(fù)雜邊界附近均需用較密的網(wǎng)格進(jìn)行過渡，既增加網(wǎng)格生成的工作量，也降低了結(jié)構(gòu)計算效率.

Matlab軟件有強(qiáng)大的計算性能，其內(nèi)置函數(shù)imread可以從圖像文件中讀取灰度或彩色圖像，進(jìn)而利用內(nèi)置函數(shù)qtdecomp對圖像進(jìn)行四叉樹分解，因此可以快速劃分網(wǎng)格，但是直接與有限元(finite element method，F(xiàn)EM)結(jié)合會產(chǎn)生懸掛節(jié)點(diǎn)[1].將其與新近提出的比例邊界有限元法(scaled finite element method，SBFEM)結(jié)合可克服有限元法對網(wǎng)格剖分不允許存在懸掛節(jié)點(diǎn)的問題[2].比例邊界有限元法是由Song和Wolf[3-5]提出的一種新型半解析數(shù)值計算方法，其應(yīng)力場和位移場在徑向是解析的，環(huán)向具有與有限元相同的精度，且僅需對結(jié)構(gòu)的邊界進(jìn)行離散，降低了一維計算維度.自提出之后，SBFEM得到快速發(fā)展.

Birk等[6]首次提出基于圖像四叉樹網(wǎng)格與SBFEM相結(jié)合的方法并模擬了地震波在不均勻土壤中的傳播；Ooi等[7]提出了基于混合多邊形和四叉樹網(wǎng)格相結(jié)合的SBFEM算法，高效地模擬了裂紋擴(kuò)展問題，且網(wǎng)格對復(fù)雜的邊界有更好的適應(yīng)性；隨后Albert Saputra等[8]基于圖像的四叉樹網(wǎng)格和八叉樹網(wǎng)格進(jìn)行了混凝土試件的靜力分析，研究中SBFEM的單元為規(guī)則的正方形或正六面體，取得了比較好的效果.鄒德高等[9]結(jié)合多邊形比例邊界有限元法和四叉樹網(wǎng)格對面板壩進(jìn)行跨尺度精細(xì)化靜動力分析，計算結(jié)果與FEM計算結(jié)果非常接近.He等[10]將圖像四叉樹與SBFEM結(jié)合擴(kuò)展到了非均質(zhì)材料的導(dǎo)電性分析，并驗證了該方法的合理性.以上研究中大多采用規(guī)則四叉樹對結(jié)構(gòu)進(jìn)行剖分，研究對象大都為規(guī)則結(jié)構(gòu)，如應(yīng)用在不規(guī)則結(jié)構(gòu)邊界，會產(chǎn)生鋸齒形狀，與實際幾何邊界不符.

基于圖像進(jìn)行四叉樹網(wǎng)格劃分與SBFEM結(jié)合可克服FEM對網(wǎng)格剖分不允許存在懸掛節(jié)點(diǎn)的問題.一方面充分利用Matlab先進(jìn)高效的圖像讀取功能，通過不同的像素對結(jié)構(gòu)的材料進(jìn)行自動分區(qū)，同時四叉樹算法非常方便快捷，這些優(yōu)勢克服了傳統(tǒng)的網(wǎng)格生成方法的缺陷，極大地提高了網(wǎng)格生成效率，使大家從繁瑣的網(wǎng)格剖分中解脫出來，更適合非均質(zhì)材料和復(fù)雜邊界的網(wǎng)格過渡.另一方面利用比例邊界有限元法半解析特點(diǎn)，可高效地進(jìn)行結(jié)構(gòu)動力分析以得到較高的計算精度.

本文采用一種基于圖像與四叉樹網(wǎng)格劃分相結(jié)合的方式以實現(xiàn)結(jié)構(gòu)的快速網(wǎng)格劃分，在邊界附近依實際幾何邊界對單元進(jìn)行切割，以滿足結(jié)構(gòu)的實際幾何邊界條件，改進(jìn)了以往大多單元模式僅為規(guī)則正方形、剖分后邊界為鋸齒形的缺陷.將其與比例邊界有限元法相結(jié)合，建立基于圖像的結(jié)構(gòu)動力計算，極大減輕了以往計算網(wǎng)格生成過程的繁瑣工作量，并通過兩個算例驗證了本文方法的先進(jìn)性和高效性.

1 基于圖像進(jìn)行四叉樹網(wǎng)格劃分

1.1 四叉樹數(shù)據(jù)結(jié)構(gòu)

四叉樹是一種樹狀數(shù)據(jù)結(jié)構(gòu)，每個節(jié)點(diǎn)最多有四個子樹.從四叉樹的根單元開始逐層劃分，每一個單元可以進(jìn)一步劃分為四個子單元，直到所有單元滿足劃分要求則停止.

四叉樹可以快速實現(xiàn)網(wǎng)格劃分，從1 m級到1 mm級僅需10層劃分(210=1 024)，且劃分網(wǎng)格以正方形為主，可自動進(jìn)行高質(zhì)量單元的離散.

1.2 基于圖像的四叉樹網(wǎng)格劃分步驟

1.2.1 圖像導(dǎo)入以及最大、最小單元尺寸和閾值的確定

可以使用Matlab內(nèi)置函數(shù)imread快速讀取灰度圖像并儲存為N×N的數(shù)字矩陣，要求圖像的像素值N必須為2的整數(shù)冪，才可以進(jìn)一步進(jìn)行四叉樹劃分.分辨率(Resolution)為結(jié)構(gòu)實際尺寸除以N，其物理意義為一個像素代表的實際尺寸大小.矩陣中數(shù)字的大小代表顏色深淺，取值范圍為0～255，0代表黑色，255代表白色，數(shù)值越大顏色越淺，因此可以用不同的顏色代表不同的材料.

Matlab內(nèi)置函數(shù)qtdecomp可以對圖像進(jìn)行四叉樹劃分，在劃分前需要確定最大、最小單元尺寸(MaxDim、MinDim)和閾值(QTthreshold).尤其在進(jìn)行時域分析時應(yīng)保證每個波長至少有10個節(jié)點(diǎn)才能識別波的傳播[11]，即一個波長內(nèi)至少包含9個兩節(jié)點(diǎn)線單元，由此確定最大單元尺寸.最大單元尺寸和最小單元尺寸必須是2的整數(shù)冪.網(wǎng)格大小要求最小的單元尺寸不能小于MinDim，最大的單元尺寸不能大于MaxDim，且要求相鄰兩個像素顏色強(qiáng)度的差小于QTthreshold，滿足則停止劃分網(wǎng)格，反之繼續(xù)劃分.

讀取一個8×8像素大小的灰度圖如圖1(a)所示，該圖像儲存為對應(yīng)的二維數(shù)字矩陣，如圖1(b)所示.

圖1 灰度圖像及對應(yīng)的圖像矩陣

1.2.2 平衡四叉樹網(wǎng)格剖分

對圖1結(jié)構(gòu)進(jìn)行四叉樹劃分，設(shè)置MaxDim=4、MinDim=1、QTthreshold=2，遵循2∶1的原則進(jìn)行四叉樹劃分，即相鄰單元尺寸滿足2∶1的平衡原則.

首先，將整個圖形區(qū)域劃分成4個大小相等的正方形單元，如圖2(a).檢查顏色強(qiáng)度的均勻性(即是否滿足閾值條件)；其次，進(jìn)一步劃分已經(jīng)劃分過的子域，如圖2(b).第3次劃分后滿足所有條件，則停止劃分，如圖2(c).

圖2 四叉樹劃分過程

經(jīng)過平衡四叉樹網(wǎng)格劃分可能產(chǎn)生的單元類型共有6種，如圖3所示，經(jīng)4次旋轉(zhuǎn)可得到全部24種單元類型，懸掛節(jié)點(diǎn)可以視為邊界上增加一個兩節(jié)點(diǎn)線單元，處理起來非常方便，因此SBFEM可允許存在懸掛節(jié)點(diǎn)，每個多邊形都可以作為SBFEM中的一個計算子域.

圖3 遵循2∶1劃分規(guī)則可能生成的多邊形單元

1.2.3 邊界單元的切割

為了適應(yīng)結(jié)構(gòu)實際幾何邊界條件，在邊界處需要進(jìn)行切割，可以將正方形單元依邊界條件切割成任意多邊形單元，具體形式如圖4所示.

圖4 經(jīng)過切割可能產(chǎn)生的單元類型

2 比例邊界有限元法基本理論

比例邊界有限元法與有限元法類似，首先需要將整個結(jié)構(gòu)的計算區(qū)域離散為一個或多個子域，然后再將所有子域的質(zhì)量矩陣、剛度矩陣按順序組合起來，最后由平衡方程求解出位移，根據(jù)幾何方程和物理方程求得應(yīng)變和應(yīng)力結(jié)果.在每個子域內(nèi)必須設(shè)置一個比例中心，比例中心的選取需要滿足對整個子域邊界可見的條件[12]，以圖5這個子域為例，將比例中心O放置在中心處，以O(shè)為原點(diǎn)建立局部比例邊界坐標(biāo)系(ξ，η)，ξ為徑向坐標(biāo)，η為環(huán)向坐標(biāo)沿著比例中心逆時針放置.對于子域內(nèi)任意一點(diǎn)，笛卡爾坐標(biāo)(x，y)可以用局部比例邊界坐標(biāo)(ξ，η)與形函數(shù)N(η)表示為

圖5 子域模型

x=ξxb(η)=ξN(η)xb

(1)

y=ξyb(η)=ξN(η)yb

(2)

式中：xb=[x1，x2]T，yb=[y1，y2]T為子域邊界上兩節(jié)點(diǎn)線單元的節(jié)點(diǎn)坐標(biāo)；N(η)=[N1(η)，N2(η)]為兩節(jié)點(diǎn)線單元的插值形函數(shù).

由彈性動力學(xué)的基本方程可以推導(dǎo)得到SBFEM位移控制方程[13]：

E2u(ξ)+ω2M0ξ2u(ξ)=0

(3)

其中：u(ξ)為位移向量；ω為頻率；E0、E1、E2、M0為系數(shù)矩陣，具體見文獻(xiàn)[13].由方程(3)可得到單元的位移徑向解析解.

單元的環(huán)向采用類似有限元的方法離散，引進(jìn)插值函數(shù)，最終可建立整體運(yùn)動方程：

(4)

3 算例分析

3.1 雙材料板動力算例

如圖6(a)所示，邊長為L=400 m的正方形雙材料板，材料參數(shù)分別為：E1=10 GPa，ν1=0.15，ρ1=2 800 kg/m3；E2=20 GPa，ν2=0.2，ρ2=2 600 kg/m3.邊界條件為：y=0處為固定端，在y=400 m處施加y方向相反方向的沖擊荷載，荷載隨時間變化如圖6(b)所示，最大值為400 kN.

圖6 雙材料板模型及其所受沖擊荷載

考慮該問題為平面應(yīng)力問題，讀取雙材料板的模型的圖像，像素值為512×512，材料1的顏色用98表示，材料2的顏色用198表示，如圖7(a)所示，基于該圖像進(jìn)行四叉樹網(wǎng)格劃分.一個像素的大小代表0.78 m，即Resolution=0.78 m，選取MaxDim=64，MinDim=16，對應(yīng)的真實尺寸分別為49.92 m，12.48 m，進(jìn)行動力分析時，兩種材料的網(wǎng)格疏密可以不相同，只需要在交界處加密，計算網(wǎng)格如圖7(b)所示，共279個節(jié)點(diǎn)，232個單元.ABAQUS有限元計算網(wǎng)格如圖7(c)所示，共289個節(jié)點(diǎn)，256個單元.采用Newmark時間積分求解動力方程，Newmark參數(shù)選擇β=0.25，γ=0.5.

圖7 雙材料板的圖像及其計算網(wǎng)格

基于上述方法對雙材料板進(jìn)行四叉樹網(wǎng)格劃分并采用SBFEM進(jìn)行時域分析，計算總時間為20 s，增量步長為0.02 s.選取坐標(biāo)為(200，300)的M點(diǎn)和(200，100)的N點(diǎn)，分別計算兩個點(diǎn)的y方向位移和y方向應(yīng)力，與有限元軟件ABAQUS計算結(jié)果進(jìn)行對比，如圖8所示.

圖8 M點(diǎn)、N點(diǎn)y方向位移及應(yīng)力結(jié)果

對于雙材料板的動力分析，為得到較高精度的結(jié)果，ABAQUS需要較密的計算網(wǎng)格，本文方法只需在兩種材料交界處進(jìn)行加密，可以實現(xiàn)粗細(xì)網(wǎng)格的過渡.本算例SBFEM僅需9.3 s，而ABAQUS需37 s，后者為前者的近4倍.結(jié)果顯示，兩種方法計算所得位移和應(yīng)力結(jié)果非常接近，驗證了該程序求解動力問題的正確性.

3.2 重力壩地震響應(yīng)分析

某混凝土重力壩剖面示意圖如圖9所示.壩高為110 m，壩頂寬15 m，壩底寬75 m，假定壩上游為直立坡.采用無質(zhì)量地基，地基范圍分別由壩底向下和上游、下游兩側(cè)延伸至兩倍壩高.壩體混凝土泊松比ν1=0.15，密度ρ1=2 643 kg/m3，彈性模量E1=31.027 GPa，假定采用Rayleigh阻尼反映混凝土的阻尼性質(zhì)，即C=a0M+a1K，每一階振型的阻尼比均取0.05，通過計算前兩階振型及自振頻率，求出系數(shù)a0=1.010 21，a1=0.002 08.壩基材料泊松比ν2=0.3，彈性模量E2=20 GPa，將壩基與壩體看作一個整體，計算坐標(biāo)系為：順河向為x，豎向為y，建立比例邊界有限元計算模型.計算中考慮壩體自重和地震荷載，重力加速度取g=9.81 m/s2，將地震力以慣性力的形式施加在壩體上.

圖9 混凝土重力壩剖面示意圖(單位：m)

輸入的地震波E1 Centro水平向加速度時程曲線如圖10所示，豎直方向地震波取水平向的2/3，地震作用總時長為30 s，計算時間增量步長選取為0.02 s，邊界條件為：壩基兩側(cè)水平方向位移固定，壩基底部固定邊界.采用Newmark法進(jìn)行動力計算，Newmark參數(shù)選擇β=0.25，γ=0.5.

圖10 水平向E1 Centro地震時程曲線

將重力壩具體尺寸的圖形導(dǎo)入程序中，生成像素值為1 024×1 024的灰度圖像如圖11(a)所示，如果有材料分區(qū)，只要不同的材料分區(qū)用不同顏色表示即可.一個像素的大小代表0.5 m，即Resolution=0.5 m，選取MaxDim=128，MinDim=1，對應(yīng)的真實尺寸分別為64 m，0.5 m，基于該圖像進(jìn)行四叉樹網(wǎng)格劃分如圖11(b)所示，為了更好地模擬重力壩的邊界，采用前述單元切割法對邊界處的四叉樹單元進(jìn)行切割形成任意多邊形單元.圖11(c)為ABAQUS有限元網(wǎng)格圖.

圖11 重力壩圖像及其計算網(wǎng)格

選取壩踵、下游壩面折坡面處、壩頂上游處(如圖9所示A、B、C)3個點(diǎn)作為關(guān)鍵點(diǎn)進(jìn)行時程觀測，輸出A、B、C3點(diǎn)的x方向位移時程曲線和A、B兩個點(diǎn)的x方向應(yīng)力時程曲線，并與有限元軟件ABAQUS計算所得結(jié)果進(jìn)行對比，3點(diǎn)的水平向位移時程曲線如圖12所示，A、B點(diǎn)的x方向應(yīng)力時程曲線如圖13所示.

圖12 A、B、C點(diǎn)x方向位移時程曲線

圖13 A、B點(diǎn)x方向應(yīng)力時程曲線

其中ABAQUS有限元的節(jié)點(diǎn)總數(shù)1 377，單元總數(shù)1 234，自由度數(shù)為2 754；SBFEM節(jié)點(diǎn)總數(shù)627，單元總數(shù)445，自由度數(shù)為1 254，這些數(shù)值均只有有限元的36%～46%.從計算時間來看，SBFEM只需201.1 s，而ABAQUS需462 s，后者是前者的2.3倍.由圖12、13可知，本文方法計算所得結(jié)果與ABAQUS計算所得結(jié)果非常接近，但相比有限元節(jié)省了大量的計算自由度.

4 結(jié) 論

基于圖像進(jìn)行四叉樹網(wǎng)格劃分，能夠?qū)崿F(xiàn)網(wǎng)格的快速劃分以及不同材料之間網(wǎng)格的粗細(xì)過渡，與比例邊界有限元法相結(jié)合，避免了有限元模型中的懸掛結(jié)點(diǎn)問題.本文將邊界切割引進(jìn)四叉樹生成的正方形單元，可以依邊界條件將正方形單元切割成任意多邊形單元，避免以往四叉樹單元生成的邊界為鋸齒形的缺陷.數(shù)值算例表明，與有限元模型相比，本文方法極大減小了結(jié)構(gòu)網(wǎng)格生成的前處理工作量，且粗細(xì)網(wǎng)格過渡方便快捷，模型的單元數(shù)、結(jié)點(diǎn)數(shù)目均大大減少，不僅有較高的計算精度，而且也有較高的計算效率，是一種非常有前景的新型結(jié)構(gòu)計算方法.