Numerical Solution of Diffusion Equation
刘浩庭 LIU Hao-ting
（北京通州区潞河中学，北京 101199）
（Beijing Luhe High School，Beijing 101199，China）
摘要：扩散现象是其初始密度不均匀分布引起的，它会对物质粒子的分布状态产生影响，最终达到物质在空间均匀分布状态。本研究通过分离变量法对给定条件的粒子浓度在一维空间分布下的扩散现象的解析解进行计算，同时结合使用欧拉法利用计算物理的方法对扩散过程进行了数值模拟。通过对比理论数据与模拟实验数据对等离子体一维扩散现象进行阐释与讨论。
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.
关键词：一维扩散；分离变量法；欧拉法
Key words: one-dimensional diffusion；the method of separation of variables；Euler's polygonal arc method
中图分类号：O122.2 文献标识码：A 文章编号：1006-4311（2019）29-0272-04
1 简介
受控热核聚变是受世人瞩目的前沿重要课题，其目的便是探索清洁可持续的新能源。其原理是两个较轻原子核在融合过程中由于质量亏损而释放出来巨大能量。通常情况下是由一个氘粒子和一个氚粒子融合成一个氦离子并释放大量能量。在核聚变研究中，等离子体的运动在很大程度上影响着整个研究过程，其中可能一个现象就是扩散过程。
扩散效应是指物质分子由于密度分布不均匀而发生离子转移，直到均匀分布的现象。扩散是由于分子、原子等粒子的热运动产生的物质迁移现象。一般情况下，扩散现象都发生在一种或几种物质于同一物态或不同物态之间，是由于不同区域之间的浓度差或温度差所引起的[2]。比方说，将两个容器中的两种气体联通。由于分子热运动，过一段时间后这两种气体就会混匀，这是两种物质的扩散。换而言之，在空间中存在一种物质，它们在初始时刻都分布在这个区域的左侧，由于密度差，这些物质就会向右扩散，逐渐充满整个区域，这是一种物质的扩散，一般发生在气体，液体，等离子体中。
本研究主要在于等离子体的扩散，我们首先明确等离子体扩散的定义。在一个具体的核聚变装置中，高温的等离子体被约束在一个较小的范围之中，而由于其在空间中的密度梯度，等离子体之间会产生碰撞，表现为粒子从密度高的地方像密度低的地方迁移，我们称其为等离子体的扩散现象[1]。
首先，为了求解扩散方程，我们引入菲克定律。菲克定律是指在不依靠宏观的混合作用发生的传质现象时，描述分子扩散过程中传质通量与浓度梯度之间关系的定律[1]。菲克定律的表达式为：
■（1）
其中q为扩散流强度，■为浓度梯度，D为扩散系数。取空间中x与x+dx，y与y+dx，z与z+dx之间的小平行六面体为代表进行研究，平行六面体浓度变化取决于穿过表面的扩散流。先考察单位时间x方向的扩散流，在左表面，流量■流入平行六面体，在右表面，流量■流出。于是我们可以得到：
■（2）
流入与流出相抵之后，单位时间内x方向净流入
量为：
■（3）
因为粒子数是守恒的，所以如果六面体里边没有源或汇，则六面体中单位时间内增加粒子数应等于单位时间内流入的粒子数，即：
■
（4）
将（4）进行化简，我们就可以得到三维的扩散方程：
■（5）
所以由此可知，一维的扩散方程为：
■（6）
于是，我们推导出了一维扩散方程。这个式子表示的是在一维空间分布下的一些等离子体，将会有的运动状态。通过确定的初始条件和边界条件就能够解出其运动和分布状态。
事实上，由于初始条件往往很复杂，常常无法得到精确的理论解，这时就需要通过数值计算才能得到其解。
随着电子计算机的发展，计算物理逐步形成的一门新兴的边缘学科。计算物理学以计算机为工具，采用数学方法解决物理问题，是物理、数学和计算机三者相结合的产物[3]。计算物理方法主要运用计算机对复杂的物理问题进行数值计算或模拟实验，包括模拟物理过程，研究物理规律，检验理论预测的正确性，核实实验数据的可靠性等等，探索新的物理规律[4]。计算物理有时被视作理论物理的重要工具，有时被看做一种“计算机实验”，同时也有人将其看作介于理论物理与实验物理之间的第三条物理学分支。
2 扩散方程解法
2.1 解析解法
我们需要讨论特殊情况下的等离子状态，假设在初始时刻粒子分布在一个一维的长为l的区域内。我们假设此区域一侧为x=0，另一侧为x=l。初始时刻时两端的粒子密度为零，且当有粒子扩散到边缘时，直接从两侧流失。即：■。假设我们将粒子从中间的位置充入，其密度分布必将体现出中间高两边低的状态，我们将这个模型稍作简化，并使用一个二次函数表达式来进行描述，即■。这样我们就得到了所有的边界条件，将此两项边界条件与（6）相结合便可以得出理论解析解。
扩散方程的解析解：
我们决定采用分离变量法解决这个问题。在这个研究之中，浓度是时间与位置的函数，所以我们对其进行分离变量的处理，得到：
■（7）
将（7）代入（6）我们可以得到：
■（8）
移项得：■（9）
此时我们定义等式中左右两项的值等于?姿。由于我们不知道?姿的正负，所以我们此时进行分类讨论：
①?姿=0：此时X（x）=0，无意义，舍去。
②?姿＜0：我们可以得到，此时：
■（10）
将此式带回到边界条件中，经过计算可得到：
■（11）
要想满足此方程组，唯一的解为：c1=c2=0，没有意义，舍去。
③?姿＜0：此时我们同样先考虑X，我们可以得到：
■（12）
我们将（12）代到边界条件中，经过计算得到：
■（13）
所以，很明显地，我们可以得到：
■（14）
以下之k皆为正整数，我们将其代回到（9）可以得到：
■（15）
经计算可得：
■（16）
至此，我们就得到了X（x）和T（t），我们将其带回到（7）中，就可以得到密度分布的公式
■（17）
在边界条件中，我们曾列出初始时刻的状态，我们将其代入，得：
■（18）
则此时的ck可求，得：
■（19）
此时，我们需要讨论k的奇偶性带来的影响，由于这里的k只代表倍数关系我们表示奇偶性时可以用2k+1与2k表示。解得：
■（20）
所以只有k取奇数时是有意义的，我们将ck代回，就可以得到此条件下的一维扩散方程：
■
（21）
虽然已经得到解析解法，但是只有在简单的，有定解条件的情况下，才可以得出解析解。在现实的使用中，情况会十分复杂，三维扩散也将会十分繁琐。所以就需要引入计算物理方法，以一维扩散为模板，通过对比两种方法的计算结果判断两种方法计算结果的一致性，并在实际应用上拓展至三维扩散的计算，进行应用。
2.2 数值方法
2.2.1 欧拉法的推导
我们已经通过分离变量法得到了解析解，接下来我们采用计算物理方法进行计算并进行两种算法的对比。我们决定采用欧拉折线法，使用软件计算并画图。
首先，我们要对欧拉折线法进行推导。
我们要解决的是常微分方程的数值解问题，首先我们设定：
■（22）
我们假设方程是连续的，满足利普希茨条件[5]，■L为利普希茨常数，随不同的函数变化。如此可以确定函数的连续性，方程的解y=y（x）是唯一的。在此基础上我们可以对它做节点离散，得到节点离散方程。
首先，我们规定两个节点之间的间距为步长，满足■。一般情况下h是一个定值。于是我们就可以得知：
xn=x0+nh，n=0，1，2…（23）
我们可以对■进行泰勒展开：
■
■（24）
由（22），我们得到：
■（25）
我们将高阶项忽略不计，将（22）（23）与其整合，可以得到：
■（26）
移项得：
■（27）
以上便是欧拉折线法的证明，（27）便是欧拉折线法公式。
我们将欧拉折线法公式与扩散方程结合起来。我们将（6）的左右两边同时用欧拉法进行后差分格式的离散化处理。得到：
■（28）
我们将此式化简，得到：
■（29）
由此我们得到了递推公式：
■（30）
虽然此种欧拉方法计算相对简单，但是精度很差，因为只有第一个点是切线，后面都是首尾相接的折线，如果后面的斜率变化频繁，变化程度大，此方法的精度就会差很多，所以我们要对此方法进行改进。
2.2.2 改进欧拉法
改进方法：由于欧拉折线法的精度低，主要问题体现在斜率的取值上，所以我们将斜率取值进行优化。我们将一点的斜率用前后两点的平均斜率表示。
■（31）
我们再将其代回到欧拉法中，可以得到改进后的欧拉法：
■
（32）
我们再利用此方法，将（6）代入，进行离散化处理，就可以得到改进欧拉法的递推公式。
■（33）
此公式的形式较为复杂，但应具有更好的精确度。
同时我们定义一个稳定系数■，在后期使用电脑软件进行编程计算时应对其进行设定。我们要想让离散化结果成立的话，就一定要满足稳定系数远远小于一。只有在此情况下，函数才会是收敛的，此时的函数才会有一个趋近值，不会因积累计算忽略而带来巨大的误差。我们的递推公式才是有意义的。
2.3 计算结果及分析
我们已经得到了解析解和数值解两种方法的解。我们采用计算机软件，对于计算结果进行分析。代码如下：
u=zeros(100,10000); %存放解的矩阵
s=(1/100000)/(1/100)^2; %稳定系数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
for i=2:99
u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);
end
end
for j =1:10000
u(100,j)=u(99,j);
end
disp(u);
x=linspace(0,1,100);
y=linspace(0,0.1,100);
a=1;
b=1;
l=1;
[x,t]=meshgrid(x,y);
uu=0;
m=2;
for k=0:m
uu=uu+8*b/pi^3/(2*k^2+1)^3*exp(-((2*k+1)*pi*a/l)^2*t).*sin((2*k+1)*pi/l*x);
end
surf(x,t,uu);
xx=linspace(0,1,100);
x1=linspace(0,1,100);
plot(xx,uu(50,:),'red')
hold on
plot(xx,u(:,5000),'blue')
tt=linspace(0,0.1,10000);
xx1=linspace(0,1,100);
tt1=linspace(0,0.1,100);
subplot(1,2,1)
imagesc(tt,xx1,u)
subplot(1,2,2)
us=uu';
imagesc(tt1,xx1,us)
我们使用软件Matlab对解析解与数值解进行了可视化。
首先我们将两种解法作整体上的对比，图1中颜色越亮代表密度越高。纵轴为空间，横轴为时间。我们发现两幅图的图像是十分相似的，随着时间的推移，粒子从一开始分布在中间，逐渐向两边扩散，并最终溢出此空间。
其次，我们取t=0.5这个时刻，观察此时的密度分布情况（图2），其中，我们发现在此时刻，理论解法与数值解法的结果同样有很高的一致性。
3 结论
我们讨论了在■，初始条件下粒子在一维空间中的扩散问题，通过理论解法，数值解法两种方法进行计算，得出了两种解，并进行了对比。在一维扩散中，我们假设了一种初始状态分布，使其呈二次函数分布，其粒子会逐渐向两侧扩散，过程如图像所示，并最终从两端流失。通过理论解与数值解的对比，我们可以看出这种数值解法与理论解法相差不大，具有一定可行性与准确性。对我们理解物质的扩散过程由一定的帮助。
对于欧拉法，以及改进欧拉法，要对比其准确度，应计算其局部截断误差。局部截断误差[6]即为，假设第n步计算精准，■。欧拉法的局部截断误差为■，而改进欧拉法的局部截断误差为■。说明改进后的欧拉法的误差应当更小。两种方法计算的扩散图像存在一些差异，导致的原因可能是各部计算中的忽略项数导致，亦有可能为欧拉法中折线精度不够导致出现误差。
4 总结与展望
此项目研究得出一维扩散方程，可以使用此方法拓展至二维三维，以此得出某一指定条件下空间中的粒子扩散问题的解，也可以通过更换其他差分格式如有限时域差分进行代码优化，使其更加精确。
参考文献：
[1]郑春开.等离子体物理[M].北京大学出版社，2009.
[2]叶伟国.大学物理[M].清华大学出版社，2012.
[3]顾昌鑫.计算物理学[M].复旦大学出版社，2010.
[4]庞涛（美）.计算物理导论[M].世界图书出版社，2011.
[5]姜玉山.数学物理方程[M].清华大学出版社，2014.
[6]贺俐.计算方法[M].武汉大学出版社，1998.
[7]冯立伟.热传导方程方程几种差分格式的MATLAB数值解法比较[J].沈阳化工大学学报，2011.
[8]史策.热传导方程有限差分法的MATLAB实现[J].咸阳师范学院学报，2009.
[9]田兵.用MATLAB解偏微分方程[J].阴山学刊，2006. |
