組合雜交三角形單元的加權(quán)能量正交關(guān)系

張玲，張偉偉

(西北工業(yè)大學(xué) 理學(xué)院，西安 710129)

組合雜交三角形單元的加權(quán)能量正交關(guān)系

張玲，張偉偉

(西北工業(yè)大學(xué) 理學(xué)院，西安 710129)

組合變分原理可以增強雜交元方法解的穩(wěn)定性。建立熱傳導(dǎo)方程基于區(qū)域分解的組合雜交有限元方法，給出單元上溫度梯度插值為線性、但溫度插值為協(xié)調(diào)線性插值與非協(xié)調(diào)二次插值之和的組合雜交三角形單元，并通過數(shù)值實驗驗證理論結(jié)果的正確性。結(jié)果表明：分片線性溫度梯度插值的散度(熱源)與非協(xié)調(diào)溫度插值是加權(quán)能量正交的；組合雜交三角形元剛度矩陣等同于協(xié)調(diào)的三角形線性元剛度矩陣，即非協(xié)調(diào)部分無溫度增強特性。

組合雜交元；三角形單元；熱傳導(dǎo)；能量正交；剛度矩陣；溫度增強

0 引 言

隨著航空航天、汽車、醫(yī)療等領(lǐng)域尖端技術(shù)的發(fā)展，僅研究材料的力學(xué)行為已不能滿足實際應(yīng)用的需求，需要關(guān)注材料的多物理場(例如熱、力、電、磁)耦合行為[1-3]，因此相關(guān)物理量(例如溫度、位移等)及其空間梯度、時間變化率等的計算精度都非常重要。

針對熱傳導(dǎo)-輻射問題，Z.Yang等[4-6]發(fā)展了周期及隨機復(fù)合材料的高階多尺度分析方法。針對彈性力學(xué)問題，T.Zhou[7-8]和聶玉峰等[9]建立了其組合雜交變分原理。作為穩(wěn)定化的變分原理，為應(yīng)力離散空間的優(yōu)化設(shè)計提供了很大便利，并成功建立了求解彈性力學(xué)問題的高性能四邊形單元、六面體單元以及板單元等[10-12]。在探索熱力耦合問題[13-15]的高性能組合雜交有限元求解算法之前，有必要先探索熱傳導(dǎo)問題的有效求解算法。

當(dāng)材料為各向同性，熱傳導(dǎo)問題的數(shù)學(xué)模型簡化為Poisson方程。不同于組合雜交矩形元[16]，本文建立Poisson方程的用非協(xié)調(diào)模式增強溫度插值函數(shù)的組合雜交三角形元，論證其分片線性溫度梯度插值的散度(熱源)與非協(xié)調(diào)溫度插值加權(quán)能量正交關(guān)系。

1 組合雜交變分原理

Poisson方程邊值問題為

(1)

(2)

式中：

a(σ,τ)=(σ,τ)Ωm

U=Uc⊕UI

UI|Ωm={span(Bubbles)}

α為組合參數(shù)，α∈(0,1)；Th為區(qū)域Ω的有限元剖分，Th={Ωm}。

對于網(wǎng)格Th，令Γh，Uh為相應(yīng)區(qū)域剖分的有限元離散空間，滿足Γh?Γ,Uh?U，則對上述問題有如下離散形式：

求(σh,uh)∈Γh×Uh，使得

αb2(σh,v)-b1(σh,vI)+(1-α)d(uh,v)=(f,v)

(?v∈Uh)

(3)

αa(σh,τ)-αb2(τ,uh)+b1(τ,uhI)=0

(?τ∈Γh)

(4)

2 場函數(shù)的插值近似

對單元Ωm，設(shè)pi(xi,yi)(i=1,2,3)為三角形單元按逆時針方向排列的三個頂點，(λ1,λ2,λ3)為三角形單元上任一點p(x,y)的面積坐標，Δ為三角形單元面積，直角坐標和面積坐標有如下關(guān)系：

λi=(ai+bix+ciy)/(2Δ)

(5)

式中：ai=xjyk-xkyj,bi=yj-yk,ci=xk-xj,i,j,k輪換,i=1,2,3。

溫度v=vc+vI的插值函數(shù)為

(6)

vI為非協(xié)調(diào)Bubble，將線性插值豐富為完全二次多項式以提高逼近精度。

溫度梯度τ的插值函數(shù)為分片線性多項式：

τ=(λ1I2,λ2I2,λ3I2)β

(7)式中：β=[τx(p1)τy(p1)τx(p2)τy(p2)τx(p3)τy(p3)]T，為結(jié)點溫度梯度參數(shù)；I2為2×2的單位矩陣。

為計算方便，定義矩陣Li=[bici]T,i=1,2。

3 剛度矩陣分析

式(4)中測試函數(shù)空間Γh是分片定義的，因此在每個單元Ωm上可以由式(4)解出溫度梯度σ，即用溫度uc和uI表示溫度梯度σ，再將σ的表達式帶入式(3)，進而可得到單元剛度矩陣。

在單元Ωm上，由式(6)～式(7)可得：

(8)

式中：

將式(8)帶入式(4)，由τ的任意性可知：

(9)

根據(jù)式(3)推導(dǎo)出單元剛度矩陣：

αb2m(σ,v)-b1m(σ,vI)

(10)

式中：

[0]3×3為零矩陣，即

(11)

式(11)說明分片線性溫度梯度插值的散度(熱源)與非協(xié)調(diào)溫度插值是能量正交的。

為了得到此單元和協(xié)調(diào)的三角形線性單元的等價性，還需計算矩陣E，根據(jù)式(3)和式(8)可得：

(12)

式中：

綜上可得：

αb2m(σ,v)-b1m(σ,vI)+(1-α)dm(u,v)

(13)

式中：

由于單元剛度矩陣的對稱性，可得單元剛度矩陣

(14)

式(14)表明三角形單元內(nèi)部自由度和頂點自由度無耦合，靜力凝聚內(nèi)部自由度后，D11保持不變，仍為協(xié)調(diào)的線性單元的剛度矩陣，由此證明溫度梯度插值為分片線性多項式，溫度插值函數(shù)為協(xié)調(diào)的線性部分和非協(xié)調(diào)部分的二次部分的三角形組合雜交元等價于基于最小勢能原理的協(xié)調(diào)線性三角形單元，非協(xié)調(diào)溫度插值部分無精度增強特性。

4 數(shù)值實驗

通過算例驗證上述分析結(jié)果，針對Poisson方程，即f=0，采用組合雜交元進行計算，并與線性元計算結(jié)果進行比較。外邊徑固定溫度為5K，內(nèi)邊徑給定熱流密度為10W/m2，熱傳導(dǎo)系數(shù)為20W/(m·K)，參考點為內(nèi)邊徑上任意一點，真解為6.648。計算區(qū)域為圓環(huán)，內(nèi)徑是3m，外徑是9m[18]。剖分如圖1所示，計算結(jié)果如表1所示。

表1 溫度絕對誤差

從表1可以看出：線性元與本文建立的組合雜交元計算溫度相同，數(shù)值結(jié)果與理論結(jié)果一致。

線性元在粗網(wǎng)格和細網(wǎng)格下的溫度分布如圖2所示。

從圖2可以看出：溫度沿半徑由內(nèi)向外降低，在內(nèi)邊徑達到最大值。

線性元在極細網(wǎng)格下沿x方向和y方向的熱流分布如圖3所示。

從圖3可以看出：線性元在極細網(wǎng)格下沿x方向和y方向的熱流分布均在內(nèi)邊徑處達到最大值。

5 結(jié) 論

對于Poisson方程，本文給出了其組合變分原理，并分析論證了組合雜交三角形元的加權(quán)能量正交關(guān)系。與四邊形單元所得結(jié)論不同，僅分片線性溫度梯度插值的散度(熱源)與非協(xié)調(diào)溫度插值是能量正交的。此時，組合雜交三角形元剛度矩陣等同于協(xié)調(diào)的三角形線性元剛度矩陣。構(gòu)造單元增強精度格式仍需著眼于突破恒等關(guān)系D12=[0]3×3。

[1]LiuW,QinY.Multi-physicscouplingmodelofcoalspontaneouscombustioninlongwallgobareabasedonmovingcoordinates[J].Fuel, 2017, 188: 553-566.

[2]SuH,RahmaniR,RahnejatH.Thermohydrodynamicsofbidirectionalgroovedrygassealswithslipflow[J].InternationalJournalofThermalSciences, 2016, 110: 270-284.

[3]BonitoA,DeVoreRA,NochettoRH.Adaptivefiniteelementmethodsforellipticproblemswithdiscontinuouscoefficients[J].SIAMJournalonNumericalAnalysis, 2013, 51(6): 3106-3134.

[4]YangZ,CuiJ,SunY.Transientheatconductionproblemwithradiationboundaryconditionofstatisticallyinhomogeneousmaterialsbysecond-ordertwo-scalemethod[J].InternationalJournalofHeatandMassTransfer, 2016, 100: 362-377.

[5]YangZ,CuiJ,SunY,etal.Multiscaleanalysismethodforthermo-mechanicalperformanceofperiodicporousmaterialswithinteriorsurfaceradiation[J].InternationalJournalforNumericalMethodsinEngineering, 2016, 105(5): 323-350.

[6]YangZ,CuiJ,ZhouS.Thermo-mechanicalanalysisofperiodicporousmaterialswithmicroscaleheattransferbymultiscaleasymptoticexpansionmethod[J].InternationalJournalofHeatandMassTransfer, 2016, 92: 904-919.

[7]ZhouT.Finiteelementmethodbasedoncombinationof“saddlepoint”variationalformulations[J].ScienceinChinaSeriesE:TechnologicalSciences, 1997, 40(3): 285-300.

[8]ZhouT.Stabilizedhybridfiniteelementmethodsbasedonthecombinationofsaddlepointprinciplesofelasticityproblems[J].MathematicsofComputation, 2003, 72(244): 1655-1673.

[9] 聶玉峰, 周天孝, 聶鐵軍. 三角形單元協(xié)調(diào)與非協(xié)調(diào)位移的能量正交關(guān)系[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 1999, 20(6): 619-624.NieYufeng,ZhouTianxiao,NieTiejun.Theenergyorthogonalrelationbetweenconformingandnon-conformingdisplacementsoftriangularelement[J].AppliedMathematicsandMechanics, 1999, 20(6): 619-624.(inChinese)

[10]ZhouTX,NieYF.Combinedhybridapproachtofiniteelementschemesofhighperformance[J].InternationalJournalforNumericalMethodsinEngineering, 2001, 51(2): 181-202.

[11] 聶玉峰, 周天孝. 高性能八節(jié)點六面體組合雜交元[J]. 數(shù)值計算與計算機應(yīng)用, 2003, 24(3): 231-240.NieYufeng,ZhouTianxiao. 8-nodehexahedroncombinedhybridelementwithhighperformance[J].JournalonNumericalMethodsandComputerApplications, 2003, 24(3): 231-240.(inChinese)

[12]ZhouTX,XieXP.Zeroenergy-errormechanismofthecombinedhybridmethodandimprovementofAllman’smembraneelementwithdrillingd.o.f.’s[J].InternationalJournalforNumericalMethodsinBiomedicalEngineering, 2004, 20(3): 241-250.

[13]RenB,QianJ,ZengX,etal.Recentdevelopmentsonthermo-mechanicalsimulationsofductilefailurebymeshfreemethod[J].ComputerModelinginEngineering&Sciences, 2011, 71(3): 253-278.

[14]YangZ,CuiJ.Thestatisticalsecond-ordertwo-scaleanalysisfordynamicthermo-mechanicalperformancesofthecompositestructurewithconsistentrandomdistributionofparticles[J].ComputationalMaterialsScience, 2013, 69: 359-373.

[15]GuanX,YuH,TianX.Astochasticsecond-orderandtwo-scalethermo-mechanicalmodelforstrengthpredictionofconcretematerials[J].InternationalJournalforNumericalMethodsinEngineering, 2016, 108(8): 885-901.

[16] 聶玉峰, 張玲, 王惠玲. 組合雜交Wilson矩形單元的加權(quán)能量正交關(guān)系[J]. 陜西師范大學(xué)學(xué)報: 自然科學(xué)版, 2014, 42(6): 26-30.

Nie Yufeng, Zhang Ling, Wang Huiling. The weighted energy orthogonal relation of combined hybrid Wilson rectangular element[J]. Journal of Shaanxi Normal University: Natural Science Edition, 2014, 42(6): 26-30.(in Chinese)

[17] Brezzi F, Fortin M. Mixed and hybrid finite element method[M]. Berlin Heidelberg: Springer-Verlag, 1991.

[18] Macneal R H, Harder R L. A proposed standard set of problems to test finite element accuracy[J]. Finite Elements in Analysis and Design, 1985, 1(1): 3-20.

(編輯：趙毓梅)

The Weighted Energy Orthogonal Relation of Combined Hybrid Triangular Element

Zhang Ling, Zhang Weiwei

(School of Natural and Applied Sciences, Northwestern Polytechnical University, Xi’an 710129, China)

Variation principles can enhance the stability of numerical solution with the combined hybrid finite element method. Combined hybrid finite element method of heat transfer equation is built on the basis of the domain decomposition technique. The combined hybrid triangular element, in which the temperature gradient is interpolated by linear polynomials on each element, but the temperature is interpolated by the sum of the linear polynomials and the non-conforming quadratic polynomials, is given. The numerical experiments are carried out to verify the accuracy of the theoretical results. The results indicate that the divergence of piecewise linear temperature gradient interpolation and the non-conforming temperature interpolation are of the weighted energy orthogonal relation. The stiffness matrix of this element is equivalent to the conforming triangular linear element, and the non-conforming parts have no contribution to temperature evaluation.

combined hybrid element; triangular element; heat transfer; energy orthogonal; stiffness matrix; enhanced temperature

2016-12-06；

2016-12-20

國家自然科學(xué)基金(11471262，11501450) 西北工業(yè)大學(xué)博士論文創(chuàng)新基金(CX201524)

張玲,shanshi211@163.com

1674-8190(2017)01-073-05

[O242.21]

A

10.16615/j.cnki.1674-8190.2017.01.011

張 玲(1987-)，女，博士研究生。主要研究方向：熱力學(xué)問題的高效有限元方法。

張偉偉(1986-)，女，博士，講師。主要研究方向：自適應(yīng)并行有限元方法。