计算流体力学作业电子版
西工大-计算流体力学大作业
计算流体力学大作业学号： 姓名：1、不可压平面流通过二维容器（如图）。
采用 简单迭代、超松弛迭代 求解 势流方程获得容器内的速势和速度分布 。
边界条件按照课本中给，即流经 A 、B 的体积流量为1。
要求： 1）推导差分方程的迭代公式；2）编写计算机程序 ； 3）绘制计算结果曲线 。
答：1）迭代公式推导对于容器中的定常流场，其支配方程为22220x yφφ∂∂+=∂∂ 求解域为下图所示矩形区域则支配方程由有限差分形式代换，得1,,1,,1,,122220()()i j i j i ji j i j i j x y φφφφφφ+-+--+-++=∆∆具有22()()x y ∆+∆的截断误差对于正方形网格，有22()()x y h ∆=∆=，则上式可改写为n=17,1,1,,1,11()4i j i j i j i j i j φφφφφ+-+-=+++若采用简单迭代公式，即Liebmann 公式，则有(1)()(1)()(1),1,1,,1,11()4n n n n n i j i j i j i j i j φφφφφ++++-+-=+++若采用超松弛迭代，即SOR 公式，则有(1)()()(1)()(1),,1,1,,1,1(1)()4n n n n n n i j i j i j i j i j i j ωφωφφφφφ++++-+-=-++++其中松弛因子12ω<<。
ω最佳值opt ω为opt ω=式中cos(/)cos(/)m n αππ=+，m ，n 分别表示在网格系统中垂直线和水平线的总数。
2）计算机程序本程序采用C 语言编写。
程序源代码如下： #include<stdio.h> #include<math.h> void main() { int m=25,n=17,ilast[17],jlast[25]; int step1,step2; double h=0.25; double psi_j[25][17],psiprv_j,vel_j[25][17],velx_j[25][17],vely_j[25][17]; double psi_c[25][17],psiprv_c,vel_c[25][17],velx_c[25][17],vely_c[25][17]; double Pi,Alpha,Omega,Error; int i,j; for(i=0;i<17;i++) jlast[i]=17; for(i=17;i<m;i++) jlast[i]=17-(i-16); for(j=0;j<9;j++) ilast[j]=25; for(j=9;j<n;j++) ilast[j]=25-(j-8); //数据初始化 for(j=0;j<n;j++) { psi_j[0][j]=1.0; psi_c[0][j]=1.0;}for(i=1;i<m;i++){psi_j[i][jlast[i]-1]=1.0;psi_c[i][jlast[i]-1]=1.0; }for(j=0;j<8;j++){psi_j[m-1][j]=1.0;psi_c[m-1][j]=1.0;}for(i=1;i<m-1;i++){if(i>6 && i<21){psi_j[i][0]=0.0;psi_c[i][0]=0.0;}else{psi_j[i][0]=1.0;psi_c[i][0]=1.0;}}for(i=1;i<m-1;i++){for(j=1;j<jlast[i]-1;j++){psi_j[i][j]=0.5;psi_c[i][j]=0.5;}}//处理右上角数据for(i=0;i<m;i++){for(j=0;j<n;j++){if(j>jlast[i]-1){psi_j[i][j]=0;vel_j[i][j]=3;psi_c[i][j]=0;vel_c[i][j]=3;}}}Pi=4.0*atan(1.0);Alpha=cos(Pi/m)+cos(Pi/n);Omega=(8.0-4*sqrt(4-pow(Alpha,2)))/pow(Alpha,2);//计算速势step1=0;step2=0;//简单迭代while(1){Error=0.0;for(i=1;i<m-1;i++){for(j=1;j<jlast[i]-1;j++){psiprv_j=psi_j[i][j];psi_j[i][j]=(psi_j[i-1][j]+psi_j[i+1][j]+psi_j[i][j-1]+psi_j[i][j+1])/4.0;Error=Error+fabs(psi_j[i][j]-psiprv_j);}}step1++;if(step1>1000)break;if(Error<=0.001)break;}//超松弛迭代while(1){Error=0.0;for(i=1;i<m-1;i++){for(j=1;j<jlast[i]-1;j++){psiprv_c=psi_c[i][j];psi_c[i][j]=(1-Omega)*psi_c[i][j]+Omega*(psi_c[i-1][j]+psi_c[i+1][j]+psi_c[i][j-1]+psi_c[i][j+1])/4.0;Error=Error+fabs(psi_c[i][j]-psiprv_c);}}step2++;if(step2>1000)break;if(Error<=0.001)break;}//计算速度for(i=0;i<m;i++){for(j=0;j<jlast[i];j++){if(j==0){vely_j[i][j]=(-3*psi_j[i][j]+4*psi_j[i][j+1]-psi_j[i][j+2])/2/h;vely_c[i][j]=(-3*psi_c[i][j]+4*psi_c[i][j+1]-psi_c[i][j+2])/2/h;}else if(j==jlast[i]-1){vely_j[i][j]=(psi_j[i][j-2]-4*psi_j[i][j-1]+3*psi_j[i][j])/2/h;vely_c[i][j]=(psi_c[i][j-2]-4*psi_c[i][j-1]+3*psi_c[i][j])/2/h;}else{vely_j[i][j]=(psi_j[i][j+1]-psi_j[i][j-1])/2/h;vely_c[i][j]=(psi_c[i][j+1]-psi_c[i][j-1])/2/h;}}}for(j=0;j<n;j++){for(i=0;i<ilast[j];i++){if(i==0){velx_j[i][j]=(-3*psi_j[i][j]+4*psi_j[i+1][j]-psi_j[i+2][j])/2/h;velx_c[i][j]=(-3*psi_c[i][j]+4*psi_c[i+1][j]-psi_c[i+2][j])/2/h;}else if(i==ilast[j]-1){velx_j[i][j]=(psi_j[i-2][j]-4*psi_j[i-1][j]+3*psi_j[i][j])/2/h;velx_c[i][j]=(psi_c[i-2][j]-4*psi_c[i-1][j]+3*psi_c[i][j])/2/h;}else{velx_j[i][j]=(psi_j[i+1][j]-psi_j[i-1][j])/2/h;velx_c[i][j]=(psi_c[i+1][j]-psi_c[i-1][j])/2/h;}}}for(i=0;i<m;i++){for(j=0;j<jlast[i];j++){vel_j[i][j]=sqrt(pow(velx_j[i][j],2)+pow(vely_j[i][j],2));vel_c[i][j]=sqrt(pow(velx_c[i][j],2)+pow(vely_c[i][j],2));}}//输出结果分布FILE *fp;fp=fopen("f:\\ESL\\YFresult.txt","w");fprintf(fp,"简单迭代结果\n");fprintf(fp,"速度势分布\n");for(j=n-1;j>=0;j--){for(i=0;i<ilast[j];i++){fprintf(fp,"%-10.6f\n",psi_j[i][j]);}}fprintf(fp,"速度分布\n");for(j=n-1;j>=0;j--){for(i=0;i<ilast[j];i++){fprintf(fp,"%-10.6f\n",vel_j[i][j]);}}fprintf(fp,"超松弛迭代结果\n");fprintf(fp,"速度势分布\n");for(j=n-1;j>=0;j--){for(i=0;i<ilast[j];i++){fprintf(fp,"%-10.6f\n",psi_c[i][j]);}}fprintf(fp,"速度分布\n");for(j=n-1;j>=0;j--){for(i=0;i<ilast[j];i++){fprintf(fp,"%-10.6f\n",vel_c[i][j]);}}fclose(fp);//输出tecplot数据FILE *fp1;fp1=fopen("f:\\ESL\\TECPLOT-result.txt","w");fprintf(fp1,"title=erwei grid\n");fprintf(fp1,"variables=x, y, psi_easy, velocity_easy, psi_SOR\n, velocity_SOR\n");fprintf(fp1,"zone t=grid,i=25,j=17,f=point\n");for(j=0;j<n;j++){for(i=0;i<m;i++){fprintf(fp1,"%-10.6f,%-10.6f,%-10.6f,%-10.6f,%-10.6f,%-10.6f\n",i*h,j*h,psi_j[i][j],vel_j[i][j],p si_c[i][j],vel_c[i][j]);}}fclose(fp1);}3）计算结果采用简单迭代，容器内的速势和速度分布速势分布（简单迭代）速度分布（简单迭代）采用超松弛迭代，容器内的速势和速度分布速势分布（SOR ） 速度分布（SOR ）2、用点源（汇）分布在对称轴的源汇模拟流体绕过NACA0012旋称体的二维轴对称势流解。
《计算流体力学》作业答案
计算流体力学作业答案问题1：什么是计算流体力学？计算流体力学（Computational Fluid Dynamics，简称CFD）是研究流体力学问题的一种方法，它使用数值方法对流体流动进行数值模拟和计算。
主要包括求解流体运动的方程组，通过空间离散和时间积分等计算方法，得到流体在给定条件下的运动和相应的物理量。
问题2：CFD的应用领域有哪些？CFD的应用领域非常广泛，包括但不限于以下几个方面：1.汽车工业：CFD可以用于汽车流场的模拟和优化，包括空气动力学性能和燃烧过程等。
2.航空航天工业：CFD可以用于飞机、火箭等流体动力学性能的预测和优化，包括机身、机翼的设计和改进等。
3.能源领域：CFD可以用于燃烧、热交换等能源领域的流体力学问题求解和优化。
4.管道流动：CFD可以用于石油、化工等行业的管道流动模拟和流体输送优化。
5.空气净化：CFD可以用于大气污染物的传输和分布模拟，以及空气净化设备的设计和改进。
6.生物医药：CFD可以用于生物流体输送和生物反应过程的模拟和分析，包括血液流动、药物输送等。
问题3：CFD的数值方法有哪些？CFD的数值方法一般包括以下几种：1.有限差分法（Finite Difference Method，FDM）：将模拟区域划分为网格，并在网格上离散化流体运动的方程组，利用有限差分近似求解。
2.有限体积法（Finite Volume Method，FVM）：将模拟区域划分为有限体积单元，通过对流体流量和通量的控制方程进行离散化，求解离散化方程组。
3.有限元法（Finite Element Method，FEM）：将模拟区域划分为有限元网格，通过对流体运动方程进行弱形式的变分推导，将流动问题转化为求解线性方程组。
4.谱方法（Spectral Method）：采用谱方法可以对流体运动方程进行高精度的空间离散，通常基于傅里叶变换或者基函数展开的方式进行求解。
5.计算网格方法（Meshless Methods）：不依赖网格的数值方法，主要包括粒子方法（Particle Methods）、网格自适应方法（Gridless Method）等。
计算流体力学课程大作业
《计算流体力学》课程大作业——基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟张伊哲 航博1011、 引言和综述2、 问题的提出，怎样使用涡量-流函数方法建立差分格式3、 程序说明4、 计算结果和讨论5、 结论1引言虽然不可压缩流动的控制方程从形式上看更为简单，但实际上，目前不可压缩流动的数值方法远远不如可压缩流动的数值方法成熟。
考虑不可压缩流动的N-S 方程：01()P t νρ∇⋅=⎧⎪∂⎨+∇⋅=-∇+∆⎪∂⎩U UUU f U (1.1)其中ν是运动粘性系数，认为是常数。
将方程组写成无量纲的形式：01()Re P t∇⋅=⎧⎪∂⎨+∇⋅=-∇+∆⎪∂⎩U UUU f U (1.2) 其中Re 是雷诺数。
从数学角度看，不可压缩流动的控制方程中不含有密度对时间的偏导数项，方程表现出椭圆-抛物组合型的特点；从物理意义上看，在不可压缩流动中，压力这一物理量的波动具有无穷大的传播速度，它瞬间传遍全场，以使不可压缩条件在任何时间、任何位置满足，这就是椭圆型方程的物理意义。
这就造成不可压缩的N-S 方程不能使用比较成熟的发展型．．．偏微分方程的数值求解理论和方法。
如果将动量方程和连续性方程完全耦合求解，即使使用显示的离散格式，也将会得到一个刚性很强的、庞大的稀疏线性方程组，计算量巨大，更重要的问题是不易收敛。
因此，实际应用中，通常都必须将连续方程和动量方程在一定程度上解耦。
目前，求解不可压缩流动的方法主要有涡量-流函数法，SIMPLE 法及其衍生的改进方法，有限元法，谱方法等，这些方法各有优缺点。
其中涡量-流函数法是解决二维不可压缩流动的有效方法。
作者本学期学习了研究生计算流体课程，为了熟悉计算流体的基本方法，选择使用涡量-流函数法计算不可压缩方腔驱动流问题，并且对于不同雷诺数下的解进行比较和分析，得出一些结论。
本文接下来的内容安排为：第2节提出不可压缩方腔驱动流问题，并分析该问题怎样使用涡量-流函数方法建立差分格式、选择边界条件。
计算流体力学大作业
计算液体力学基础及应用课程期末作业-----程序调试最终版学号：134212059 姓名：徐影ContentsCFD模型示意图一、拟一维喷管理论解求解二、拟一维喷管的CFD求解三、理论值与CFD解的对比CFD模型示意图两圆弧直径为10米，喉部直径为0.59米，长为3米clear all;I=imread('xuying.png'); imshow(I)一、拟一维喷管理论解求解喷管内马赫数的变化公依赖于面积比A/A0，所以可以将Ma作为x的函数1.2.采用隐函数绘图给出理论的马赫数解gamma=1.4;h0=59/100;% 取学生学号后两位数的十分之一作喉部直径syms x Ma A_x y;% xz为x坐标,Ma为马赫数A_x=((10.59-2*sqrt(25-(x-1.5)^2))/0.59)^2;% A_x为面积系数figure('Color',[1 1 1]);set(gcf,'position',[0,0,1.5*468,468]);plot_Ma=A_x^2-(2/(gamma+1)+(gamma-1)/(gamma+1)*y^2)^((gamma+1)/(gamma-1))/y^2;subplot(1,2,1);gca=ezplot(plot_Ma,[0,3]);xlabel('x');ylabel('马赫数');title('采用隐函数求解的马赫数结果');grid on; % 得到两条曲线，由递增规律选取上升曲线段，从该曲线上得到一系列点的坐标为[x0,Ma0]load tk.mat;x_0=tk(:,1);Ma_0=tk(:,2);% 这里load的数据采用某算法从上面出的图取点拟合得到，用到polyval和polyfit函数subplot(1,2,2);plot(x_0,Ma_0);xlabel('x');ylabel('马赫数');title('马赫数的理论解');grid on;求出马赫数后，压力、密度、温度的变化都是Ma的函数，求出理论值并绘图1.2.3.p_0=(1+(gamma-1)/2*Ma_0.^2).^(-gamma/(gamma-1));rho_0=(1+(gamma-1)/2*Ma_0.^2).^(-1/(gamma-1));t_0=(1+(gamma-1)/2*Ma_0.^2).^-1;figure('Color',[1 1 1]);set(gcf,'position',[0,0,1.5*468,1.5*468]);subplot(3,1,1);plot(x_0,p_0);title('压力比理论值');xlabel('x');ylabel('p');grid on; subplot(3,1,2);plot(x_0,rho_0);title('密度比理论值');xlabel('x');ylabel('rho');grid on; subplot(3,1,3);plot(x_0,t_0);title('温度比理论值');xlabel('x');ylabel('T');grid on;二、拟一维喷管的CFD求解clear all;L=3;N=31;dx=L/(N-1);x=linspace(0,L,N);C=0.5;n=2000;student_num=59;A=((10+student_num/100-2*((25-((x-1.5).^2))).^0.5)/(student_num/100)).^2;%面积比A/A_0与x坐标的关系第一步，密度比、温度比、速度比的初始条件设定1.2.3.Rou=1-0.3146*x;rhobi=zeros(1,n);T=1-0.2314*x;V=(0.1+1.09*x).*sqrt(T);P_rou_t=zeros(size(Rou));P_v_t=zeros(size(Rou));P_T_t=zeros(size(Rou));P_rou_t_2=zeros(size(Rou));P_v_t_2=zeros(size(Rou));P_T_t_2=zeros(size(Rou));第二步，预估步第三步，并求Δt,求rou, V, T的预测量1.2.3.第四步，修正步第五步，求平均时间导数1.2.3.最后，得到t+Delta t时刻流动参数的修正值为1.2.3.第七步，边界条件处理for j=1:ntemp=Rou(16);% 第二步，预估步for i=2:30P_rou_t(i)=-V(i)*((Rou(i+1)-Rou(i))/dx)-Rou(i)*((V(i+1)-V(i))/dx)-Rou(i)*V(i)*((log(A(i+1))-log(A(i)))/dx);P_v_t(i)=-V(i)*((V(i+1)-V(i))/dx)-((T(i+1)-T(i))/dx+((Rou(i+1)-Rou(i))/dx)*T(i)/Rou(i))*1/1.4;P_T_t(i)=-V(i)*((T(i+1)-T(i))/dx)-0.4*T(i)*(((V(i+1)-V(i))/dx)+V(i)*((log(A(i+1))-log(A(i)))/dx));end% 第三步，并求Δt,求rou, V, T的预测量dt=C*(dx./(V(2:30)+sqrt(T(2:30))));dt=min(dt);Rou1(2:30)=Rou(2:30)+P_rou_t(2:30).*dt;V1(2:30)=V(2:30)+P_v_t(2:30).*dt;T1(2:30)=T(2:30)+P_T_t(2:30).*dt;V1(1)=V(1);T1(1)=T(1);Rou1(1)=Rou(1);% 第四步，修正步%for i=2:30P_rou_t_2(i)=-V1(i)*((Rou1(i)-Rou1(i-1))/dx)-Rou1(i)*((V1(i)-V1(i-1))/dx)-Rou1(i)*V1(i)*((log(A(i))-log(A(i-1)))/dx); P_v_t_2(i)=-V1(i)*((V1(i)-V1(i-1))/dx)-((T1(i)-T1(i-1))/dx+((Rou1(i)-Rou1(i-1))/dx)*T1(i)/Rou1(i))*1/1.4;P_T_t_2(i)=-V1(i)*((T1(i)-T1(i-1))/dx)-0.4*T1(i)*(((V1(i)-V1(i-1))/dx)+V1(i)*((log(A(i))-log(A(i-1)))/dx));end% 第五步，求平均时间导数P_rou_av=(P_rou_t+P_rou_t_2)/2;P_v_av=(P_v_t+P_v_t_2)/2;P_T_av=(P_T_t+P_T_t_2)/2;% 最后，得到t+Delta t时刻流动参数的修正值为Rou(2:30)=Rou(2:30)+P_rou_av(2:30).*dt;T(2:30)=T(2:30)+P_T_av(2:30).*dt;V(2:30)=V(2:30)+P_v_av(2:30).*dt;P(2:30)=Rou(2:30).*T(2:30);% 第七步，边界条件处理V(1)=2*V(2)-V(3);V(31)=2*V(30)-V(29);Rou(31)=2*Rou(30)-Rou(29);T(31)=2*T(30)-T(29);p=Rou.*T;Ma=V./sqrt(T);rhobi(j)=abs((temp-Rou(16))/temp); % 计算后一次时间步与前一时间步之间的密度比的变化情况，以此检验CFD过程收敛性质end最终结果的绘图figure('Color',[1 1 1]);set(gcf,'position',[0,0,1.2*468,1.5*468]);subplot(3,1,1);plot(1:n,rhobi);xlabel('x');ylabel('Ma');title('相对密度比');grid on;% 密度比收敛情况绘图subplot(3,1,2);plot(x,Ma);title('喷管内马赫数分布');xlabel('x');ylabel('Ma');grid on;% 马赫数CFD值绘图subplot(3,1,3);plot(x,p);title('喷管内压力分布');xlabel('x');ylabel('p');grid on; % 压力分布CFD值绘图shu=[x;A;Ma;V;T;p;Rou];显示各参量最终计算结果fprintf('%6s\t%12s\t%12s\t%12s\t%12s\t%12s\t%12s\r\n','x','A/A_0','Ma','v/v_0','T/T_0','p/p_0','rho')% 依次显示坐标点、形状参数、马赫数、速度、温度、压力的结果fprintf('%6.1f\t%12.4f\t%12.4f\t%12.4f\t%12.4f\t%12.4f\t%12.4f\r\n',shu)x A/A_0 Ma v/v_0 T/T_0 p/p_0 rho0.0 3.1709 0.1859 0.1859 1.0000 1.0000 1.00000.1 2.8156 0.2124 0.2121 0.9975 0.9915 0.99390.2 2.5056 0.2389 0.2383 0.9956 0.9847 0.98900.3 2.2361 0.2711 0.2700 0.9922 0.9728 0.98050.4 2.0030 0.3056 0.3038 0.9885 0.9602 0.97140.5 1.8022 0.3451 0.3422 0.9834 0.9433 0.95910.6 1.6303 0.3882 0.3838 0.9775 0.9234 0.94470.7 1.4844 0.4364 0.4298 0.9700 0.8989 0.92670.8 1.3617 0.4891 0.4794 0.9611 0.8701 0.90540.9 1.2600 0.5469 0.5331 0.9502 0.8362 0.88001.0 1.1771 0.6096 0.5903 0.9374 0.7974 0.85071.1 1.1116 0.6776 0.6508 0.9224 0.7536 0.81701.2 1.0620 0.7507 0.7142 0.9051 0.7053 0.77921.3 1.0273 0.8289 0.7800 0.8855 0.6532 0.73761.4 1.0068 0.9119 0.8475 0.8636 0.5982 0.69271.5 1.0000 0.9998 0.9160 0.8394 0.5416 0.64521.6 1.0068 1.0921 0.9849 0.8132 0.4847 0.59601.7 1.0273 1.1887 1.0534 0.7853 0.4288 0.54611.8 1.0620 1.2893 1.1210 0.7559 0.3753 0.49641.9 1.1116 1.3934 1.1869 0.7255 0.3250 0.44802.0 1.1771 1.5009 1.2507 0.6943 0.2788 0.40152.1 1.2600 1.6113 1.3119 0.6629 0.2371 0.35762.2 1.3617 1.7245 1.3705 0.6315 0.2001 0.31682.3 1.4844 1.8398 1.4258 0.6006 0.1678 0.27952.4 1.6303 1.9576 1.4782 0.5702 0.1400 0.24552.5 1.8022 2.0764 1.5269 0.5408 0.1163 0.21512.6 2.0030 2.1983 1.5732 0.5122 0.0962 0.1879。
《计算流体力学》作业
《计算流体力学》作业西安交通大学航天学院1（第二章）如图所示，在一个二维平行板通道内的中心线上置有一个温度均匀的正方形柱体，计算区域入口流体温度为T in = c ，流速已达充分发展，上下平板绝热，出口边界离开柱体比较远。
试写出层流、稳态对流换热的控制方程组，并对所取定的计算区域写出流速及温度的边界条件。
u (y)T in = c绝热绝热T h进口出口xy2（第四章）推导二阶偏导的四阶精度差分（均分网格）考虑函数在(x, y )=(1,1) 处，(1)计算的精确值；(2)分别采用一阶前差，一阶后差及中心差分近似，其中∆x =∆y= 0.1，计算(1,1)处的估算值。
计算它们与(1)结果的相对误差；(3)重复(2)，只是∆x =∆y= 0.01。
将此时得到的有限差分结果与(2)对比。
(),=+x y x y e e φ/,/∂∂∂∂x y φφ/,/∂∂∂∂x y φφ3（第四章）4（第四章）试证明一维非稳态对流-扩散方程的隐式中心差分（对流、扩散项）格式（1）无条件稳定；（2）是相容的。
数值求解一维线性对流方程：给定初始分布：5（第四章）编程题格式稳定性验证1. 用LAX 格式求解：（1）CFL=0.5时，t=0，t=0.8，1.2时刻的波形；（2）CFL=1.2时，t=0，t=0.8，1.2时刻的波形；2. 用FTCS 格式求解：（1）CFL=0.5时，t=0，t=0.4，t=0.6时刻的波形；（2）CFL=1.2时，t=0，t=0.4时刻的波形；网格数=100网格数=100数值求解一维线性对流方程：给定初始分布：6（第四章）编程题数值粘性影响验证1. 用LAX格式求解：（1）CFL=0.5时，网格数=100，t=0，t=1，2时刻的波形；（2）CFL=0.5时，网格数=100，200，400，800，t=2时刻的波形；有一个二维稳态无源项的对流-扩散问题，已知。
四个边上的Φ值如图所示，其中Δx=Δy=2。
计算流体力学作业习题
2014级西安理工大学计算流体力学作业1.写出通用方程，并说明其如何代表各类守恒定律。
由守恒型对流-扩散方程：()()()div U div T grad S t φφρφρφφ∂+=+∂ 其中φ为通用变量；T φ为广义扩散系数；S φ为广义原项。
若令1;1;0T S φφφ===时，则得到质量守恒方程（mass conservation equation ）()()()()0u v w t x y zρρρρ∂∂∂∂+++=∂∂∂∂ 若令;i u φ=时，则得动量守恒方程（momentum conservation equation ） 以x 方向为例分析，设;u Pu S S x φφ∂==-∂，通用方程可化为：()()()()(2)u uu vu wu P udivU t x y z x x x ρρρρλη∂∂∂∂∂∂∂+++=-++∂∂∂∂∂∂∂z v u u w F y x y z z x ηηρ⎡⎤⎛⎫∂∂∂∂⎡∂∂⎤⎛⎫+++++⎢⎥ ⎪ ⎪⎢⎥∂∂∂∂∂∂⎝⎭⎣⎦⎝⎭⎣⎦同理可证明y 、z 方向的动量守恒方程式 若令;;T pT T S S C φφλφ===时，则得到能量守恒方程(energy conservationequation)()()()()hh div Uh div U div gradT S t ρρρλφ∂+=-+++∂()()()Tp h div Uh div gradT S t C ρλρ∂+=+∂证毕2.用控制体积法离散0)(=+++s dxdT k dx d dx dT u dt dT ，要求对S 线性化，据你的理解，谈谈网格如何划分？交界面传热系数何如何计算？边界条件如何处理？根据守恒型对流-扩散方程： ()()()u T S t x x x ρφρϕφ∂∂∂∂'+=+∂∂∂∂，对一维模型进行分析，则有：0)(=+++s dx dTk dx d dx dT u dt dT将该一维模型的守恒形式在图A 所示的控制容积P 在△t 时间内做积分。
计算流体力学作业电子版
1. 已知有限体积法求解的通用控制方程为()()()u div div grad S tρφρφφ∂+=Γ+∂ 其中φ为通用变量，可代表u 、v 、w 、T 等求解变量。
（1）试说明通用控制方程中各项的物理意义；（2）对于特定的方程，φ、Γ、S 具有特定的形式，对应于质量守恒方程、动量守恒方程、能量守恒方程、（多种化学组分的）组分质量守恒方程，试写出φ、Γ、S 的具体表达式。
解答：（1） 方程中的各项（从左到右）分别是瞬态项、对流项、扩散项及源项。
（2）2. 简述有限体积法建立离散方程时应遵守的四条基本原则。
解答：1) 控制体积界面上的连续性原则当一个面为相邻的两个控制体积所共有时，在这两个控制体积的离散方程中，通过该界面的通量（包括热通量、质量通量、动量通量）的表达式必须相同。
即：通过某特定界面从一个控制体积所流出的热通量，必须等于通过该界面进入相邻控制体积的热通量，否则，总体平衡就得不到满足。
2) 正系数原则中心节点系数aP 和相邻节点系数anb 必须恒为正值。
该原则是求得合理解的重要保证。
当违背这一原则，结果往往是得到物理上不真实的解。
例如，如果相邻节点的系数为负值，就可能出现边界温度的增加引起的相邻节点温度降低。
3) 源项的负斜率线性化原则源项斜率为负可以保证正系数原则。
从式(C2)中看到，当相邻节点的系数皆为正值，但有源项Sp 的存在，中心节点系数aP 仍有可能为负。
当我们规定Sp ≤0，便可以保证aP 为正值。
4) 系数aP 等于相邻节点系数之和原则当源项为0时，我们发现中心节点系数等于相邻节点系数之和，而当有源项存在时也应该保证这一原则，如果不能满足这个条件，可以取SP 为0。
3．什么是对流质量流量F 、扩散传导量D 以及Pelclet 数Pe ，试用定义式表达之。
解答：F 表示通过界面上单位面积的对流质量通量(convective mass flux)，简称对流质量流量；D 表示界面的扩散传导性(量)(diffusion conductance)。
计算流体力学大作业
1 提出问题[问题描述]Sod 激波管问题是典型的一类Riemann 问题。
如图所示，一管道左侧为高温高压气体，右侧为低温低压气体，中间用薄膜隔开。
t=0 时刻，突然撤去薄膜，试分析其他的运动。
Sod 模型问题：在一维激波管的左侧初始分布为：0 ,1 ,1111===u p ρ，右侧分布为：0 ,1.0 ,125.0222===u p ρ，两种状态之间有一隔膜位于5.0=x 处。
隔膜突然去掉，试给出在14.0=t 时刻Euler 方程的准确解，并给出在区间10≤≤x 这一时刻u p , ,ρ的分布图。
2 一维Euler 方程组分析可知，一维激波管流体流动符合一维Euler 方程，具体方程如下： 矢量方程：0U ft x∂∂+=∂∂ (0.1)分量方程：连续性方程、动量方程和能量方程分别是：222,,p u ρ()()()()2000u tx u u pt x x u E p E tx ρρρρ∂⎧∂+=⎪∂∂⎪⎪∂∂∂⎪++=⎨∂∂∂⎪⎪∂+⎡⎤∂⎣⎦+=⎪∂∂⎪⎩ (0.2)其中 22v u E c T ρ⎛⎫=+ ⎪⎝⎭对于完全气体，在量纲为一的形式下，状态方程为：()2p T Ma ργ∞=(0.3)在量纲为一的定义下，定容热容v c 为：()211v c Ma γγ∞=- (0.4)联立(1.2)，(1.3)，(1.4)消去温度T 和定容比热v c ，得到气体压力公式为：()2112p E u γρ⎛⎫=-- ⎪⎝⎭(0.5)上式中γ为气体常数，对于理想气体4.1=γ。
3 Euler 方程组的离散3.1 Jacibian 矩阵特征值的分裂Jacibian 矩阵A 的三个特征值分别是123;;u u c u c λλλ==+=-，依据如下算法将其分裂成正负特征值：()12222k k k λλελ±±+=(0.6)3.2 流通矢量的分裂这里对流通矢量的分裂选用Steger-Warming 分裂法，分裂后的流通矢量为()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭＋＋＋＋＋＋＋＋＋＋＋(0.7)()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭－－－－－－－－－－－(0.8)其中：()()()223321c w γλλγ±±±-+=- c 为量纲为一的声速：22Tc Ma ∞=(0.9)联立(1.3)，(1.9)式，消去来流马赫数得：ργp c =3.3 一阶迎风显示格式离散Euler 方程组 10n n i i x i x i U U f f t xδδ+-++--++=∆∆ (0.10)得到()()n+1nj j 11U =U j j j j t f f f f x++---+∆⎡⎤--+-⎣⎦∆ 算法如下：① 已知初始时刻t=0的速度、压力及密度分布000,,j j j u P ρ，则可得到特征值分裂值0k λ±，从而求出流通矢量0j f ±；② 应用一阶迎风显示格式可以计算出1t t =∆时刻的组合变量1j U ，从而得到1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ；③ 利用1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ可得特征值分裂值1k λ±，从而求出流通矢量1j f ±；④ 按照步骤2的方法即可得到2t t =∆时刻的速度、压力及密度分布222,,j j j u P ρ；⑤ 循环以上过程即可得到()1t n t =+∆时刻的速度、压力及密度分布n+1n+1n+1,,j j j u P ρ。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1. 已知有限体积法求解的通用控制方程为()()()u div div grad S tρφρφφ∂+=Γ+∂ 其中φ为通用变量，可代表u 、v 、w 、T 等求解变量。
（1）试说明通用控制方程中各项的物理意义；（2）对于特定的方程，φ、Γ、S 具有特定的形式，对应于质量守恒方程、动量守恒方程、能量守恒方程、（多种化学组分的）组分质量守恒方程，试写出φ、Γ、S 的具体表达式。
解答：（1） 方程中的各项（从左到右）分别是瞬态项、对流项、扩散项及源项。
（2）2. 简述有限体积法建立离散方程时应遵守的四条基本原则。
解答：1) 控制体积界面上的连续性原则当一个面为相邻的两个控制体积所共有时，在这两个控制体积的离散方程中，通过该界面的通量（包括热通量、质量通量、动量通量）的表达式必须相同。
即：通过某特定界面从一个控制体积所流出的热通量，必须等于通过该界面进入相邻控制体积的热通量，否则，总体平衡就得不到满足。
2) 正系数原则中心节点系数aP 和相邻节点系数anb 必须恒为正值。
该原则是求得合理解的重要保证。
当违背这一原则，结果往往是得到物理上不真实的解。
例如，如果相邻节点的系数为负值，就可能出现边界温度的增加引起的相邻节点温度降低。
3) 源项的负斜率线性化原则源项斜率为负可以保证正系数原则。
从式(C2)中看到，当相邻节点的系数皆为正值，但有源项Sp 的存在，中心节点系数aP 仍有可能为负。
当我们规定Sp ≤0，便可以保证aP 为正值。
4) 系数aP 等于相邻节点系数之和原则当源项为0时，我们发现中心节点系数等于相邻节点系数之和，而当有源项存在时也应该保证这一原则，如果不能满足这个条件，可以取SP 为0。
3．什么是对流质量流量F 、扩散传导量D 以及Pelclet 数Pe ，试用定义式表达之。
解答：F 表示通过界面上单位面积的对流质量通量(convective mass flux)，简称对流质量流量；D 表示界面的扩散传导性(量)(diffusion conductance)。
定义表达式如下：F u D x ρΓ=⎧⎪⎨=⎪δ⎩在此基础上，定义一维单元的Peclet 数Pe 如下：/F uPe D x ρΓ==δ4．二维计算域中节点间距分别为x δ和y δ，控制体积P 的四个界面为e 、w 、n 、5．Pelclet 数Pe 是如何定义的？试根据Pe 的取值范围分析一个流场的流动状况。
解答：当Pe =0，对流-扩散演变为纯扩散问题，即流场中没有流动，只有扩散； 当Pe >0时，流体沿正 x 方向流动；当Pe 数很大时，对流-扩散问题演变为纯对流问题； 当Pe <0时，流体沿负 x 方向流动。
6.什么是假扩散与人工粘性？如何消除或减弱数值计算中假扩散的影响？ 解答：假扩散：对流-扩散方程中，一阶导数项(对流项)离散格式的截断误差，小于二阶而引起的较大的数值计算误差的现象称为假扩散(false diffusion)。
因为低阶离散格式截差的首项包含二阶导数，在数值计算结果中，使扩散的作用被人为地放大了，相当于引入了人工粘性(artificial viscosity)或数值粘性(numerical viscosity)。
消除或减轻数值计算假扩散的方法：采用高阶格式(如二阶迎风格式，QUICK 格式)；采取自适应网格技术。
7.若动量方程和连续性方程分别为()()W P w P E e w w e e D D F F φφφφφφ---=-=-w e F F图1扩散项采用中心差分格式进行离散，试推导二阶迎风格式的对流-扩散方程的离散方程。
解答：1）*当流动沿正方向，即>w u ，>e u （>w F ，>e F ）时，有WW W w φφφ5.05.1-=，WP e φφφ5.05.1-=，此时离散方程(A5)变为（2）*流动沿负方向，即<w u ，<e u （<w F ，<e F ）时，有有E P w φφφ5.05.1-=，EE E e φφφ5.05.1-=，此时离散方程(A5)变为其中，当流动沿正方向0>w F 及0>e F 时，1=α； 当流动沿负方向<w F 及<e F 时，0=α.8.一维稳态问题控制方程为d()d d d d d u S x x x ρφφΓ⎛⎫=+ ⎪⎝⎭，其中广义变量φ（速度、温度、浓度等）在端点A 和B 的边界值为已知，如下图2所示。
采用中心差分格式，试推导控制方程的离散方程及其系数方程组。
图2解答：中心差分格式(Central differencing scheme)：界面物理量利用节点物理量通过线性9.一维稳态问题的控制方程为d()d d d d d u S x x x ρφφΓ⎛⎫=+ ⎪⎝⎭（1）在分段线性插值基础上引入曲率修正（如下图所示）C EP e 812-+=φφφ曲率修正C 为：⎩⎨⎧<+->+-=0,20,2u u C EE E P W P E φφφφφφ试推导此控制方程（1）的QUICK 格式。
解答：当流动沿正方向，即>w u ，>e u （>w F ，>e F ）时，有*流动沿负方向，即<w u ，<e u （<w F ，<e F ）时，得离散方程其中，当流动沿正方向>w F 时，1=w ；>e F 时，1=e α；当流动沿负方向<w F 时，0=w α；10.若变量P φ对时间的积分可表为()[]t f f dt PP tt tP ∆∆01φφφ-+=⎰+ 其中0P φ为t 时刻的值，P φ为t+∆t 时刻的值。
试根据加权因子f 的取值，分析瞬态问题对时间的几种积分方案，并写出相应的表达式。
解答：式中：上标0代表t 时刻，0P φ是t 时刻的值；P φ是t+∆t 时刻的值；f 为0与1之间的加权因子：当f =0时，变量取老值进行时间积分， 当f =1时，变量取新值进行时间积分。
11同812.一维瞬态问题的控制方程为()()S x x x u t +⎪⎭⎫⎝⎛∂∂∂∂=∂∂+∂∂φφρρφΓ 式中φ为广义变量，Γ是相应于φ的广义扩散系数，S 是广义源项。
空间离散采用中心差分格式，时间离散采用全隐式方案，试推导其离散过程，写出离散方程的通用形式以及各项系数的表达式。
解答：假定变量P φ对时间的积分为式中：上标0代表t 时刻，0P φ是t 时刻的值；P φ是t+∆t 时刻的值；f 为0与1之间的加权因子：当f =0时，变量取老值进行时间积分， 当f =1时，变量取新值进行时间积分。
同理，将P φ、E φ、W φ及C P P S S φ+采用类似式(3-5)进行时间积分，式(3-4)可写成P a W a E a 0P +(3-10)(3-11)13.同714.试根据通过界面上单位面积的对流质量流量F 和扩散传导性(量) D 的定义完成下表的填写： 解答：b a a nb nb P P +∑=φφ (1)式中下标nb 表示相邻节点，对于三维问题，相邻节点包括W 、E 、N 、S 、T 、B ；离散方程中各项系数为：⎪⎪⎪⎩⎪⎪⎪⎨⎧∆∆=∆-+∆+∑=∆+=t V a V S a F a a V S a b P P P P nb P C P P 00000ρφ 试给出中心差分格式下离散方程（1）的系数nb a 的表达式，完成下表的填写： 解答：中心差分格式下离散方程的系数nba 的表达式16.简述交错网格与同位网格的区别。
解答：交错网格(staggered grid)：就是将标量(如压力p、温度T和密度等)在正常的网格节点上存储和计算，而将速度的各分量分别在错位后的网格上存储和计算，错位后网格的中心位于原控制体积的界面上。
同位网格（collocated grid）：把速度u、v及压力p同时存储在同一个网格节点上。
即：系统中只有一种类型的控制体积，所有的变量均在此控制体积的中心点处定义和存储，所有控制方程均在此控制体积上进行离散。
17.交错网格有哪些特点？解答：●主控制体积为求解压力p的控制体积[称为标量控制体积(scale controlvolume)或p控制体积]，控制体积的节点P称为主节点或标量节点(scalenode) [图2(a) ]。
速度u在主控制体积的东、西界面e和w上定义和存储，速度v在主控制体积的南、北界面s和n上定义和存储。
u和v各自的控制体积：分别以速度所在位置(界面e和界面n)为中心的，分别称为u控制体积(u-control volume)和v控制体积(v-control volume)●使用交错网格，生成离散方程的方法和过程与原来的基于普通网格的方法和过程一致，只是需要注意所使用的控制体积有所变化。
●交错网格中，由于所有标量(如压力、温度、密度等)仍然在主控制体积上存储，因此，以这些标量为因变量的输运方程的离散过程及离散结果与前面的一样，主要是在交错网格中生成的u和v两个动量方程的离散方程时，积分用的控制体积不再是原来的主控制体积，而是u和v各自的控制体积，同时压力梯度项从源项中分离出来。
18.简述流动控制方程数值分离解法的基本过程。
解答：离解法中应用广泛的是压力修正法，其求解过程为： 对于当前时间步：(1) 假定初始压力场。
(2) 利用压力场求解动量方程，得到速度场。
(3) 利用速度场求解连续方程，得到压力场修正值。
(4) 根据需要，求解湍流方程及其它标量方程。
(5) 判断当前时间步上的计算是否收敛。
若不收敛，返回第二步，迭代计算； 若收敛，重复上述步骤，计算下一时间步的物理量。
压力修正法有很多实现方式，其中，压力耦合方程组的半隐式方法(SIMPLE 算法Semi-Implicit Method for Pressure-Linked Equations ，求解压力耦合方程组的半隐式方法)应用最为广泛，也是各种商用CFD 软件普遍采纳的算法。
19.简述流动控制方程数值耦合解法的基本过程。
解答：● 耦合解法同时求解离散方程组，联立求解出各变量等，其求解过程： 对于当前时间步：(1) 假定初始压力和速度等变量，确定离散方程的系数及常数项等。
(2) 联立求解连续方程、动量方程、能量方程。
(3) 求解湍流方程及其它标量方程。
(4) 判断当前时间步上的计算是否收敛。
若不收敛，返回第二步，迭代计算； 若收敛，重复上述步骤，计算下一时间步的物理量。
20.简述SIMPLE 算法的基本思想和基本步骤。
解答：（1）动量方程的离散 在同位网格上（图1），对于稳态问题动量方程对控制体积P 作积分：()JI J I J i J i nb nb J I J I b A p p u a u a ,,,1,,,+-+=+∑(1)JI A ,为P 点压力作用的面积，实际为控制体积在y 向的高度jj y y -+1.在已知压力场的情况下，Ji p ,和Ji p ,1+可通过线性插值得出。