中科大计算流体力学CFD

1、CFD 实验报告姓名：学号、题目： 利用中心差分格式近似导数d2y/dx2，数值求解常微分方程(0 x 1)如=sin2 (0 x 1)dx2x=x=1sin 2T步长分别取 Ax =0.05, 0.01, 0.001, 0.0001。二、报告要求：列出全部计算公式和步骤；表列出程序中各主要符号和数组意义；绘出数值计算结果的函数曲线,并与精确解比较；比较不同差分格式和不同网格步长计算结果的精度和代价附源程序。三、相关差分格式二阶导数d2y/dx2的三点差分格式有向前差分、向后差分和中心差分，表达式分别如下：一阶向前差分： 竺上= t2j+i *j + O(Ax) TOC o 1-5 h z d

2、x2Ax2一阶向后差分：列=.j-1 + j + O (Ax)dx2Ax2二阶中心差分：纽=+厂乞+1 + O (Ax 2)dx2Ax2代入微分方程可以得到差分方程,表达式分别如下：一阶向前差分一阶向后差分二阶中心差分u 2u + u一阶向前差分一阶向后差分二阶中心差分/1j = sin 2 x HYPERLINK l bookmark12 o Current Document Ax2ju 2u + ui2jj = sin 2 x HYPERLINK l bookmark16 o Current Document Ax2ju 2u + ujji = sin 2 x HYPERLINK l bo

3、okmark20 o Current Document Ax2j对于三种差分格式，差分格式可以改写成AY = b的形式，其中A是相同的,非齐次项b不同，如下所示:b=1b=3-211-2 11 -211b=1b=3-211-2 11 -211-2系数矩阵ysin（2 xAx 21sin（2 x） Ax 2sin（2x） Ax2 - yk - 21sin（2x ）Ax2 - ysin（2x ）Ax23sin（2 x） Ax 2sin （2x ）-Ax2 - y I k1sin（2x ） Ax2 - ysin（2x ）Ax22sin（2 x） Ax 2kv 2sin（2x ）Ax2 - yk-11

4、一阶向前差分一阶向后差分二阶中心差分yk - 2人丿。求解AY = byk - 2人丿。12四、计算公式和步骤；关于精确解的推导：已知士 = si 2 x，对X进行两次积分得到y = -4sin2x + C1x + C2，再结合x=1边界条件yx=0 = 0和y“ =1 -晋得到相对应的-和Cx=1关于数值求解方法：对于方程组AY = b可直接求解，也可以使用追赶法求解，下面介绍简单追 赶法求解三对角方程组的过程。对于三对角方程组：dcyb1111adcb22 2 22adcybn1n1n1n1n1adybn n n n系数矩阵A中仅三对角线上的数值不全为0,其余位置上的数值全为0,是典 型的

5、对角占优的三对角矩阵，利用高斯消去法，经过n-1次消元可以化为同解方 程组：1 u11u2 y1 y2q1巳21 uyqn1n1n11ynqn如上所示，求这些值的过程（即消元过程）称为追：u1uiqiu1uiqi=c / d , q = b / d ,1 / 1 1 =c / Id u a 丿=(b q a )/ (dii 1 ii(i = 2,., n 1)u a )i 1 i(i = 2,., n )再利用回代过程求出方程组的各变量：y = q ,(i (i = n 1,n 2,.,1)y = q u yi i i i+1这一逆序求变量的过程（回代过程）称为赶。五、计算结果与分析根据题意需

6、要分别取步长Ax =0.05, 0.01, 0.001, 0.0001进行计算，因此采用MATLAB 进行编程运算，然后将数值解与精确解进行比较，如下图所示：步长 Ax =0.05步 长 Ax =0.01步长 Ax =0.001步长 Ax =0.0001从上面四个图可以看出，这几个格式的精确解和数值解之间符合度很好，其 中中心差分格式精确度更高，并且随着步长的减小与精确解的符合度越高。下面 将计算结果用表格的形式列出：追赶法求解Ax =0.05Ax =0.01Ax =0.001Ax =0.0001计算耗时、/.、八前差0.0068530.0148030.0489311.111667后差0.00

7、66450.0141160.0486621.103086中心差0.0055260.0158000.0596941.072931表1不同格式不同步长计算耗时（单位s）追赶法求解Ax =0.05Ax =0.01Ax =0.001Ax =0.0001计算精度、/.、八前差4.6976e-048.5578e-058.3727e-068.3542e-07后差3.6746e-048.1482e-058.3317e-068.3501e-07中心差6.8119e-085.4443e-105.4441e-135.3257e-16表2 不同格式不同步长计算精度由不同格式不同步长计算耗时和计算精度的结果对比可知：1

8、、从时间代价来看，三种格式采用追赶法求解所需要的时间基本一样，从计算 精度来看，中心差分格式的计算精度明显优于向前差分和向后差分格式。2、随着步长的减小，虽然求解的计算精度增高，但所需要的时间却大大增加。3、对于这三种格式而言，中心差分格式明显优于另两种格式。4、从上可知，提高计算精度可以通过减少步长和采用高阶格式的方法，由于减 少步长需要用时间作为代价，所以采用高阶格式相对而言比较好。六、源程序及其说明符号说明：long变量取值范围dx步长n+1节点数y0始端边界值yl末端边界值A系数矩阵b非齐次项向量Y数值解向量Y e精确解向量T E节点截断误差Q精度源程序如下：clear; %清空lon

9、g=l;%变量取值范围dx=0.05; %设置步长n=long/dx; % 此处n=节点数一1 y0=0;%设定始端边界条件 yl=1-sin(2)/4; %设定末端边界条件A=zeros(n-1,n-1);%设置系数矩阵 AA(1,1:2)=-2 1;%始端边界for i=2:n-2A(i,i-1:i+1)=1,-2,1;endA(n-1,n-2:n-1)=1,-2;% 末端边界b=zeros(n-1,1);%设置矩阵 bfor i=1:n-1%b(i)=dxA2*sin(i-1)*2*dx); %阶向前差分%b(i)=dxA2*sin(i+1)*2*dx); %阶向后差分b(i)=dxA2

10、*sin(i*2*dx); %中心差分endb(1)=b-yO;%设置边界条件b(n-1)=b(n-1)-yl;x=dx*(0:n);%求解精确解Y_e=x-sin(2*x)/4;);disp(SC16OO5O41 李文建);disp( );tic; % 追赶法求解并开始计时a=1; c=-2;u(1)=a/c; q(1)=b(1)/c; for i=2:n-2u(i)=a/(c-u(i-1)*a);endfor i=2:n-1q(i)=(b(i)-q(i-1)/(c-u(i-1)*a);end y(n-1)=q(n-1); for i=n-2:-1:1y(i)=q(i)-u(i)*y(i+1);endY(2:n)=y(l:n-1);