物理分析课程设计---一维扩散方程求解

课程设计报告

课程名称：核反应堆物理分析题目：一维扩散方程求解院系：核科学与工程学院班级：

学号：

姓名：

指导教师：

成绩：

教师签名：

日期：2011 年6月日

目录

摘要 (1)

课程设计的目的与要求 (1)

设计正文 (1)

课程设计总结或结论 (3)

参考文献 (4)

摘要和关键词

摘要

这个设计用微分方程的差分数值求解方法，运用MATLAB编程计算出一维扩散方程中子通量密度的离散解。

关键词：一维扩散方程

一．课程设计的目的与要求

学习使用微分方程的数值解法（差分方法）来近似求解一维扩散方程，掌握差分方法的核心思想，熟练使用matlab数据处理，origin绘图软件。通过给定的微分方程及边界条件，计算平板型，圆柱形，球形反应堆中子通量密度分布。

二．设计正文

通过查找有关资料，根据二阶线性微分方程

○1

转换为差分方程的一般公式

其中

○2

h为给定步长,

我们把原方程化简为○3

对比方程○1和○3得出

○4

把○4代入○2

等式右端向量

差分方程其实就是一个线性方程组，此线性方程组的系数矩阵为：

则有

这是一个三对角阵，故可用追赶法解式○3。

下面通过matlab程序来计算变换后的差分方程的解。

所编程序如下：

clear;

N=input('请输入参数：');

alpha=input('请输入alpha值：');

if alpha==0

rmax=input('请输入平板的厚度：');

f0=input('请输入平板中心的中子通量密度：');

elseif alpha==1

rmax=input('请输入堆芯半径：');

f0=input('请输入圆柱中心的中子通量密度：');

elseif alpha==2

rmax=input('请输入堆芯半径：');

f0=input('请输入球形中心的中子通量密度：');

end

h=rmax/N;

D=0.8*10^(-2)

for i=1:1:N-1

a(1,i)=2*D*(i-1/2)^alpha*h^(alpha-1);

c(1,i)=2*D*(i+1/2)^alpha*h^(alpha-1);

b(1,i)=a(1,i)+c(1,i)+2*h*8.5*10^(-28)*(i*h)^2;

g(1,i)=2*i^2*h^3*10^14*cos(pi*i*h/2);

end

newa=a(:,2:N-1);

newc=c(:,1:N-2);

Hb=diag(b);

Hc=diag(newc,1);

Ha=diag(newa,-1);

H=-Ha+Hb-Hc;

G=g;

G(1,1)=g(1,1)+a(1,1)*f0;

p(1,1)=b(1,1);

for k=1:1:N-2

q(1,k)=c(1,k)/p(1,k);

p(1,k+1)=b(1,k+1)-a(1,k+1)*q(1,k);

end

for k=1:1:N-2

y(1,1)=G(1,1)/p(1,1);

y(1,k+1)=(G(1,k+1)-a(1,k+1)*y(1,k))/p(1,k+1);

end

for k=N-2:-1:1

u(1,N)=y(1,N-1);

u(1,k+1)=y(1,k)-q(1,k)*u(1,k+2);

end

u(1,1)=f0;

u(1,N+1)=0;

X=0:h:rmax;

P=polyfit(X,u,5)

U=polyval(P,X);

plot(X,U)

三．课程设计总结或结论

本次课程设计加深了我对中子扩散理论的认识，充分的将理论和实践结合起来。基本了解源迭代，中子通量密度迭代的方法，充分认识到有限差分法在扩散临界计算中的重要性。

四．参考文献

微分方程数值方法（第二版）胡建伟汤怀民编MATLAB程序设计与应用（第二版）刘卫国主编

核反应堆物理分析谢仲生等编

附录（设计流程图、程序、表格、数据等）

真实的离散数据

5次多项式拟合所得图形

微分方程在网点