计算流体力学大作业
西工大-计算流体力学大作业
计算流体力学大作业学号： 姓名：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）等。
流体力学大作业
“水流动力学基本原理的应用”大作业
姓名：
学号：
专业班级：
成绩：
教师评语：
年月日
1、对水流流向问题有如下一些说法：“水一定从高处向低处流”,“水一定从压强大的地方向压强小的地方流”，“水一定从流速大的地方向流速小的地方流”。
这些说法是否正确？为什么？
2、在写总流能量方程221112221222p p z z h g g g g
ωαυαυρρ++=+++ 时，过水断面上的计
算点、基准面、压强标准是否可以任意选取？为什么？
3、液流通过如图所示管道流入大气中，已知：U 形测压管中水银柱高差m h Hg 2.0=∆，
10.72h m =水柱高，管径m d 1.01=，管嘴出口直径m d 05.02=，不计管中水头损失，试
求：管中流量Q 。
4、如图所示为一水平面上的渐变弯管，已知：
断面1-1处的压强2
3
1/1098m N p ⨯=，流速14/m s υ=，管径mm d 2001=，管径mm d 1002=，转角
45=α，不计弯管的水头损失。
试求：水流作用在弯管上的力
5、如图所示为闸下底板上的消力墩，已知：跃前断面水深h 1=0.6m ，流速v 1=15m/s ，跃后断面水深h 2=4m ，墩宽b=1.6m ，试求：水流对消力墩的作用力。
题2.9图
6、学过恒定总流能量方程及动量方程及其应用这部分内容以后，你觉得有些什么收获？有什么疑惑或者模糊的地方？。
计算流体力学课程大作业
《计算流体力学》课程大作业——基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟张伊哲 航博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 提出问题[问题描述]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.什么叫控制方程？常用的控制方程有哪几个？各用在什么场合？流体流动要受物理守恒定律的支配，基本的守恒定律包括：质量守恒定律、动量守恒定律、能量守恒定律。
如果流动包含有不同组分的混合或相互作用，系统还要遵守组分守恒定律。
如果流动处于湍流状态，系统还要遵守附加的湍流输运方程。
控制方程是这些守恒定律的数学描述。
常用的控制方程有质量守恒方程、动量守恒方程、能量守恒方程、组分质量守恒方程。
质量守恒方程和动量守恒方程任何流动问题都必须满足，能量守恒定律是包含有热交换的流动系统必须满足的基本定律。
组分质量守恒方程，在一个特定的系统中，可能存在质的交换，或者存在多种化学组分，每种组分都需要遵守组分质量守恒定律。
4.研究控制方程通用形式的意义何在？请分析控制方程通用形式中各项的意义。
建立控制方程通用形式是为了便于对各控制方程进行分析，并用同一程序对各控制方程进行求解。
各项依次为瞬态项、对流项、扩散项、源项。
6.CFD商用软件与用户自行设计的CFD程序相比，各有何优势？常用的商用CFD软件有哪些？特点如何？由于CFD的复杂性及计算机软硬件条件的多样性，用户各自的应用程序往往缺乏通用性。
CFD商用软件的特点是功能比较全面、适用性强。
具有比较易用的前后处理系统和其他CAD及CFD软件的接口能力，便于用户快速完成造型、网格划分等工作。
具有比较完备的容错机制和操作界面，稳定性高。
可在多种计算机、多种操作系统，包括并行环境下运行。
常用的商用CFD软件有PHOENICS、CFX、SRAR-CD、FIDAP、FLUENT。
PHOENICS除了通用CFD软件应该拥有的功能外，PHOENICS软件有自己独特的功能：开放性、CAD接口、运动物体功能、多种模型选择、双重算法选择、多模块选择。
计算流体力学大作业报告
课程综合作业课程名称： _________ 计算流体力学 ___________专业班级： _______________ 研究方向：_______________ 学生姓名： ________________ 学号：________________完成日期： _______________________________________计算流体力学课程综合报告1. 简介计算流体动力学（Computational Fluid Dynamics ，简称CFD是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析。
其基本思想为: 把原来在时间域及空间域上连续的物理量的场，如速度场和压力场，用一系列有限个离散点上的变量值的集合来代替，通过一定的原则和方式建立起关于这些离散点上场变量之间关系的代数方程组，然后求解代数方程组获得场变量的近似值。
CFD可以看作是在流动基本方程（质量守恒方程、动量守恒方程、能量守恒方程）控制下对流动的数值模拟。
通过这种数值模拟，我们可以得到极其复杂问题的流场内各个位置上的基本物理量（速度、压力、温度、浓度等）的分布，以及这些物理量随时间的变化情况，确定旋涡分布特性、空化特性及脱流区等。
还可据此算出相关的其他物理星，如旋转式流体机械的转矩、水力损失和效率等。
此外，与CAD联合，还可进行结构优化设计等。
2. 计算流体动学的特点:①流动问题的控制方程一般是非线性的，自变量多，计算域的几何形状和边界条件复杂，很难求得解析解，而用CFD方法则有可能找出满足工程需要的数值解。
②可利用计算机进行各种数值试验，例如，选择不同流动参数进行物理方程中各项有效性和敏感性试验，从而进行方案比较。
③它不受物理模型和实验模型的限制，省钱省时，有较多的灵活性，能给出详细和完整的资料，很容易模拟特殊尺寸、高温、有毒、易燃等真实条件和实验中只能接近而无法达到的理想条件。
④数值解法是一种离散近似的计算方法，依赖于物理上合理、数学上适用、适合于在计算机上进行计算的离散的有限数学模型，且最终结果不能提供任何形式的解析表达式，只是有限个离散点上的数值解，并有一定的计算误差。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
南京理工大学动力工程学院计算流体力学大作业题目基于Fluent的小口径炮弹流体动力学分析专业姓名学号电话成绩年月日基于Fluent的小口径炮弹流体动力学分析摘要小口径火炮武器系统广泛应用于陆军、海军和空军，用于野战防空、要地防空、舰船防空和飞机空中近距格斗。
本文以小口径炮弹为研究对象，对其进行了飞行过程中的流体动力学分析，对其控制方程进行了分析，最后利用ANSYS软件的Fluent模块对其在来流马赫数为2.5，迎角为5度的情况时的空气绕流情况进行了仿真分析，得到了炮弹的阻力系数和升力系数变换图、速度矢量图、流线绕流图和弹的压力分布图，并对所得到的结果进行了分析，得出了一些结论。
这对以后小口径炮弹的改进有很大的帮助。
关键词：小口径火炮仿真 Fluent1、引言小口径速射火炮是抗击中低空飞机、直升机、巡航导弹、战役战术导弹的重要武器装备，是形成弹幕、终端毁伤来袭武器以保卫重要目标的最后一道屏障。
随着战场条件和目标特性的变化，对近程防空反导武器提出了新的需求，在国内外现有小口径速射火炮武器系统的基础上，分析高射速发射火炮武器系统的特点，分析炮弹在出炮口后的飞行流体动力学特性有非常重要的意义。
小口径速射火炮【1】，涵盖23mm、25mm、30mm、35mm、37mm等口径，发射方式涵盖转管发射（多管转管自动机、多转管自动机共架）、转膛发射、双管联动、并行发射及电控串行发射（“金属风暴”）等。
随着技术的进步，小口径速射火炮性能突飞猛进，瞬时射速达到几万～几十万发/min。
其中，射速为1000～8000发/min的小口径火炮发射、弹药技术等技术群称为“高射速发射技术”；而发射速度达到8000发/min以上的小口径火炮发射技术、弹药技术等技术群则称为“超高射速发射技术”。
高射速发射技术，由小口径火炮武器系统的雷达、光电等传感器跟踪来袭目标，计算机解算，指挥火炮，发射密集弹丸形成弹幕，击落穿过中远程防空火力的“漏网者”，有效保卫重要目标、战略要地、机动部队和二次打击能力，是抗击巡航导弹、空地导弹、反舰导弹、制导炸弹以及无人飞机等攻击的有效屏障。
在近年来的历次防空反导作战中都发挥了重要作用。
1973年10月的第4次中东战争，以色列共损失飞机120架，其中被23mm高炮武器系统击落的占55%，在埃军突破“巴列夫”防线的一次战斗中，前3小时以军飞机就被埃军23mm高炮击落18架。
1982年英阿马岛战争中，马岛阿方的地面防空力量只有一个防空营，但它的30mm和40mm高炮却击落了英军11架飞机，占英军损失飞机总数的32%。
1998年12月17日，美英对伊拉克又发动了代号为“沙漠之狐”的空中打击，美军在头两天即发射了305枚巡航导弹，但其中的77枚被23mm高炮武器系统的密集炮火击落实践表明，小口径速射火炮在未来的信息化战争中，大有用武之地，是陆海空三军和二炮的必须装备。
根据小口径速射火炮武器系统不断扩展的应用范围，需要加强应用基础的研究，促进自主创新；需要新的理论支持，推动炮口流场理论发展；开发智能化供输弹、补弹技术，提高持续射击时间；开发大威力、全天候传感器，网络化传感器，提高传感器威力。
未来战争，小口径速射火炮武器系统打击、毁伤目标的种类、性能都有较大发展。
相应地，其预警、跟踪、瞄准、毁伤等性能必须有较大的飞跃，才能适应“建设信息化部队，打赢信息化军队”战略目标的需要。
而对小口径炮弹进行空气动力学分析，对以后的改进有很大的帮助。
2、物理模型本文研究的小口径速射炮弹模型如下图所示，弹头母线为圆弧，船尾角为8度，为了求解方便将炮弹做了如下假设：炮弹是无尾翼的，出炮口后不考虑旋转，由于炮弹的流场是三维的、有粘的，来流速度大部分为超音速的，因此采用理想可压缩气体模型，为了简化计算采用层流模型，控制方程采用无量纲N-S方程，流动为定常流动。
以上工作均在软件中设置完成。
图1 平面图图2 三维图3、控制方程在Fluent 中设置求解器为耦合可压缩模型，这个模型是专门用来求解高速可压缩流体N-S 方程的。
主要用到的控制方程有：质量守恒方程：()0=+∂∂u div tp ρ （1） 式中矢量符号： ()zw y v x u u div ∂∂+∂∂+∂∂=（2） 动量守恒方程【2】： ()()()()()()()()()()()()w v u S w zp z w w y w v x w u t w S v yp z v w y v v x v u t v S u xp z u w y u v x u u t u +∇+∂∂-=∂∂+∂∂+∂∂+∂∂+∇+∂∂-=∂∂+∂∂+∂∂+∂∂+∇+∂∂-=∂∂+∂∂+∂∂+∂∂222υρρρρυρρρρυρρρρ （3）式中符号u S 、v S 、和w S 是动量守恒方程的广义源项。
能量守恒方程：()()T p S gradT c k div uT div t T +⎪⎪⎭⎫ ⎝⎛=+∂∂ρρ （4）式中，p c 是比热容，T 为温度，k 为流体的传热系数，T S 为流体的内热源及由于粘性作用流体机械能转换为热能的部分，grad （）=zy x ∂∂+∂∂+∂∂()()()。
由于是定常流动，所以0()=∂∂t。
4、数值计算方法对于在求解域内所建立的偏微分方程，理论上是有真解的。
但由于所处理的问题自身的复杂性，一般很难得到方程的真解。
因此，需要通过数值方法把计算域内有限数量位置上的因变量当作基本未知量来处理，从而建立方程，然后通过求解代数方程组来得到这些节点值【3】。
有限体积法【4】是计算域划分为一系列的控制体积，将待解微分方程对每一个控制体积积分得出离散方程。
它是目前CFD 应用最广的一种方法。
使用有限体积法建立离散方程，很重要的一步是将控制体积界面上的物理量及其导数通过节点物理量差值求出。
常用的离散格式有中心差分格式、一阶迎风格式、指数格式、二阶迎风格式和QUICK 格式。
本模型用的是二阶迎风格式【5】。
5、边界条件边界条件是指在求解域的边界上所求解的变量或其一阶导数随地点及时间变化的规律。
本文的边界条件为，静压为1个大气压，物面上边界条件为绝热，对速度取无滑移边界条件。
为了保证计算的精度，采用结构网格划分求解区域。
初始条件的给定是不能随意的，不然有可能造成大的偏差，对于本文来说，速度、压强、密度和温度的初始条件给定为常数，即等于来流条件。
轴向流动边界【6】：ii A x r A v dx dR⎪⎪⎭⎫ ⎝⎛+=∞11ϕϕ （5） 横向流动边界：)cos(2θαϕ∞-=∂∂v r（6）6、数值计算结果及分析本文主要研究小口径炮弹在非旋转条件下，不同马赫数下弹丸的空气动力特性，主要分析其外部绕流。
Ma=2.5，攻角为5度是的流场分析由于本文是分析炮弹的外部绕流，所以要建立计算域，建立计算域后，对其进行网格划分，其结果如下图所示。
图3 计算域网格划分图4计算域内部网格网格划分好后，其扭曲度只有0.6，小于0.66，说明网格划分并不好，这是由于电脑内存不足的原因，无法改变，但勉强也还可用。
将网格划分好，并将边界条件设置好后，对其进行求解，求解到120步左右时，计算就达到收敛准则了，如下图所示。
图5通过对阻力系数和升力系数的检测，我们可得到小口径炮弹的阻力系数和升力系数变换图，如下图所示：图6 阻力系数图图7 升力系数图从上图可看出，阻力系数和升力系数在几十步后可使稳定，阻力系数稳定在0.35左右，升力系数稳定在0.29左右。
计算完成后，对其结果进行后处理，本文中使用ANSYS Workbench中的后处理模块进行后处理【7】，得到了XY平面和XZ 平面的速度矢量图，如下图所示。
图8 XY平面速度矢量图图9 XZ平面速度矢量图通过速度矢量图，我们可以看出速度矢量在弹头部有一个明显下降，然后在船尾处有一个明显增强，从而，我们可以知道在弹头处和船尾处形成了激波。
在炮弹底部，可看到速度明显低于弹头和船尾处速度，形成了回流，甚至形成了真空区，从而此处会有压强差，从而会产生底部阻力。
另外，再利用Workbench得到了小口径炮弹的三维流线图和弹体的压力云图和密度等值线图如下图所示。
图10 三维流线图图11 压力云图图12 密度等值线图由上图可看出，弹头部受到的压力较大，所以我们在设计弹体外壳时，要充分考虑到弹头部的受压能力。
从图12可看出，弹头部和船尾部其密度都有较大改变，从此图，我们可推测在这两个位置，形成了激波，从而导致了空气密度有较大的改变。
7、结论本文对小口径炮弹射出后的流场的空气动力特性进行了分析与研究，主要得到了以下结论：1、超声速来流会在锥头部产生一道脱体激波，在椎体与圆柱体转折处一般会产生一道斜激波，在弹尾，空腔的气流与绕炮弹壳体的来流相汇，形成尾迹，弹尾有时会有回流产生。
2、炮弹阻力主要来源于头部所产生的波阻和弹的底部产生的底阻。
3、弹头部受到的压力最大，设计弹体时要考虑弹头的受压能力。
虽然在这次仿真中，我们得到了一些有用的结论，但不得不承认，这里面还有一些不足需要我们改进。
首先，我们的假设有一些不符合实际情况，如认为弹是无旋的，环境为标准大气压下；另外，由于用的电脑为自己的笔记本电脑，在画网格时可能画的不是很好，以至于影响了结果。
最后，感谢任登凤老师的悉心教导，通过您的课程，我学到了很多东西，掌握了流体动力学分析的基础理论，对Fluent软件也有了一定的了解，这对我以后的学习和生活帮助很大。
参考文献[1]朱森元 .小口径速射火炮武器系统发展展望[J].兵工自动化，2008，27（6）：1-8.[2]钱翼稷 .空气动力学[M].北京：北京航空航天大学出版社，2004.9[3]姜根柱 .固冲增程炮弹内外流场统一计算研究[D].南京：南京理工大学硕士论文[4]有限元法[M][5]江春波，张永良，丁则平.计算流体力学[M].中国电力出版社[6]臧国才，李树常.弹箭空气动力学[M].兵器工业出版社，1989.9[7]凌桂龙，丁金滨，温正 .ANSYS Workbench 13.0从入门到精通[M].北京：清华大学出版社，2012.1。