劉浩庭

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

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ù)的新能源。其原理是兩個較輕原子核在融合過程中由于質(zhì)量虧損而釋放出來巨大能量。通常情況下是由一個氘粒子和一個氚粒子融合成一個氦離子并釋放大量能量。在核聚變研究中，等離子體的運動在很大程度上影響著整個研究過程，其中可能一個現(xiàn)象就是擴散過程。

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

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

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

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

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

量為：

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

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

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

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

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

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

2? 擴散方程解法

2.1 解析解法

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

擴散方程的解析解：

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

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

移項得：（9）

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

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

②？姿<0：我們可以得到，此時：

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

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

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

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

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

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

經(jīng)計算可得：

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

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

則此時的ck可求，得：

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

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

（21）

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

2.2 數(shù)值方法

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

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

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

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

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

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

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

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

由（22），我們得到：

我們將高階項忽略不計，將（22）（23）與其整合，可以得到：

移項得：

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

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

們將此式化簡，得到：

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

雖然此種歐拉方法計算相對簡單，但是精度很差，因為只有第一個點是切線，后面都是首尾相接的折線，如果后面的斜率變化頻繁，變化程度大，此方法的精度就會差很多，所以我們要對此方法進(jìn)行改進(jìn)。

2.2.2 改進(jìn)歐拉法

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

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

（32）

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

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

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

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

我們已經(jīng)得到了解析解和數(shù)值解兩種方法的解。我們采用計算機軟件，對于計算結(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