计算流体力学作业2
西工大-计算流体力学大作业
计算流体力学大作业学号： 姓名：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旋称体的二维轴对称势流解。
计算流体力学作业习题
计算流体⼒学作业习题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 P u 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 ρφρ?φ'+=+，对⼀维模型进⾏分析，则有：k dx d dx dT u dt dT将该⼀维模型的守恒形式在图A 所⽰的控制容积P 在△t 时间内做积分。
图A[]()()()()()et tt tet ttew e w wttwTT TT dx uT uT dt Kdt Sdsdt x x +?+?+-+-=--+（1）⾮稳态项选定T 随x 变化且为阶梯式，既有：()()et t t t t t P P wT T dx T T x +?+?-=-??（2）对流项[]()()()()t tt tew e w tuT uT dt uT uT t+-=-?（3）扩散项()()()()t tt t e w e w tT T T T dt t x x xx +-=-?????()()E P e e T T Tx x δ-?=? ()()P w w T T T w xx δ-?=?（4）原项令S 对t 和x 呈阶梯式变化，既有：t tettwSdsdt S x t+?=综上所述，可以推导出下式：2()()22t t t t tt t t t E wE P w P P u u T T T T T K S t x x φφ+?--+-+=-+由图A 可知，本次⽹格划分采⽤的是外节点法结构化⽹格划分。
流体力学第二章
姓名： 作业成绩： 姓名： 班级： 评阅成绩： 班级： 学号： 学号：流体力学第二章作业1、一矩形平板闸门AB ，门的转轴位于A 端，已知闸门宽3m ，门重9.8kN ，门与水平面夹角030α=，闸门左右水位分别为11m h =，23 1.73m h h ==，若不计门轴的摩擦，在门的B 点用钢索起吊，试求闸门启动时所需的拉力F 。
2、如图示，一直角形闸门，已知闸门的宽度1m B =， 1m h =，试求关闭闸门受需的垂直作用力。
3、圆柱形锅炉直径为3m D =，右端为焊接的盖板，水深为1.5m ，锅炉内气体压力可通过U 型管水银压差计来量测。
若测得压差计两只管内水银面高差为0.6m h ∆=，试求作用在锅炉盖板上的静水总压力大小。
（注：①水银密度3313.610kg/m Hg ρ=⨯；②水面上压强不为零）4、某竖直隔板上开有矩形孔口，如图所示，高1m a =，宽3m b =。
现用一直径2m d =的圆柱筒将其堵塞。
隔板两侧充水，左侧2m h =，水位高差0.6m z =。
求作用于该圆柱筒的静水总压力。
注：（1）计算水平压力时左侧投影面积的抵消关系；（2）计算铅垂压力时的压力体（浮力定律）5、如图所示，闸门AB 宽 1.2m b =，A 点为铰链连接，闸门可以绕A 点转动。
左边容器压力表G 的读数为014.7kPa p =-，右边容器中油的密度为'3850kg/m ρ=，问在B 点作用多大的水平力F 才能使闸门AB 保持平衡。
（16分）注：液面存在气压0p 时，压力中心公式为：0sin sin xC D C C g I y y p A gy A ραρα=++6、如图示，在一蓄水容器垂直壁的下面有一1/4圆柱面部件AB ，其长度为1m l =，半径为0.4m R =，容器中水深为 1.2m h =，试求作用在部件AB 上的总压力。
（请大家作图，并在图上标注压力体）答案是按l=1.5m 做的。
计算流体力学课后题作业
课后习题第一章1.计算流体动力学的基本任务是什么计算流体动力学是通过计算机数值计算和图像显示，对包含有流体流动和热传导等相关物理现象的系统所做的分析。
2.什么叫控制方程？常用的控制方程有哪几个？各用在什么场合？流体流动要受物理守恒定律的支配，基本的守恒定律包括：质量守恒定律、动量守恒定律、能量守恒定律。
如果流动包含有不同组分的混合或相互作用，系统还要遵守组分守恒定律。
如果流动处于湍流状态，系统还要遵守附加的湍流输运方程。
控制方程是这些守恒定律的数学描述。
常用的控制方程有质量守恒方程、动量守恒方程、能量守恒方程、组分质量守恒方程。
质量守恒方程和动量守恒方程任何流动问题都必须满足，能量守恒定律是包含有热交换的流动系统必须满足的基本定律组分质量守恒方程，在一个特定的系统中，可能存在质的交换，或者存在多种化学组分，每种组分都需要遵守组分质量守恒定律。
4.研究控制方程通用形式的意义何在？请分析控制方程通用形式中各项的意义。
建立控制方程通用形式是为了便于对各控制方程进行分析，并用同一程序对各控制方程进行求解。
各项依次为瞬态项、对流项、扩散项、源项。
6.CFD商用软件与用户自行设计的CFD程序相比，各有何优势？常用的商用CFD软件有哪些？特点如何？由于CFD的复杂性及计算机软硬件条件的多样性，用户各自的应用程序往往缺乏通用性。
CFD商用软件的特点是功能比较全面、适用性强。
具有比较易用的前后处理系统和其他CAD及CFD软件的接口能力，便于用户快速完成造型、网格划分等工作。
具有比较完备的容错机制和操作界面，稳定性高。
可在多种计算机、多种操作系统，包括并行环境下运行。
常用的商用CFD软件有PHOENICS、CFX、SRARCD、FIDAP、FLUENT。
PHOENICS除了通用CFD软件应该拥有的功能外，PHOENICS软件有自己独特的功能：开放性、CAD接口、运动物体功能、多种模型选择、双重算法选择、多模块选择。
流体力学作业2
第二章 流体的平衡2.2 给出如下体力场，分别在均质或正压流体斜压流体情况下说明流场能否静止：()()()();1222222k y xy x j x zx z i z yz y f++++++++=()k r r kf,23-=为常数，r 是内径。
[解答] （1）流体平衡的基本方程为P F b ∇=ρ对于正压流体（含均质流体）平衡的必要条件是体力有势，即体力无旋0=⨯∇b F 。
对于斜压流体在有势力作用下不可能处于平衡态。
()()()k y x j x z i z y y xy x x zx z zyz y z y x kji F b -+-+-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++++++∂∂∂∂∂∂=⨯∇222222222则当z y x ==时，正压流体平衡；当不满足此条件时，斜压流体才有可能平衡。
（2）0003=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-∂∂∂∂∂∂=⨯∇r k r k ji F b ϑθ，∴正压流体平衡，斜压流体不能平衡。
2.3 图示一种酒精和水银的双液测压计，细管上端为大气压时酒精液面高度为零，当细管上端的表压为p 时，酒精的液面下降h,试用321,,d d d 和h 来表示p 。
酒精和水银密度均为已知。
[解答] 设开始时酒精液面高度（即液面上端距酒精水银的分界面）为H ，则有流体静力学规律知，分界面即等压面，从而有H g P gH P atm atm ∆+=+mercury alcohol ρρ（1）施加表压P 后，同理可得 )()(21mercury 1alcohol h h H g P h h H g P P atm atm ++∆+=+-++ρρ（2）其中1h 为分界面下降高度，2h 为水银面上升高度，由液体排开体积相等得：()2222312221444h d d h d h d -==πππ（3）解之得h d d d h h d d h 222321222211,-==（4）()(),12-并将（4）代入得)()(21mercury 1alcohol h h g h h g P ++-=ρρ)11()1(22232221mercury 2221alcohol d d d ghd d d gh -++-=ρρ 2.4 如图的微测压计用来测量两容器E 和B 中的气体压强差，试用21,,,ρρδd 表示,B E p p -并说明当横截面积,A a 〈〈，而且两种液体1ρ和2ρ相近似，很小的B E p p -就能引起很大的液面高差d ，从而提高了测量精度。
计算流体力学课程大作业
《计算流体力学》课程大作业——基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟张伊哲 航博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节提出不可压缩方腔驱动流问题，并分析该问题怎样使用涡量-流函数方法建立差分格式、选择边界条件。
重庆大学高等流体作业流体作业2(DOC)
研究生课程考核试卷科目：高等流体力学教师：何川姓名：苗闪闪学号：20111002060专业：动力工程及工程热物理类别：学术考生成绩：卷面成绩平时成绩课程综合成绩阅卷评语：阅卷教师（签名）重庆大学研究生院制四、设的气体以的速度以零攻角的速度定常绕流长度为的大平板，试用数值解讨论边界层内的流动规律。
解：设以匀速U运动的流体沿切线方向绕过一静止平板，由于流体的粘性作用，板面上的流体流速被降为零，板壁附近极薄区域的流体沿板的法向存在很大的速度梯度，该区域被称为流动边界层。
，为大雷诺数绕流。
沿板面取x坐标，板的发现方向取y坐标舍边界层的厚度为δ=δ(x)，因厚度很薄，既有。
YV ∞x根据问题的特性，忽略质量力，边界层内的速度矢量,压强，流体力学基本方程具体化为：C.E. (4.1) M.E. (4.2)(4.3)下面，应数量级比较的方法分析方程中各项的数量级关系。
首先，选特征参数U、L，奖方程组中个变物无量纲化，并在无量纲化过程中注意使用各无量纲变量的取值范围为(0，1)或取1的数量级（记作(~1)）。
即令：注：大雷诺数流动时，压力采用流动造成的动压作为参照具有相当的数量级关系。
将各无量纲变量代入式（4.1），并比较其各项的数量级关系，有：将各无量纲变量代入式（4.2），并比较其各项的数量级关系，有：1将各无量纲变量代入式（4.3），并比较其各项的数量级关系，有：1根据物理依据和工程事实，在连续型方程中，认为A=B；在x方向的大量方程中，去掉数量级较低的第一项，保留第二项，即边界层中的粘性作用，认为该项于压力梯度项具有相同的数量级，则有；分析y方向的动量方程，对流项与粘性扩散项具有相同的数量级，而压力梯度项却具有高得多的数量级。
得到简化的近似方程组：C.E. (4.4) M.E. (4.5)(4.6)由式（4.6）可以认为，在任一过断流面上，边界层内各点的压力与其外边界上的流势压力p e相等，即，而边界层外势流区满足这里，即有边界层方程即可表示为C.E. (4.4) M.E. (4.7)对于平板绕流，V e=U=c，则平板然刘的边界层发成可以简化为C.E. (4.4) M.E. (4.8)其定界边界条件为：y=0 : u=0 , v=0(4.9a) y→∞: u=U(4.9b)下面采用无量纲相似性解法。
流体力学第二章作业答案
2.3 如图，用U 型水银测压计测量水容器中某点压强，已知H 1=6cm ，H 2=4cm ，求A 点的压强。
解：选择水和水银的分界面作为等压面得11222()γγ++=+a A p H H p H故A 点压强为：511212() 1.1410Pa γγγ=++-=⨯A a p p H H2.5 水压机是由两个尺寸不同而彼此连通的，以及置于缸筒内的一对活塞组成，缸内充满水或油，如图示：已知大小活塞的面积分别为A 2，A 1，若忽略两活塞的质量及其与圆筒摩阻的影响，当小活塞加力F 1时，求大活塞所产生的力F 2。
帕斯卡定律：加在密闭液体上的压强，能够大小不变地由液体向各个方向传递。
根据静压力基本方程(p=p 0+ρgh),盛放在密闭容器内的液体，其外加压强p 0发生变化时，只要液体仍保持其原来的静止状态不变，液体中任一点的压强均将发生同样大小的变化。
这就是说，在密闭容器内，施加于静止液体上的压强将以等值同时传到各点。
这就是帕斯卡原理，或称静压传递原理。
解：由111F p A =，222Fp A =，根据静压传递原理：12p p =1221F A F A ⇒=2.6如图示高H =1m 的容器中，上半装油下半装水，油上部真空表读数p 1=4500Pa ，水下部压力表读数p 2=4500Pa ，试求油的密度ρ。
解：由题意可得abs1a 1p p p =-，abs2a 2p p p =+abs1abs222H H p gp ργ++= 解得abs2abs1213()()22836.7kg/m 22a a H Hp p p p p p gH gH γγρ--+---===2.7 用两个水银测压计连接到水管中心线上，左边测压计中交界面在中心A 点之下的距离为Z ，其水银柱高度为h 。
右边测压计中交界面在中心A 点之下的距离为Z+∆Z ，其水银柱高为h+∆h 。
（1）试求∆h 与∆Z 的关系。
（2）如果令水银的相对密度为13.6，∆Z=136cm 时，求∆h 是多少？解：（1）分别取测压计中交界面为等压面得，a 12AA 2a 1()()p h z p p z z p h h γγγγ+=+⎧⎨++∆=++∆⎩ 解得∆h 与∆Z 的关系为：h z ∆=∆12γγ （2）当∆Z=136cm 时，cm 1012=∆=∆γγzh 2.9 如图示一铅直矩形平板AB 如图2所示，板宽为1.5米，板高h =2.0米，板顶水深h 1=1米，求板所受的总压力的大小及力的作用点。
计算流体力学大作业
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 ρ。
流体力学作业2答案
作业2答案（第3章、第4章）第3章一、选择题1、流体运动的连续性方程是根据（C）原理导出的。
A、动量守恒B、质量守恒C、能量守恒D、力的平衡2、流线和迹线重合的条件为（ C ）A、恒定流B、非恒定流C、非恒定均匀流二、判断题1、以每个流体质点运动规律为研究对象的方法称为拉格朗日法。
（正确）2、恒定流一定是均匀流。
（错误）3、涡流是指流体质点在运动中不绕自身轴旋转的流动。
（正确）4、无旋流就是无涡流。
（正确）5、非均匀流一定是非恒定流。
（错误）三、简答题1、述流体运动的两种方法是什么？简述其内容。
答：研究流体运动有两种不同的观点，因而形成两种不同的方法：一种方法是从分析流体各个质点的运动着手，即跟踪流体质点的方法来研究整个流体的运动，称之为拉格朗日法；另一种方法则是从分析流体所占据的空间中各固定点处的流体的运动着手，即设立观察站的方法来研究流体在整个空间里的运动，称其为欧拉法2. 流体微团体运动分为哪几种形式?答：①平移②线变形③角变形④旋转变形。
3. 写出恒定平面势流中流函数、势函数与流速的关系。
（改为：写出恒定平面势流中流函数具有的性质，流函数与流速势的关系。
）答：流函数具有的性质（1）流函数相等的点组成的线即流线，或曰，同一流线上个点的流函数为常数。
（2）两流线间的函数值之差为为单宽流量。
（3）平面势流的流函数为一个调和函数。
答：流函数与流速势的关系（1）流函数与势函数为调和函数。
（2）等势线与等流函数线正交。
4．什么是过流断面和断面平均流速？为什么要引入断面平均流速？答：与流线正交的断面叫过流断面。
过流断面上点流速的平均值为断面平均流速。
引入断面平均流速的概念是为了在工程应用中简化计算。
5．如图所示，水流通过由两段等截面及一段变截面组成的管道，试问：(1)当阀门开度一定，上游水位保持不变，各段管中，是恒定流还是非恒定流？是均匀流还是非均匀流？(2)当阀门开度一定，上游水位随时间下降，这时管中是恒定流还是非恒定流？(3)恒定流情况下，当判别第II 段管中是渐变流还是急变流时，与该段管长有无关系？答：（1）是恒定流。
流体力学作业(2)
流体力学作业（2）
1、在封闭管端完全真空的情况下，水银柱差Z2=50mm，求盛水容器液面绝对压强p1和水面高度Z1。
（水银重度133375N/m3，水重度9807 N/m3）（p1 =6668.75Pa，Z1=0.68m）
2、封闭容器水面的绝对压强p0 =107.7KPa，当地大气压强p a =98.07Kpa，试求（1）水深h1=0.8m 时，A点的绝对压强和相对压强。
；（2）压力表M和酒精（重度7.944KN/m3）测压计h的读数。
（p Aabs=115545.6 Pa，p A=17475.6 Pa，p M=9630 Pa，h=1.212m）
3、已知水深h=1.2m，水银柱高度h p=240mm，大气压强p a =730mmHg，连接橡皮软管中全部是空气，求封闭水箱水面的绝对压强及真空度。
（p abs=53585.35Pa，p v=43778.4 Pa）
4、已知图中Z=1m，h=2m，求A的相对压强及测压管中液面气体压强的真空度。
（p A= - 9807Pa，p v=19614 Pa）
5、已知水箱真空表M的读数为0.98 KN/m2，水箱与油箱的液面差H=0.15m，水银柱差h2=0.2m，油的重度为7.85 KN/m3，求h1。
（h1=12.377m）。
流体力学课后作业答案
流 体 力 学
解： h p p0 117.7 98.07 2mH O 2
g
9800
P pC A ghC A 9800 (2.8 0.3) 0.5 0.6 9.114kN
ye I Cx yC A (1/12) 0.5 (0.6) 3.1 0.5 0.6
流 2-7 测压管中水银柱差Δh=100mm，在水深h=2.5m处安 体 装测压表M，求其读数，并图示测压管水头线的位置。 力 学
测压管水头线 p0 h’ h M
解： p p0 w gh
Hg g h w gh (13.6 0.1 1 2.5) 9.8 37.83kPa
1 1 1C 1
油 h1 水 h2 θ
800 9.8 0.5 (1/ sin 60) 1 4.52kN
1 gh1 2 gh1
'
800 1 1000 h1
'
F2 2 gh2C A2
h1 0.8m
'
1000 9.8 (0.8 1) (2 / sin 60) 1 40.74kN
l
hf
13.6 0.92
2
h 1.24m
又 h l v f
0.92
0.2
d 2g 设为层流 Re 64 320 2000 成立
h
vd Re
1 0.025 320
7.8 10 m /s
2
5
若反向流动，Q不变，Re不变，λ不变，hf不变， 所以h不变，只是反向高差为9cm。
P Px Pz 121.85kN
CFD-作业题
试利用Fourier方法，分析其稳定性
4.1 构造高分辨率差分格式，并进行理论分析及数值实验
针对单波方程： u u 0
t x
对于空间导数，构造出一种不超过6点格式；并进行Fourier误差分析， 画出kr,ki的曲线。 形如： u 要求：精度不限； u a u a u a u a u a u x 网格基架点数不超过6个； u u a u a u a u a u a u x 能够分辨的波数范围尽量宽； …… （即kr,ki曲线近可能接近准确解） 给出差分的具体表达式， 画出kr,ki的曲线； 说明构造格式的阶数，并采用本PPT第5页的方法给出的精度验证；
如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。 我们要求其有4阶精度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。 6个未知数（a1-a6）， 5个方程； 有1个自由参数。 调整这个自由参数，使得kr,ki曲 线最为理想。 如何调整？ 1) 可以人工调整，观察kr,ki曲线，选取满意的。 2）可自动调整，设立一个优化目标函数。 例如 * : * ki ( * ) 0.05 调整自由参数，使得该目标函数取最大值。 思路：牺牲精度，提高分辨率
t 0:
x [0,1]
x 0.5 0,1,1 (u, , p) x 0.5 0,0.125,0.1
计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进行比较 （要求画在一张图上）。 要求： 1） 空间网格数100， 时间推进格式选用3阶Runge-Kutta,时间步长自选。 2） 可选用逐点分裂，也可选用特征分裂。 3） 建议采用本讲作业题2（或作业题1）自行构造的差分格式计算。 （作业题2是激波捕捉格式，效果应当会好些）。 如果作业题1和作业题2遇到困难，也可采用现有的差分格式。
国开电大流体力学形考作业2参考答案
A.该点的相对压强为负值
B.该点的绝ቤተ መጻሕፍቲ ባይዱ压强为正值
C.该点的相对压强为正值
D.该点的绝对压强为负值
【答案】：该点的相对压强为负值；该点的绝对压强为正值
题目2.2．流体静压强的（）。
A.大小与受压面积方位无关
B.方向与受压面有关
C.大小与受压面积有关
【答案】：大小与受压面积方位无关
【答案】：
解：用水柱高表示
1.该点绝对压强：8.16mH2o
2.该点相对压强：-1.84mH2o
3.该点真空压强：1.84mH2o
用水银柱表示
1.该点绝对压强：599.1 mm Hg
2.该点相对压强：-135.4 mm Hg
3.该点真空压强：135.4 mm Hg
题目3.1．等压面与质量力垂直。
对
错
【答案】：对
题目4.2．某点存在真空，是指该点的压强小于大气压强。
对
错
【答案】：错
题目5.1．等压面是水平面的条件是什么？
【答案】：
1.连续介质2.同一介质3.单一重力的作用下
题目6.1．已知某点绝对压强为80kN/m2，当地大气压强pa=98kN/m2。试将该点绝对压强、相对压强和真空压强用水柱及水银柱表示。
计算流体力学作业
计算流体力学课程作业任课教师：***姓名：学号：指导老师：目录1.写出通用方程，并说明如何代表各类守恒方程。
(1)2.推导流体运动的质量、动量守恒方程。
(2)3.简述源项线性化、网格划分问题。
(5)4.用ddx (KðTðx)+S=0，谈谈边界条件如何处理。
(7)5.用有限体积法离散ρc ðTðt=ððx(KðTðx)，并推广到二维、三维问题，写出过程。
(8)6.从不同角度对流体运动分类。
(11)7.谈谈物理模型试验与计算流体力学方法的关系。
(12)8.讨论离散对流项时离散格式的进化过程。
(13)9.利用幂函数格式离散二维、三维通用方程的离散方程。
(14)10.解释交错网格的概念。
(15)11.简述压力校正法解N-S方程的过程。
(15)12.思考∑a nb v nb′为什么可以省去。
(17)1.写出通用方程，并说明如何代表各类守恒方程。
答：(1)写出通用方程。
在Cartesian坐标系下单位体积黏性流动N-S方程组微分形式如下：{ðρðt+∇∙(ρV)=0 (1)ðρuðt+∇∙(ρuV)=∇∙(μ∇u)+13μ[ððx(∇∙V)]−ðpðx+F x+S mx(2a)ðρvðt+∇∙(ρvV)=∇∙(μ∇v)+13μ[ððy(∇∙V)]−ðpðy+F y+S my(2b)ðρw ðt +∇∙(ρwV)=∇∙(μ∇w)+13μ[ððz(∇∙V)]−ðpðz+F z+S mz(2c)ðρeðt+∇∙(ρeV)=∇∙(k∇T)−p(∇∙V)+Φ+Q (3)上述微分形式黏性流动N-S方程组中，式(1)为连续性方程，式(2a)、(2b)、(2c)分别为x、y、z方向上的动量方程，式(3)为能量方程。
计算流体力学大作业
南京理工大学动力工程学院计算流体力学大作业题目基于Fluent的小口径炮弹流体动力学分析专业姓名学号电话成绩年月日基于Fluent的小口径炮弹流体动力学分析摘要小口径火炮武器系统广泛应用于陆军、海军和空军，用于野战防空、要地防空、舰船防空和飞机空中近距格斗。
本文以小口径炮弹为研究对象，对其进行了飞行过程中的流体动力学分析，对其控制方程进行了分析，最后利用ANSYS软件的Fluent模块对其在来流马赫数为2.5，迎角为5度的情况时的空气绕流情况进行了仿真分析，得到了炮弹的阻力系数和升力系数变换图、速度矢量图、流线绕流图和弹的压力分布图，并对所得到的结果进行了分析，得出了一些结论。
这对以后小口径炮弹的改进有很大的帮助。
关键词：小口径火炮仿真 Fluent1、引言小口径速射火炮是抗击中低空飞机、直升机、巡航导弹、战役战术导弹的重要武器装备，是形成弹幕、终端毁伤来袭武器以保卫重要目标的最后一道屏障。
随着战场条件和目标特性的变化，对近程防空反导武器提出了新的需求，在国内外现有小口径速射火炮武器系统的基础上，分析高射速发射火炮武器系统的特点，分析炮弹在出炮口后的飞行流体动力学特性有非常重要的意义。
小口径速射火炮【1】，涵盖23mm、25mm、30mm、35mm、37mm等口径，发射方式涵盖转管发射（多管转管自动机、多转管自动机共架）、转膛发射、双管联动、并行发射及电控串行发射（“金属风暴”）等。
随着技术的进步，小口径速射火炮性能突飞猛进，瞬时射速达到几万～几十万发/min。
其中，射速为1000～8000发/min的小口径火炮发射、弹药技术等技术群称为“高射速发射技术”；而发射速度达到8000发/min以上的小口径火炮发射技术、弹药技术等技术群则称为“超高射速发射技术”。
高射速发射技术，由小口径火炮武器系统的雷达、光电等传感器跟踪来袭目标，计算机解算，指挥火炮，发射密集弹丸形成弹幕，击落穿过中远程防空火力的“漏网者”，有效保卫重要目标、战略要地、机动部队和二次打击能力，是抗击巡航导弹、空地导弹、反舰导弹、制导炸弹以及无人飞机等攻击的有效屏障。
流体运动学II
u(x x, y y, z z,t) v(x x, y y, z z,t) w(x x, y y, z z,t)
• 以x方向速度为例，由泰勒级数得，
u(x x, y y, z z,t)
u(x, y, z,t) u x u y u z （18）
• 角变形速率：过任意点作正交微元流体面，每个面上任作 两条过Ｏ点的正交流体线，在时间间隔 t 前后有两条相
应的角平分线，则每条流体线与这两条角平分线的夹角的 减小值对时间的导数称为该平面上的角变形速率。
• 由图形几何关系，
1 2
1
2
（4）
1
tan 1
u zt z
z
u z
t
（5）
2
tan 2
有旋与无旋运动
• 流体质点的涡量定义为，
i jk
Ω = V = rotV = = 2ω （25）
x y z u vw
• 涡量表示流体质点绕自身轴旋转角速度的2 倍。并由涡量是否为零，定义无旋流动与 有旋运动。
有旋流动的一般性质
• 有旋流动又叫旋涡运动，它是流体运动的 一种重要类型。流动究竟是有旋还是无旋， 是根据流体微团本身是否旋转来决定的， 而不是根据流体微团的轨迹形状来决定的。
• 涡量场的散度为0，
Ω = V 0
涡线
• 涡线是这样一条曲线，曲 线上任意一点的切线方向 与在该点的流体的涡量方
向一致。因为 Ω 2ω ，
所以涡线也可看作流体微 团的瞬时转动轴线。涡线 是对同一时刻而言的，不 同时刻涡线可能不同。
涡线方程
• 根据涡线定义，
Ω dr = 0 （ dr 为涡线切线方向的向量）
实验流体力学作业-
实验流体力学作业1、试验流体力学作业答案第一题假设水翼所受到的升力用L表示，由题目可知，升力与速度v，水的气化压力流体密度p,水翼弦长攻角6T,水翼吃水深度/z,以上这些因素都会影响升力的大小，此题采纳7T定理的方法来解决，首先写出以下函数关系式：L=f、v，Pv，p，b，a，h、①取pv6三个量作为基本单位，其他物理量用无量纲数;T来表示：Lpv2b2兀，—b’71=a式屮为欧拉数，由于角度原来就是无量纲量，故不需要做无量纲化处理,由;T定理得到下式：习惯上用代替去，[emailprotected]，得到下式.?L1?12—2pV~lr由上式可知，欧拉数〃是与水翼所受到的升力相关的相像准则数。
第二题流体力2、学的试验设备品种繁多，主要试验设施冇两大类型：一类是进行水动力学试验讨论的水槽、水池、水洞等；另一类是进行空气动力学试验研宄的风洞。
第一类试验设施以循环水槽为例，其主要特点如下：便于观测，测量吋间不受限制，可以对物体四周的流态进行充分的观看。
不需要占用大量的建筑面积。
存在的主要问题是水流的匀称度和水面的平滑性不简单满足精确测定阻力的要求。
试验安装方便，试验模型小，有利于进行原理性、探究性的教学和科研试验,工作效率高。
第二类试验设施以风洞为例，其主耍的特点如下：风洞与循环水槽的结构组成类似，都是由收缩段、试验段等组成，唯一的不同是风洞的试验段可以是敞开式，空气不像循环水槽里的水一样3、要循环利用，而是直接排出到室外，其次流淌条件简单掌握。
＞第三题和似理论在现代科技中的最主要价值在于它指导模型试验上，确定“模型”与“原型”的相像程度、等级等。
因为一般原型都比较大，要讨论它的特性，必需将其缩小肯定的尺度，通过讨论模型的流淌特点，将流淌的参数通过相像理论还原到原型上，从而得到相应的结论。
相像理论是试验的理论，用以指异试验的根本布局问题，它为模拟试验提供指导，尺度的缩小或放太，参数的提高或降低，介质性能的转变等，H的在于以最低的成木和在最短的运转周期内摸清所讨论模型的内部规律性。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Programc-solving 1D burgers equation--------------------------------------------------c spatial discretization------------c time discretization---------------program mainc implicit noneparameter(nn=400)dimension u1(nn),u2(nn),x(nn),a(nn),b(nn),c(nn),d(nn),c2(nn)& ,d2(nn), e(nn),a2(nn),b2(nn),d3(nn),c3(nn),f(nn),f3(nn)Real u1 !velocity vectorReal u2 !velocity vectorReal dt !time stepReal dx !spatial stepreal x,l !position vector, total lengthReal tf !end timeReal t !current timeReal nu !viscosityReal r,k !r=nu*dt/dx**2, k=0.5*dt/dxReal a,b,c,d,c2,d2 !matrix coefficents for implicit schemeReal e,f,f3,c3,d3,a2,b2 !Real q,delta !coefficients for implicit schemeREAL lambda, mu !coefficients for u and u**3 termsREAL flxd,flxg !boundaries fluxes or velocity valueInteger n !number of spatial discretization pointsInteger boun,init !boundary type index, initialization type indexInteger scheme !scheme indexInteger i,j !loop indexesINTEGER it !current iterationInteger passauv !saving to file every passauv iterations!End declaration *************************************************************************! Reading Data from don file ************************************************************** ! write(*,*) 'good1'open(8,file='input.txt',status='old')read (8,*) lread(8,*) nread(8,*) dtread(8,*) tfread(8,*) passauvread(8,*) nuread(8,*) schemeread(8,*) initread(8,*) bounread (8,*) lambdaread (8,*) muread (8,*) flxgread (8,*) flxd! if (scheme==3) then! read(8,*) q! read(8,*) delta! end ifclose(8)! write(*,*) 'good1'! End reading data *********************************************************************** ! Initialization **************************************************************************dx=1./real(n)r=nu*dt/(dx*dx)k=0.5*dt/dxopen(7,file='res.txt')do i=1,n+1x(i)=(i-1)*dx*la(i)=0.b(i)=0.c(i)=0.d(i)=0.e(i)=0.f(i)=0.end dot=0.it=0! write(*,*) 'good1'! Setting the initialization type ******************************************************! init=0 sinusoidal! init=1 shock! init=2 1/(1+x)if(init .eq. 0) thendo i=1,n+1u1(i)=sin(2*3.1415927*(i-1)/n)end doendifif(init .eq. 1) thendo i=1,n+1IF (X(i)>l/2.) THENu1(i)=flxdelseu1(i)=flxgEND IFend doendifif(init .eq. 2) thendo i=1,n+1u1(i)=-2*nu/(x(i)+1)end doendif! write(*,*) 'good1'! End Initialization *******************************************************************!**************************************************************************************** ! TIME LOOP ******************************************************************************doif (t .gt. tf) exit !exit when end time is reachedt=t+dtit=it+1! Start scheme computation **************************************************************!scheme=1 explicit!scheme=2 Crank Nicholson!scheme=3 select one! write(*,*) 'good1'!**************************************************************************************! case (1) !scheme=1 explicitif (scheme .eq.1 ) thendo i=2,nu2(i)=r*u1(i+1)+(1-2*r)*u1(i)+r*u1(i-1)& -k*u1(i)*(u1(i+1)-u1(i-1))end dodo i=1,nu1(i)=u2(i)end doendif!***************************************************************************************! write(*,*) 'good1'! case(2) !scheme=2 Crank Nicholsonif (scheme .eq.2 ) then! Matrix coefficientsdo j=2,n+1a(j)=-0.25*dt*u1(j-1)/dx-0.5*rend dodo j=1,n+1b(j)=1+rend dodo j=1,nc(j)=0.25*dt*u1(j+1)/dx-0.5*rend dodo j=2,nd(j)=0.5*r*u1(j-1)+(1-r)*u1(j)+0.5*r*u1(j+1)&+lambda*u1(j)+mu*u1(j)**3end do! Boundariesif (boun==0) then !Dirichlet boundariesb(1)=1.c(1)=0.d(1)=flxgb(n+1)=1.a(n+1)=0.d(n+1)=flxdelse !Neumann boundariesd(1)=(1-r)*u1(1)+2*r*flxg+r*u1(2)+lambda*u1(1)+mu*u1(1)**3d(n+1)=(1-r)*u1(n+1)+r*u1(n)+2*r*flxd+lambda*u1(n+1)+mu*u1(n+1)**3endif! Thomas Algorithm (tridiagonal matrix equation)c2(1)=c(1)/b(1)d2(1)=d(1)/b(1)do i=2,n+1c2(i)=c(i)/(b(i)-a(i)*c2(i-1))d2(i)=(d(i)-a(i)*d2(i-1))/(b(i)-a(i)*c2(i-1))end dou1(n+1)=d2(n+1)do i=1,nu1(n+1-i)=d2(n+1-i)-u1(n-i+2)*c2(n+1-i)end doendif!scheme=3 **************************************************** write(*,*) 'good1'! case(3)if (scheme .eq.3 ) then! add one scheme! Matrix coefficientsdo j=2,n+1a(j)=-0.5*dt*u1(j-1)/dx-rend dodo j=1,n+1b(j)=1+2*rend dodo j=1,nc(j)=0.5*dt*u1(j+1)/dx-rend dodo j=2,nd(j)=u1(j)end do! Boundariesif (boun==0) then !Dirichlet boundariesb(1)=1.c(1)=0.d(1)=flxgb(n+1)=1.a(n+1)=0.d(n+1)=flxdelse !Neumann boundariesd(1)=(1-r)*u1(1)+2*r*flxg+r*u1(2)+lambda*u1(1)+mu*u1(1)**3d(n+1)=(1-r)*u1(n+1)+r*u1(n)+2*r*flxd+lambda*u1(n+1)+mu*u1(n+1)**3endif! Thomas Algorithm (tridiagonal matrix equation)c2(1)=c(1)/b(1)d2(1)=d(1)/b(1)do i=2,n+1c2(i)=c(i)/(b(i)-a(i)*c2(i-1))d2(i)=(d(i)-a(i)*d2(i-1))/(b(i)-a(i)*c2(i-1))end dou1(n+1)=d2(n+1)do i=1,nu1(n+1-i)=d2(n+1-i)-u1(n-i+2)*c2(n+1-i)end doendif! end select!End of scheme computation **************************************************************** !Saving Data to file res.txt intermediate step*******************************************! if (MOD(it,passauv)==0.) then! write(*,*) 't=',t,'s',it! write(7,*) 't=',t,'s',it! do i=1,n+1! write(7,*) x(i),u1(i)! end do! END ifwrite(*,*) 'i=',i, 'u(i)=', u1(25)!End Saving Data***************************************************************************end do! END TIME LOOP *************************************************************************** *******************************************************************************************!Saving last computation to file ********************************************************* write(7,*) 't=',t,'s'do i=1,n+1write(7,*) x(i),u1(i)END doclose (7)write(*,*) 'good1'!End saving *******************************************************************************! Tecplot saving-------------------------------------------------------------------------- open(45,file='plot.plt',status='unknown')WRITE(45,*) 'TITLE = " variables" 'WRITE(45,*) 'VARIABLES = x, u'WRITE(45,1021) n+1do i=1,n+1write(45,*)x(i),u1(i)enddoclose (45)1021 FORMAT(' ZONE, T="u", I=',I6)!------------------------------------------------------------------------------------------stopend!========================================================================================= tf分别为0.1s、0.15s、0.2s、0.3s时的速度曲线，其中：红色线是tf=0.1s时的速度曲线，橙色线是tf=0.15s时的速度曲线，绿色线是tf=0.2s时的速度曲线，蓝色线是tf=0.3s时的速度曲线。