阿妮柯孜·奧斯曼，馮新龍，劉德民

(新疆大學數(shù)學與系統(tǒng)科學學院，新疆烏魯木齊 830017)

0 引言

本文考慮了求解非定常微極流體方程的速度校正投影方法, 設區(qū)域Ω ?Rd(d=2 或3) 是具有光滑邊界?Ω的凸連通區(qū)域, 非定常微極流體方程為:

方程中未知函數(shù)是速度u, 角速度w, 壓力p.ν ＞0 是通常的黏度系數(shù), νr, c0, ca, cd是與非對稱應力張量有關的系數(shù), ca, cd和c0滿足不等式c0+cd-ca＞0.方程中ν0=ν+νr, c1=ca+cd, c2=c0+cd-ca.當空間維數(shù)d=2時, 速度u=(u1,u2,0), 角速度w=(0,0,w3).當d=3 時, u=(u1,u2,u3), w=(w1,w2,w3).

微極流體在工業(yè)和工程中有重要的應用價值, 許多學者對其進行了研究.Ortega-Torres 等用全離散加罰有限元法證明了速度、壓力和角速度的最優(yōu)誤差估計[1].Nochetto 等提出了微極流體方程的一階半隱式全離散有限元方法, 使速度和角速度解耦求解, 并給出了算法的無條件穩(wěn)定性和最優(yōu)收斂階[2].Salgado 研究了微極流體方程的全離散分步投影有限元方法, 在連續(xù)性方程中加入grad-div 穩(wěn)定項可以改善質量守恒[3].Galdi 等研究了微極流體方程弱解的存在唯一性[4].非定常微極流體方程數(shù)值方法研究的主要困難包括強非線性、無散度限制及多場耦合性等, 使得變分問題成為典型的鞍點問題, 理論分析和數(shù)值離散時速度變量和壓力變量需滿足inf-sup 條件.投影方法在類似于文獻[5-6]所提的偏微分方程數(shù)值解中有著廣泛的應用,就是把原來的鞍點問題轉化成橢圓問題進行求解, 從而降低了求解規(guī)模并改善了算法穩(wěn)定性.Jiang 等提出了求解非定常微極流體方程的基于壓力校正的投影方法并給出了算法的穩(wěn)定性和誤差估計[7], 第一步先求速度和角速度, 第二步是把速度投影到無散度空間.本文在文獻[8] 的基礎上給出了微極流體方程的速度校正方法, 算法的第一步是投影, 第二步是校正, 第三步求角速度.與文獻[7] 的結果相比, 速度校正方法是一種全解耦的方法, 它不需要壓力初值,也不需要速度和壓力滿足inf-sup 條件[9－10], 并推廣到了非齊次Dirichlet 邊界條件.

1 預備知識

這里引入一般的標量Sobolev 空間Hm(Ω) (m=0,1,2) 和向量Sobolev 空間Hm(Ω)d, 在不引起混淆的情形下記它們的范數(shù)為‖·‖m.特別的, 當m=0 時, 令L2(Ω)=H0(Ω) 和L2(Ω)d=H0(Ω)d, 對應的內積和范數(shù)為(·,·), ‖·‖, 記:

2 時間半離散速度校正方法

2.1 一階時間半離散速度校正方法

下面估計式(10) 右邊的項:

將上述估計代入式(10), 并把式(9) 和式(10) 相加, 則有:

當k=0,1,···,m-1 時, 對上面的不等求和并用離散形式的Gronwall 引理, 則結論成立.

為了進行收斂性分析, 引入以下符號:

定理2 令(uk+1,pk+1;～uk+1,wk+1) 是一階時間半離散速度校正法的解, 0 ≤m ≤N, 則成立:

證明將t=tk+1代入式(1), 則:

將式(13) 的第一個方程減去式(3), 式(13) 的第二個方程減去式(6) 得到:

用u(tk+1)-ν0ΔtΔu(tk+1) 與式(5) 作差, 得到：

將誤差方程式(14) 改寫為:

2.2 二階時間半離散速度校正方法

步驟1: 求uk+1, pk+1滿足:

步驟2: 求～uk+1滿足:

式(25) 也可以改寫為:

證明與定理1 類似.

3 全離散速度校正方法

設τh={K} 是Ω 的一致正則三角網格, 且網格大小h=maxK∈τh{diam(K)}, 定義以下離散子空間:

3.1 一階全離散速度校正方法

證明與定理1 類似.

3.2 二階全離散速度校正方法

4 數(shù)值算例

本節(jié)對二維微極流體方程結構化和非結構化網格上協(xié)調有限元逼近問題的收斂性進行試驗,并對三維微極流體方程進行了結構化網格上的收斂性試驗.分別給出了能量穩(wěn)定性和微極流體在軸承潤滑中的應用.

4.1 收斂性試驗

設Ω=[0,1]d(d=2, 3), 終止時刻T=0.1 并且ν=νr=ca=cd=1.一階和二階速度校正投影方法分別選取Δt=h2, Δt=h.

4.1.1 二維微極流體方程在結構化網格上的收斂性試驗

考慮具有光滑解析解的流動問題, 方程(1) 的右端由精確解決定.

表1 和表2 分別給出了二維微極流體方程的一階和二階速度校正法的速度、角速度和壓力變量在結構化網格上的數(shù)值結果.由表1 和表2 可知,一階和二階格式得到的速度、角速度、壓力都能達到最優(yōu)收斂階.由此可知,速度校正方法是有效和可靠的.

表1 一階速度校正法在Δt=h2 時的誤差和收斂階

表2 二階速度校正法在Δt=h 時的誤差和收斂階

4.1.2 二維微極流體方程在非結構化網格上的收斂性試驗

考慮非結構化網格上一階和二階速度校正方法的誤差和收斂階.其中, 網格特征尺寸分別為h=(h1,h2,h3,h4,h5)=(0.270 282, 0.146 33, 0.073 319, 0.035 511, 0.018 824 5).選取精確解為:

表3 和表4 分別給出了微極流體方程一階和二階速度校正法的速度、角速度和壓力變量在非結構化網格上的數(shù)值結果.由表3 和表4 可知, 一階和二階格式得到的速度、角速度、壓力都能達到最優(yōu)收斂階.由此可知, 速度校正方法是有效和可靠的.

表3 一階速度校正方法在Δt=0.000 1 時的誤差和收斂階

表4 二階速度校正方法在Δt=0.000 1 時的誤差和收斂階

4.1.3 三維微極流體方程在結構化網格上的收斂性試驗

為了更好地說明速度校正投影法的有效性和可靠性, 給出了三維微極流體方程對應的收斂性檢驗.精確解如下:

如表5 和表6 所示, 給出了三維微極流體方程一階和二階速度校正法在非結構化網格上的數(shù)值結果.由表5 和表6 可知, 一階和二階格式得到的速度、角速度、壓力都能達到最優(yōu)收斂階.

表6 二階速度校正法在Δt=h 時的誤差和收斂階

4.2 能量穩(wěn)定性

考慮含有周期邊界條件的求解區(qū)域Ω=[0,1]2, 分析在給定初值條件下, 外力為f =g=0 時系統(tǒng)的能量耗散過程.選擇初始條件

方程(1) 中, 當ν=νr=ca=cd=c0=1, T=2, h=1/60 時.給出了一階和二階速度校正投影方法的能量耗散圖.

考慮系統(tǒng)能量En:=‖un‖2+‖vn‖2+‖wn‖2, 圖1 給出了系統(tǒng)在不同時間步長下的能量耗散過程.四個能量曲線在所有時間步長上都呈現(xiàn)單調衰減.其中Δt=0.1,0.05,0.025,0.012 5.

圖1 在不同時間步長下一階(左)和二階(右)格式的系統(tǒng)能量

4.3 微極性流體在軸承潤滑中的應用

微極流體方程在軸承潤滑問題中有著極其重要的應用.在本節(jié)中, 將具有非線性項的微極流體方程(u·?)u和j(u·?)w 加到方程(1) 的動量方程和角動量方程中, 并將對應的離散形式(uk·?)uk+1和j(uk·?)wk+1加到所有算法的第一步和第三步.顯然, 非線性項的引入不會給數(shù)值分析帶來額外的困難.其中d=2, ν=νr=ca=cd=c0=1.選取P1b-P1-P1b 有限元對.給出了不同轉速Ωr下水平速度u1,垂直速度u2,壓力和角速度用一階速度校正法求解微極流體方程, 并給出微極流體方程在軸承潤滑中應用的演化圖.設外圓半徑r1=1, 內圓半徑r2=0.3, Γ1為外邊界, Γ2為內邊界.Γ1滿足齊次邊界條件u=0, Γ2滿足u=(-Ωrr2sinθ,-Ωrr2cosθ), w=0,這里θ=arctan(y/x).分別取Ωr=100,500,1 000, Δt=0.01, T =1.

如圖2 所示,隨著轉速Ωr的增大,流體的水平速度u1、垂直速度u2、角速度、壓力也增大,同時流體的非對稱性也逐漸明顯.將線速度分量的大小與角速度的大小進行比較, 可知流體的極性效應仍然很小.當Ωr=1 000時我們的算法依然有效.

圖2 Ωr=100, 500, 1 000時的水平速度、垂直速度、角速度和壓力的演化

5 結論

本文提出了微極流體方程的速度校正投影法,針對線性的微極流體方程給出了穩(wěn)定性分析和有限元逼近的誤差估計結果,數(shù)值算例中, 將格式推廣到非線性微極流體方程.數(shù)值算例中分別給出了解析解問題、周期邊界的初值問題以及非齊次邊界的同心軸承潤滑問題.速度校正法可以有效降低求解規(guī)模, 數(shù)值算例驗證了相應的結論.在實際問題的模擬中還需要根據(jù)具體的物理工況進行參數(shù)設置和格式調整, 這將在后期進一步研究.