張玲，張偉偉

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

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

張玲，張偉偉

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

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

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

0 引 言

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

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

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

1 組合雜交變分原理

Poisson方程邊值問(wèn)題為

(1)

(2)

式中：

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

U=Uc⊕UI

UI|Ωm={span(Bubbles)}

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

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

求(σ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 場(chǎng)函數(shù)的插值近似

對(duì)單元Ωm，設(shè)pi(xi,yi)(i=1,2,3)為三角形單元按逆時(shí)針?lè)较蚺帕械娜齻€(gè)頂點(diǎn)，(λ1,λ2,λ3)為三角形單元上任一點(diǎn)p(x,y)的面積坐標(biāo)，Δ為三角形單元面積，直角坐標(biāo)和面積坐標(biāo)有如下關(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，將線性插值豐富為完全二次多項(xiàng)式以提高逼近精度。

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

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

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

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

3 剛度矩陣分析

式(4)中測(cè)試函數(shù)空間Γh是分片定義的，因此在每個(gè)單元Ωm上可以由式(4)解出溫度梯度σ，即用溫度uc和uI表示溫度梯度σ，再將σ的表達(dá)式帶入式(3)，進(jìn)而可得到單元?jiǎng)偠染仃嚒?/p>

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

(8)

式中：

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

(9)

根據(jù)式(3)推導(dǎo)出單元?jiǎng)偠染仃嚕?/p>

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

(10)

式中：

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

(11)

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

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

(12)

式中：

綜上可得：

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

(13)

式中：

由于單元?jiǎng)偠染仃嚨膶?duì)稱性，可得單元?jiǎng)偠染仃?/p>

(14)

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

4 數(shù)值實(shí)驗(yàn)

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

表1 溫度絕對(duì)誤差

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

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

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

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

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

5 結(jié) 論

對(duì)于Poisson方程，本文給出了其組合變分原理，并分析論證了組合雜交三角形元的加權(quán)能量正交關(guān)系。與四邊形單元所得結(jié)論不同，僅分片線性溫度梯度插值的散度(熱源)與非協(xié)調(diào)溫度插值是能量正交的。此時(shí)，組合雜交三角形元?jiǎng)偠染仃嚨韧趨f(xié)調(diào)的三角形線性元?jiǎng)偠染仃?。?gòu)造單元增強(qiáng)精度格式仍需著眼于突破恒等關(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é)點(diǎn)六面體組合雜交元[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(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é)報(bào): 自然科學(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

國(guó)家自然科學(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é)問(wèn)題的高效有限元方法。

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