擴(kuò)散方程的數(shù)值差分解法

劉浩庭

摘要：擴(kuò)散現(xiàn)象是其初始密度不均勻分布引起的，它會(huì)對物質(zhì)粒子的分布狀態(tài)產(chǎn)生影響，最終達(dá)到物質(zhì)在空間均勻分布狀態(tài)。本研究通過分離變量法對給定條件的粒子濃度在一維空間分布下的擴(kuò)散現(xiàn)象的解析解進(jìn)行計(jì)算，同時(shí)結(jié)合使用歐拉法利用計(jì)算物理的方法對擴(kuò)散過程進(jìn)行了數(shù)值模擬。通過對比理論數(shù)據(jù)與模擬實(shí)驗(yàn)數(shù)據(jù)對等離子體一維擴(kuò)散現(xiàn)象進(jìn)行闡釋與討論。

Abstract： The diffusion phenomenon is caused by the difference of the initial density. It affects the motion of the particles and finally make all the particles into the uniformly distribution. In this study， the analytical solution of the diffusion phenomenon of a given particle concentration in a one-dimensional space is calculated by using the method of separation of variables， and the diffusion process is simulated numerically by using the computational physics method combined with Euler's polygonal arc method. The one-dimensional plasma diffusion phenomenon is explained by comparing the theoretical data with the simulated experimental data.

關(guān)鍵詞：一維擴(kuò)散;分離變量法;歐拉法

Key words： one-dimensional diffusion;the method of separation of variables;Euler's polygonal arc method

中圖分類號：O122.2? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?文獻(xiàn)標(biāo)識碼：A? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文章編號：1006-4311（2019）29-0272-04

1? 簡介

受控?zé)岷司圩兪鞘苁廊瞬毮康那把刂匾n題，其目的便是探索清潔可持續(xù)的新能源。其原理是兩個(gè)較輕原子核在融合過程中由于質(zhì)量虧損而釋放出來巨大能量。通常情況下是由一個(gè)氘粒子和一個(gè)氚粒子融合成一個(gè)氦離子并釋放大量能量。在核聚變研究中，等離子體的運(yùn)動(dòng)在很大程度上影響著整個(gè)研究過程，其中可能一個(gè)現(xiàn)象就是擴(kuò)散過程。

擴(kuò)散效應(yīng)是指物質(zhì)分子由于密度分布不均勻而發(fā)生離子轉(zhuǎn)移，直到均勻分布的現(xiàn)象。擴(kuò)散是由于分子、原子等粒子的熱運(yùn)動(dòng)產(chǎn)生的物質(zhì)遷移現(xiàn)象。一般情況下，擴(kuò)散現(xiàn)象都發(fā)生在一種或幾種物質(zhì)于同一物態(tài)或不同物態(tài)之間，是由于不同區(qū)域之間的濃度差或溫度差所引起的[2]。比方說，將兩個(gè)容器中的兩種氣體聯(lián)通。由于分子熱運(yùn)動(dòng)，過一段時(shí)間后這兩種氣體就會(huì)混勻，這是兩種物質(zhì)的擴(kuò)散。換而言之，在空間中存在一種物質(zhì)，它們在初始時(shí)刻都分布在這個(gè)區(qū)域的左側(cè)，由于密度差，這些物質(zhì)就會(huì)向右擴(kuò)散，逐漸充滿整個(gè)區(qū)域，這是一種物質(zhì)的擴(kuò)散，一般發(fā)生在氣體，液體，等離子體中。

本研究主要在于等離子體的擴(kuò)散，我們首先明確等離子體擴(kuò)散的定義。在一個(gè)具體的核聚變裝置中，高溫的等離子體被約束在一個(gè)較小的范圍之中，而由于其在空間中的密度梯度，等離子體之間會(huì)產(chǎn)生碰撞，表現(xiàn)為粒子從密度高的地方像密度低的地方遷移，我們稱其為等離子體的擴(kuò)散現(xiàn)象[1]。

首先，為了求解擴(kuò)散方程，我們引入菲克定律。菲克定律是指在不依靠宏觀的混合作用發(fā)生的傳質(zhì)現(xiàn)象時(shí)，描述分子擴(kuò)散過程中傳質(zhì)通量與濃度梯度之間關(guān)系的定律[1]。菲克定律的表達(dá)式為：

其中q為擴(kuò)散流強(qiáng)度為濃度梯度，D為擴(kuò)散系數(shù)。取空間中x與x+dx，y與y+dx，z與z+dx之間的小平行六面體為代表進(jìn)行研究，平行六面體濃度變化取決于穿過表面的擴(kuò)散流。先考察單位時(shí)間x方向的擴(kuò)散流，在左表面，流量流入平行六面體，在右表面，流量流出。于是我們可以得到：

流入與流出相抵之后，單位時(shí)間內(nèi)x方向凈流入

量為：

因?yàn)榱Ｗ訑?shù)是守恒的，所以如果六面體里邊沒有源或匯，則六面體中單位時(shí)間內(nèi)增加粒子數(shù)應(yīng)等于單位時(shí)間內(nèi)流入的粒子數(shù)，即：

將（4）進(jìn)行化簡，我們就可以得到三維的擴(kuò)散方程：

所以由此可知，一維的擴(kuò)散方程為：

于是，我們推導(dǎo)出了一維擴(kuò)散方程。這個(gè)式子表示的是在一維空間分布下的一些等離子體，將會(huì)有的運(yùn)動(dòng)狀態(tài)。通過確定的初始條件和邊界條件就能夠解出其運(yùn)動(dòng)和分布狀態(tài)。

事實(shí)上，由于初始條件往往很復(fù)雜，常常無法得到精確的理論解，這時(shí)就需要通過數(shù)值計(jì)算才能得到其解。

隨著電子計(jì)算機(jī)的發(fā)展，計(jì)算物理逐步形成的一門新興的邊緣學(xué)科。計(jì)算物理學(xué)以計(jì)算機(jī)為工具，采用數(shù)學(xué)方法解決物理問題，是物理、數(shù)學(xué)和計(jì)算機(jī)三者相結(jié)合的產(chǎn)物[3]。計(jì)算物理方法主要運(yùn)用計(jì)算機(jī)對復(fù)雜的物理問題進(jìn)行數(shù)值計(jì)算或模擬實(shí)驗(yàn)，包括模擬物理過程，研究物理規(guī)律，檢驗(yàn)理論預(yù)測的正確性，核實(shí)實(shí)驗(yàn)數(shù)據(jù)的可靠性等等，探索新的物理規(guī)律[4]。計(jì)算物理有時(shí)被視作理論物理的重要工具，有時(shí)被看做一種“計(jì)算機(jī)實(shí)驗(yàn)”，同時(shí)也有人將其看作介于理論物理與實(shí)驗(yàn)物理之間的第三條物理學(xué)分支。

2? 擴(kuò)散方程解法

2.1 解析解法

我們需要討論特殊情況下的等離子狀態(tài)，假設(shè)在初始時(shí)刻粒子分布在一個(gè)一維的長為l的區(qū)域內(nèi)。我們假設(shè)此區(qū)域一側(cè)為x=0，另一側(cè)為x=l。初始時(shí)刻時(shí)兩端的粒子密度為零，且當(dāng)有粒子擴(kuò)散到邊緣時(shí)，直接從兩側(cè)流失。即：。假設(shè)我們將粒子從中間的位置充入，其密度分布必將體現(xiàn)出中間高兩邊低的狀態(tài)，我們將這個(gè)模型稍作簡化，并使用一個(gè)二次函數(shù)表達(dá)式來進(jìn)行描述，即■。這樣我們就得到了所有的邊界條件，將此兩項(xiàng)邊界條件與（6）相結(jié)合便可以得出理論解析解。

擴(kuò)散方程的解析解：

我們決定采用分離變量法解決這個(gè)問題。在這個(gè)研究之中，濃度是時(shí)間與位置的函數(shù)，所以我們對其進(jìn)行分離變量的處理，得到：

將（7）代入（6）我們可以得到：

移項(xiàng)得：（9）

此時(shí)我們定義等式中左右兩項(xiàng)的值等于？姿。由于我們不知道？姿的正負(fù)，所以我們此時(shí)進(jìn)行分類討論：

①？姿=0：此時(shí)X（x）=0，無意義，舍去。

②？姿<0：我們可以得到，此時(shí)：

將此式帶回到邊界條件中，經(jīng)過計(jì)算可得到：

要想滿足此方程組，唯一的解為：c1=c2=0，沒有意義，舍去。

③？姿<0：此時(shí)我們同樣先考慮X，我們可以得到：

我們將（12）代到邊界條件中，經(jīng)過計(jì)算得到：

所以，很明顯地，我們可以得到：

以下之k皆為正整數(shù)，我們將其代回到（9）可以得到：

經(jīng)計(jì)算可得：

至此，我們就得到了X（x）和T（t），我們將其帶回到（7）中，就可以得到密度分布的公式

在邊界條件中，我們曾列出初始時(shí)刻的狀態(tài)，我們將其代入，得：

則此時(shí)的ck可求，得：

此時(shí)，我們需要討論k的奇偶性帶來的影響，由于這里的k只代表倍數(shù)關(guān)系我們表示奇偶性時(shí)可以用2k+1與2k表示。解得：（20）

所以只有k取奇數(shù)時(shí)是有意義的，我們將ck代回，就可以得到此條件下的一維擴(kuò)散方程：

（21）

雖然已經(jīng)得到解析解法，但是只有在簡單的，有定解條件的情況下，才可以得出解析解。在現(xiàn)實(shí)的使用中，情況會(huì)十分復(fù)雜，三維擴(kuò)散也將會(huì)十分繁瑣。所以就需要引入計(jì)算物理方法，以一維擴(kuò)散為模板，通過對比兩種方法的計(jì)算結(jié)果判斷兩種方法計(jì)算結(jié)果的一致性，并在實(shí)際應(yīng)用上拓展至三維擴(kuò)散的計(jì)算，進(jìn)行應(yīng)用。

2.2 數(shù)值方法

2.2.1 歐拉法的推導(dǎo)

我們已經(jīng)通過分離變量法得到了解析解，接下來我們采用計(jì)算物理方法進(jìn)行計(jì)算并進(jìn)行兩種算法的對比。我們決定采用歐拉折線法，使用軟件計(jì)算并畫圖。

首先，我們要對歐拉折線法進(jìn)行推導(dǎo)。

我們要解決的是常微分方程的數(shù)值解問題，首先我們設(shè)定：

我們假設(shè)方程是連續(xù)的，滿足利普希茨條件[5]，L為利普希茨常數(shù)，隨不同的函數(shù)變化。如此可以確定函數(shù)的連續(xù)性，方程的解y=y（x）是唯一的。在此基礎(chǔ)上我們可以對它做節(jié)點(diǎn)離散，得到節(jié)點(diǎn)離散方程。

首先，我們規(guī)定兩個(gè)節(jié)點(diǎn)之間的間距為步長，滿足。一般情況下h是一個(gè)定值。于是我們就可以得知：

xn=x0+nh，n=0，1，2…（23）

我們可以對■進(jìn)行泰勒展開：

由（22），我們得到：

我們將高階項(xiàng)忽略不計(jì)，將（22）（23）與其整合，可以得到：

移項(xiàng)得：

以上便是歐拉折線法的證明，（27）便是歐拉折線法公式。

我們將歐拉折線法公式與擴(kuò)散方程結(jié)合起來。我們將（6）的左右兩邊同時(shí)用歐拉法進(jìn)行后差分格式的離散化處理。得到：

們將此式化簡，得到：

由此我們得到了遞推公式：

雖然此種歐拉方法計(jì)算相對簡單，但是精度很差，因?yàn)橹挥械谝粋€(gè)點(diǎn)是切線，后面都是首尾相接的折線，如果后面的斜率變化頻繁，變化程度大，此方法的精度就會(huì)差很多，所以我們要對此方法進(jìn)行改進(jìn)。

2.2.2 改進(jìn)歐拉法

改進(jìn)方法：由于歐拉折線法的精度低，主要問題體現(xiàn)在斜率的取值上，所以我們將斜率取值進(jìn)行優(yōu)化。我們將一點(diǎn)的斜率用前后兩點(diǎn)的平均斜率表示。

我們再將其代回到歐拉法中，可以得到改進(jìn)后的歐拉法：

（32）

我們再利用此方法，將（6）代入，進(jìn)行離散化處理，就可以得到改進(jìn)歐拉法的遞推公式。

此公式的形式較為復(fù)雜，但應(yīng)具有更好的精確度。

同時(shí)我們定義一個(gè)穩(wěn)定系數(shù)，在后期使用電腦軟件進(jìn)行編程計(jì)算時(shí)應(yīng)對其進(jìn)行設(shè)定。我們要想讓離散化結(jié)果成立的話，就一定要滿足穩(wěn)定系數(shù)遠(yuǎn)遠(yuǎn)小于一。只有在此情況下，函數(shù)才會(huì)是收斂的，此時(shí)的函數(shù)才會(huì)有一個(gè)趨近值，不會(huì)因積累計(jì)算忽略而帶來巨大的誤差。我們的遞推公式才是有意義的。

2.3 計(jì)算結(jié)果及分析

我們已經(jīng)得到了解析解和數(shù)值解兩種方法的解。我們采用計(jì)算機(jī)軟件，對于計(jì)算結(jié)果進(jìn)行分析。代碼如下：

u=zeros（100，10000）;? ? %存放解的矩陣

s=（1/100000）/（1/100）^2;? ? %穩(wěn)定系數(shù)dt/dx^2

disp（s）;

for i =1：100

u（i，1）=i*（100-i）/100^2;? ? %邊界條件

end

for j=1：10000

u（1，j）=0;

end

for j=1：9999