二維不可壓縮層流模型數值解與穩定性

肖潤梅,明亞東,李秀蘭

(山西大同大學數學與計算機科學學院,中國 大同 037009)

層流流體穩定性問題的研究長期以來一直倍受數學家的關注．當雷諾數在某個特定范圍內逐步增大時，流體將經歷一個由層流向紊流轉換的特殊過程[1]，即由一個穩定的、亞臨界的層流流體轉化成為一個超臨界的紊流流體．為此，通過對流體雷諾數和波長的計算來預測層流的穩定性問題變得尤為重要．目前，人們通過對層流流體加入微小外部擾動來預測層流流體穩定性的方法已被廣泛地接受[2]．如果在此過程中，外部擾動能夠逐漸減弱直至消失，那么流體將繼續保持穩定．反之，如果這些外部的擾動隨著時間的推移逐漸增大，那么流體將經歷一個由層流向紊流流動轉化的過渡過程．

層流流體的穩定性問題可歸結為Orr-Sommerfeld方程求解問題，該方程的求解在數學界引起了廣泛的興趣．Davey[3]通過有限元的方法來求解Orr-Sommerfeld方程，取得了良好的結果．然而由于有限元方法本身的局限性，運算精度會受到很大的限制．如在計算誤差(Δx)p漸進減小的過程中，運用有限元方法只能得到一些有限的p值，同時計算時間過長，不適合廣泛運用于求解此類問題．切比雪夫多項式在逼近理論中有重要的應用．第一類切比雪夫多項式的根(被稱為切比雪夫節點)可以用于多項式插值．相應的插值多項式能最大限度地降低龍格現象，并且提供多項式在連續函數的最佳一致逼近．Orszag[4]運用切比雪夫多項式以及QR矩陣特征值算法所計算的結果有著較高的計算精度以及較少的計算時間．

當前，譜方法已成為計算流體力學領域的三大主流數值方法之一，并受到廣大學者的關注．它是上世紀70年代發展起來的一種數值求解偏微分方程的方法，該方法在理論上具有“無窮階”收斂性，可采用快速算法，現已被廣泛用于氣象、物理、力學等諸多領域，成為繼差分法和有限元法之后又一種重要的數值方法[5-7]．由于譜方法的“無窮階”收斂性且具備高的計算精度，相對于有限元方法和切比雪夫多項式的方法，它能夠對誤差形成指數型收斂，提高了計算精度并縮短了計算時間，為求解Orr-Sommerfeld方程提供了新的思路.

本文運用譜方法的數值分析，以Poiseuille流動為例對Orr-Sommerfeld方程進行了展開與數值計算，并利用計算結果對層流流體的穩定性進行了分析討論，并驗證了譜方法的高精度性和快速收斂性．

1 Orr-Sommerfeld方程數值計算

關于影響層流流動穩定性的因素，主要認為是當流體流經渦輪壓氣機和管道時，由于入口處粗燥的壁面所造成的邊界層流的影響以及外部諸因素對流體的干擾．因此，研究層流流動例如Poiseuille流動的穩定性問題時，可歸結為求解Orr-Sommerfeld方程的數值解問題．忽略加速流體過渡過程的外部干擾以及由于壓力差造成的過渡過程，并假設其為零壓力梯度場的穩定流體，即限制研究對象為：具有一定雷諾數的，零壓力場的，二維不可壓縮穩定流體．采用線化小擾動理論，可得到一個四階線性的Orr-Sommerfeld常微分方程

(1)

φ(1)=0；φ′(1)=0，

φ(-1)=0；φ′(-1)=0，

(2)

無量綱Poiseuille流動的速度分布為U(y)=1-y2.

對于上述方程，我們并不能求得它的解析解，只能通過數值計算的方法來求出它的數值解．譜方法的主旨思想是將方程在固定的物理空間進行了離散化．通過這樣的方式，就可以更加輕松地去處理方程中的一些非線性的變量．首先，通過使用多項式Chebyshev-Gauss-Lobatto nodes[8]對方程的一階導數進行插值可以得到

(3)

由于在區間x∈[-1,1]上φ(x)為連續函數，所以通過構建N階多項式Lk(xj)對φ(x)進行插值Lk(xj)=δjk.

(4)

結合邊界條件，可得

(5)

最后得到方程

(6)

在此，需要對切比雪夫微分矩陣D引入三角恒等式，可得

(7)

(8)

(9)

(10)

把以上等式代入Orr-Sommerfeld 方程可以得到2個不同的矩陣

A=(D4-2D2+I)/Re-2iaI-idiag(1-x2)(aD2-a3I),

B=-iaD2+a3I.

(11)

D4、D2分別為矩陣D的四階和二階導數矩陣；I為單位矩陣．那么從方程(11)得到特征值c滿足

Aφk=cBφk.

(12)

2 數值計算結果分析

圖1 Poiseuille流動穩定性計算圖

求解方程(12)的特征值c，需要給定α和Re值．如果ci<0，那么層流流動將保持穩定，相反地，如果ci>0，那么層流流體就會向紊流進行轉化，ci=0表示臨界穩定狀態．圖1表示了ci隨Re和α值變化的曲線(具體計算程序見附錄1)．

圖中分別顯示了ci=0.0，ci=0.04，ci=0.06的等值線．從圖中可以看出由等值線ci=0.0所包圍的中間區域為不穩定區, 也就是說當Re和α的值在這個范圍內時, 層流流體是不穩定的，極易向紊流進行轉化．等值線ci=0.0外圍被視為穩定區域．同時可以得到臨界的雷諾數值為Recrit=5 773，也就是說如果當流體的雷諾數小于5 773時，無論α取值多少，流體始終保持層流流動．同時在臨界雷諾數處，首先產生不穩定流體的αc=1.02．對于本文得到的臨界數值可以與其他文獻所得到的值進行比對，Thomas[9]運用有限元的方法得到Recrit≈5 780，αc=1.026．Orszag[4]通過切比雪夫多項式的方法得到Recrit≈5 772，αc=1.021．由此驗證了通過譜方法求解Orr-Sommerfeld equation有著很高的計算精度以及較少的計算時間，為解決此類方程提供了新的思路．

附錄1

TheComputerprogramming

function [x, DL] = cheb2nddif(N, Dem)

I = eye(N);

Lk = logical(I);

N1 = floor(N/2);

N2 = ceil(N/2);

k = [0:N-1]′;

theta = k*pi/(N-1);

x = sin(pi*[N-1:-2:1-N]′/(2*(N-1)));

t = repmat(theta/2,1,N);

Dx = 2*sin(t′+t).*sin(t′-t);

Dx = [Dx(1:N1,:); -flipud(fliplr(Dx(1:N2,:)))];

Dx(Lk) = ones(N,1);

c = toeplitz((-1).^k);

c(1,:) = c(1,:)*2; c(N,:) = c(N,:)*2;

c(:,1) = c(:,1)/2; c(:,N) = c(:,N)/2;

delta = 1./Dx;

delta(Lk) = zeros(N,1);

D = eye(N);

for i = 1:Dem

D = i*delta.*(c.*repmat(diag(D),1,N) - D);

D(Lk) = -sum(D′);

DL(:,:,i) = D;

End

function [x, D4] = cheb4thdif(N)

I = eye(N-2);

Lk = logical(I);

N1 = floor(N/2-1);

N2 = ceil(N/2-1);

k = [1:N-2]′;

theta = k*pi/(N-1);

x = sin(pi*[N-3:-2:3-N]′/(2*(N-1)));

r = [sin(theta(1:N1)); flipud(sin(theta(1:N2)))];

gamma = r.^4;

lamda1 = -4*r.^2.*x./gamma;

lamda2 = 4*(3*x.^2-1)./gamma;

lamda3 = 24*x./gamma;

lamda4 = 24./gamma;

B = [lamda1′; lamda2′; lamda3′; lamda4′];

t = repmat(theta/2,1,N-2);

Dx = 2*sin(t′+t).*sin(t′-t);

Dx = [Dx(1:N1,:); -flipud(fliplr(Dx(1:N2,:)))];

Dx(Lk) = ones(N-2,1);

rr = r.^2.*(-1).^k;

R = rr(:,ones(1,N-2));

c = R./R′;

delta = 1./Dx;

delta(Lk) = zeros(size(x));

Z = delta′;

Z(Lk) = [];

Z = reshape(Z,N-3,N-2);

W = ones(N-3,N-2);

D = eye(N-2);

for i = 1:4

W = cumsum([B(i,:); i*W(1:N-3,:).*Z]);

D = i*delta.*(c.*repmat(diag(D),1,N-2)-D);

D(Lk) = W(N-2,:);

DL(:,:,i) = D;

end

D4 = DL(:,:,4);

N=32;

i = sqrt(-1);

a=0.1:0.01:1.2;

R=1:500:120000;

for k=1:length(a)

for j=1:length(R)

[x,DM] = cheb2nddif(N+2,2);

D2 = DM(2:N+1,2:N+1,2);

[x,D4] = cheb4thdif(N+2);

I = eye(size(D4));

A = (D4-2*D2+I)/R(j)-2*a(k)*i*I-i*diag(1-x.^2)*(a(k)*D2-a(k)^3*I);

B = -i*a(k)*D2+i*a(k)^3*I;

e = eig(A,B);

[m,l] = min(abs(imag(e)));

res(k,j)=imag(e(l));

end

參考文獻：

[1] SCHOBEIRI M T. Fluid mechanics for engineers[M]. Berlin:Springer-Verlag, 2010.

[2] DAVEY A, DRAZIN P G. The stability of poiseuille flow in the pipe[J]. J Fluid Mech, 1969,36(02):209-218.

[3] DAVEY A, NGUYEN H P F. Finite-amplitude stability of pipe flow[J]. J Fluid Mech, 1971,45(04):701-720.

[4] ORSZAG S A. Accurate solution of the orr-sommerfeld stability equation[J]. J Fluid Mech, 1971,50(04):689-703.

[5] 尚新民. 譜方法的數值分析[M]. 北京：科學出版社, 2000.

[6] CANUTO C, HUSSAINI M Y, QUARTERONI A,etal. Spectral methods in fluid dynamics[M]. Berlin:Springer-Verlag, 1993.

[7] BALTENSPERGER R, TRUMMER M R. Spectral differencing with a twist[J]. Siam J Sci Comput, 2002,24(5):1465-1487.

[8] DON W, SOLOMONOFF S. Accuracy and speed in computing the chebyshev collocation derivative[J]. Siam J Sci Comput, 1995,16(6):1253-1268.

[9] THOMAS L H. The stability of plane poiseuille flow[J]. Phys Rev, 1953,91(4):780-783.