摘? 要：斯托克斯方程是流體力學(xué)中描述粘性牛頓流體運(yùn)動(dòng)的方程組，在航空航天、流體力學(xué)、化工工程、生物工程等與流體有關(guān)的所有領(lǐng)域都有著非常重要的作用和廣泛的應(yīng)用。有限差分方法通過離散物理問題，從而以代數(shù)方程式的形式來求解出誤差較小的近似數(shù)值解。文章給出斯托克斯第一問題使用有限差分法和MATLAB矩陣運(yùn)算的計(jì)算方法，通過對網(wǎng)格進(jìn)行二維剖分細(xì)化，借用計(jì)算機(jī)的強(qiáng)大算力得出穩(wěn)定性條件下的速度曲線，并使用相關(guān)算例進(jìn)一步確定擴(kuò)散數(shù)的閾值。

關(guān)鍵詞：有限差分法；斯托克斯問題；MATLAB

中圖分類號：TP391? 文獻(xiàn)標(biāo)識(shí)碼：A? 文章編號：2096-4706（2023）20-0058-04

MATLAB Implementation of the Finite Difference Method for Stokes First Problem

CHEN Geng

（North Minzu University， Yinchuan? 750030， China）

Abstract： The Stokes equations are a set of equations describing the motion of viscous Newtonian fluids in fluid mechanics and have a very important role and wide application in all fields related to fluids such as aerospace， fluid mechanics， chemical engineering and bioengineering. The finite difference method solves approximate numerical solutions with small errors in algebraic equation form by discretizing the physical problem. This paper gives a computational approach to the Stokes first problem using the finite difference method and MATLAB matrix operations to derive the velocity profile under stability conditions by refining the grid with a two-dimensional profile， borrowing the power of computer arithmetic， and using relevant arithmetic examples to further determine the threshold value of the diffusion number.

Keywords： finite difference method; Stokes problem; MATLAB

0? 引? 言

隨著科技的進(jìn)步和計(jì)算機(jī)技術(shù)的快速發(fā)展，在自然科學(xué)和工程計(jì)算領(lǐng)域誕生了許多實(shí)際問題，人們發(fā)現(xiàn)這些問題可以通過微分方程恰當(dāng)?shù)乇硎境鰜韀1]。Stokes方程是流體力學(xué)中描述粘性牛頓流體運(yùn)動(dòng)的方程組，在航空航天，流體力學(xué)，化工工程，生物工程等與流體有關(guān)的所有領(lǐng)域都有著非常重要的作用和廣泛的應(yīng)用。同時(shí)Stokes問題在不可壓縮流體問題的穩(wěn)定性的分析上起著至關(guān)重要的作用?，F(xiàn)如今，對于微分方程的解，可以歸結(jié)為兩大類，一類是精確解，另一類是數(shù)值解。由于精確解只有在一些特殊的情況下才可以求出，因此研究微分方程的數(shù)值解就變得非常重要。對于偏微分方程數(shù)值解的研究已經(jīng)有很長的歷史，目前應(yīng)用比較廣泛的方法主要有有限差分法、有限元法、譜方法等[2，3]。

近年來，Stokes問題的數(shù)值求解方法一直是國內(nèi)外學(xué)者研究的熱點(diǎn)，有限差分法是求解微分方程的一種成熟、高效的方法，其主要思想為以差商近似導(dǎo)數(shù)，具體步驟可分為：對方程的求解區(qū)域進(jìn)行網(wǎng)格剖分；對微分算子進(jìn)行離散；建立以網(wǎng)格節(jié)點(diǎn)上的值為未知數(shù)的代數(shù)方程組等[2]。這種方法在規(guī)范幾何形狀區(qū)域內(nèi)誤差在合理范圍內(nèi)且運(yùn)算速度快，更具優(yōu)勢。MATLAB基于矩陣運(yùn)算，具有強(qiáng)大的數(shù)值運(yùn)算能力。使用MATLAB求解微分方程，近些年有大量的研究[4-7]，王憶鋒等人運(yùn)用有限差分法直接求解二維下的泊松方程[8]，廖臣等人則在有限差分法對偏微分方程求解下探究五點(diǎn)差分格式的算法[9]，巫衡竹等人使用有限差分法來求解Stokes問題的理論證明。因此探索Stokes問題新的數(shù)值計(jì)算方法及應(yīng)用具有重要的現(xiàn)實(shí)意義。

本文探討通過MATLAB使用有限差分法構(gòu)造差分格式，從而在剖分細(xì)化的網(wǎng)格上求解Stokes第一問題的誤差在合理范圍內(nèi)的數(shù)值解，并使用相關(guān)算例進(jìn)行結(jié)果檢驗(yàn)和穩(wěn)定性驗(yàn)證。

1? 斯托克斯第一問題在網(wǎng)格節(jié)點(diǎn)上的有限差分

假設(shè)有一塊無限大的平板浸沒在無界的靜止流體中，突然以速度u0沿其自身所在的平面運(yùn)動(dòng)起來，并且一直保持速度的大小和方向不變。請求解平板起動(dòng)后流體運(yùn)動(dòng)隨時(shí)間的變化過程，簡稱為平面壁突然加速問題。

考慮Navier-Stokes方程：

選擇y軸與平面壁重合并指向u0方向，由于平板是無限大的，我們假定在相同的xy平面內(nèi)速度是恒定的，故?u / ?x = ?u / ?y = 0，且無外力項(xiàng)。由文獻(xiàn)[10]知流場壓力為常數(shù)且處于同一水平高度，故p = 0。綜上簡化后的Navier-Stokes方程為：

有限差分方法通過離散物理問題，對網(wǎng)格進(jìn)行二維剖分細(xì)化，從而以代數(shù)方程式形式求解出近似數(shù)值解。接下來進(jìn)行網(wǎng)格剖分，h表示水平步長，Δt表示垂直步長，等間距劃分N×N個(gè)網(wǎng)格，網(wǎng)格節(jié)點(diǎn)標(biāo)記如圖1所示。

對函數(shù)u （z， t）做二維泰勒級數(shù)展開，寫作：

將上式中的前兩式相加，可得同時(shí)間t下的u對t的二階導(dǎo)數(shù)，同時(shí)選取向前差分格式寫出節(jié)點(diǎn)i的u對t的一階導(dǎo)數(shù)：

代入簡化后Navier-Stokes方程可得：

舍去上式中的o（Δt， h2）后整理可得：

其中j和j + 1表示時(shí)間t時(shí)刻的速度和t + Δt時(shí)刻的速度；i、i -1和i + 1表示網(wǎng)格細(xì)分后的節(jié)點(diǎn)速

度。由穩(wěn)定性準(zhǔn)則可定義擴(kuò)散數(shù)條件應(yīng)滿足：

如果擴(kuò)散數(shù)超過穩(wěn)定性標(biāo)準(zhǔn)，則解的震蕩誤差會(huì)增大。

2? 有限差分格式的求解

斯托克斯方程的邊值問題的有限差分法是利用網(wǎng)格狀模型上離散節(jié)點(diǎn)的數(shù)值解來逼近在時(shí)間序列下的真實(shí)解。在計(jì)算機(jī)環(huán)境允許的情況下，網(wǎng)格劃分越細(xì)密，則節(jié)點(diǎn)越多，離散化模型越精確。運(yùn)用MATLAB求解問題，則需將差分代數(shù)方程轉(zhuǎn)化為矩陣格式：KU = B，故需將上式整理為：3? 數(shù)值實(shí)驗(yàn)

假設(shè)在t = 0時(shí)刻之前的平板和流體均處于靜止?fàn)顟B(tài)，隨后平板開始沿著自身的平面以u0 = 15 m/s的速度運(yùn)動(dòng)，流體的運(yùn)動(dòng)粘度為α = 0.000 217 m2/s，設(shè)置步長分別為h = 0.001 m，Δt = 0.02 s或Δt = 0.002 s，由此計(jì)算斯托克斯第一問題的速度。

由已知數(shù)據(jù)可計(jì)算出擴(kuò)散數(shù)D在兩組不同的Δt步長下分別為：D1 = 4.34，D2 = 0.434，由文獻(xiàn)[10]知D1超過了穩(wěn)定性標(biāo)準(zhǔn)。本程序?yàn)楦阌谥庇^表現(xiàn)擴(kuò)散數(shù)是否震蕩，故編寫圖像展示如圖2和圖3所示。

計(jì)算結(jié)果如表2所示。

分析表3中的數(shù)值可以再次驗(yàn)證D1 = 4.34時(shí)的數(shù)值震蕩誤差較大，迭代后速度之間跳躍過大；D2 = 0.434時(shí)的數(shù)值結(jié)果符合穩(wěn)定性準(zhǔn)則。如果要減少D1的震蕩誤差，可適當(dāng)增大網(wǎng)格步長h或是進(jìn)一步細(xì)分時(shí)間步長Δt。

4? 結(jié)? 論

本文給出了斯托克斯第一問題方程的有限差分格式，利用MATLAB編制求解算法相關(guān)程序，并繪制出直觀圖像觀察數(shù)值解是否存在震蕩誤差。對MATLAB程序進(jìn)行數(shù)值實(shí)驗(yàn)，驗(yàn)證該方法的有效性并運(yùn)用計(jì)算機(jī)的巨大算力實(shí)現(xiàn)有限差分法計(jì)算流體力學(xué)，大大提高計(jì)算效率。

參考文獻(xiàn)：

[1] 吳澤康，韓文靜，耿勇，等.二維二階雙曲型方程隱式差分格式的穩(wěn)定性分析 [J].齊魯工業(yè)大學(xué)學(xué)報(bào)，2023，37（1）：75-80.

[2] 陳艷萍，魯祖亮，劉利斌.信息與計(jì)算科學(xué)叢書67：偏微分方程數(shù)值解法 [M].北京：科學(xué)出版社，2015.

[3] 戴嘉尊，邱建賢.微分方程數(shù)值解法：第2版 [M].南京：東南大學(xué)出版社，2012.

[4] 唐洪浪，桂現(xiàn)才.用MATLAB符號工具箱編程求常微分方程的通解 [J].洛陽師范學(xué)院學(xué)報(bào)，2005（2）：81-84.

[5] 趙德奎，劉勇.MATLAB在有限差分法數(shù)值計(jì)算中的應(yīng)用 [J].四川理工學(xué)院學(xué)報(bào)：自然科學(xué)版，2005（4）：61-64.

[6] 馮立偉.熱傳導(dǎo)方程幾種差分格式的MATLAB數(shù)值解法比較 [J].沈陽化工大學(xué)學(xué)報(bào)，2011，25（2）：179-182+191.

[7] 馮立偉，王選鶴，馬瑩.一種非定常不可壓N-S方程的不等階插值FDSD解法 [J].沈陽化工學(xué)院學(xué)報(bào)，2009，23（1）：80-84+96.

[8] 王憶鋒，唐利斌.利用有限差分和MATLAB矩陣運(yùn)算直接求解二維泊松方程 [J].紅外技術(shù)，2010，32（4）：213-216+230.

[9] 廖臣，祝大軍，劉盛綱.五點(diǎn)差分格式求解泊松方程并行算法的研究 [J].電子科技大學(xué)學(xué)報(bào)，2008（1）：81-83+127.

[10] MAGRAB E B，AZARM S，BALACHANDRAN B.MATLAB原理與工程應(yīng)用 [M].高會(huì)生，李新葉，胡智奇，等譯.北京：電子工業(yè)出版社，2002.

作者簡介：陳庚（2000—），男，滿族，河北承德人，本科在讀，研究方向：計(jì)算流體力學(xué)。

收稿日期：2023-03-20

基金項(xiàng)目：寧夏回族自治區(qū)大學(xué)生創(chuàng)新訓(xùn)練項(xiàng)目（S202211407040）