计算流体力学大作业sod激波管

计算流体力学大作业目录1.一维无源稳态对流-传热稳态的计算流体力学分析问题介绍与数学方程方程离散化内节点法与外节点法计算结果分析结论2.有机玻璃圆管内的爆燃爆轰计算流体力学分析背景介绍模型简化差分格式间断问题精确解相同网格下不同格式的结果对比同一格式在不同网格下的结果对比问题介绍与数学模型一维无源稳态对流-扩散：在一个长L的计算区域内，流体的流速为u；物理左边界A处（x=0）通量Φ0=1，在物理左边界B处（x=L）通量ΦL=0；流体的密度为1.0kg/m3；扩散系数Γ=0.1kg/（m·s）。求解在不同流速下、不同网格划分下域内Φ的分布情况。数学模型：SΓVt无源：SΦ=0；稳态：非稳态项=0，简化为：dxdΓdxddxud方程离散化4dxdxdΓdxddxdxudewewWPePEeweWPPExΓxΓuu对简化后的控制方程在目标节点所在的控制容积内积分：扩散项采用中心差分，得：记对流强度F=ρu，扩散传导性D=Γ/Δx。Fe=(ρu)e，Fw=(ρu)w；De=Γe/ΔxPE，Dw=Γw/ΔxWP。针对本问题，采用均匀网格，Fe=Fw=ρu，De=Dw=Γ/Δx。代入上式，得到：WPPEweDDFF方法1：对流项中的界面值Φe、Φw采用中心差分：PWw21PEe21整理得离散节点的通用形式：EEWWPPaaaFDaW21FDaE21EWPaaa方法2：对流项中的界面值Φe、Φw采用一阶上风：Wwuw˃0时，ue˃0时，PeEEWWPPaaaFDFDaW0max，DFDaE0max，EWPaaa整理得离散节点的通用形式：方法3：对流项中的界面值Φe、Φw采用混合格式，根据贝克力数选择界面值的取值：Pe˃2：WwWPw21PePe2：EPe21EEWWPPaaa整理得离散节点的通用形式：021max，，FDFaW021max，，FDFaEEWPaaa内节点法与外节点法1.内节点法即在构造控制容积时先确定出控制容积界面的位置，然后将节点取在每一个控制容积的几何中心。•对节点1的控制容积进行积分：•对节点n的控制容积进行积分：不同的离散格式将对Φe1、Φwn进行不同的处理，得到针对该两个节点的各自的离散方程。2.外节点法即在构造控制容积时先确定出节点的位置，然后使控制容积界面位于两个节点的正中间。不需要单独对节点1与节点n进行离散，离散格式获得的离散方程适用于计算域内所有节点。dxdxdΓdxddxdxudeAeA11AAeAAeexΓxΓuu1121112/dxdxdΓdxddxdxudBwnBwn12/nnwnMBBwnwnBBxΓxΓuu计算结果分析运用Fortran编程，将三种格式与精确解糅合在一起。输入计算区间的左右端点A、B的通量值，运用Jacobi迭代求解每种格式对应的线性方程组。在计算之前先需要给计算域赋一个初场，通过比较新一轮迭代结果与上一轮迭代的结果差值不大于0.00001作为迭代收敛的判据。（1）网格数n=5；u=0.1m/s（Pe=0.2，F=0.1，D=0.5）1.内节点法：节点位置解析解数值解与误差率中心差分误差率一阶上风误差率混合误差率1x=0.10.93880.94210.00350.9348-0.00430.94210.00352x=0.30.79640.80100.00580.7914-0.00630.80060.00533x=0.50.62250.62760.00820.6193-0.00510.62760.00824x=0.70.41000.41620.01510.41290.00710.41630.01545x=0.90.15050.15790.04920.16520.09770.15790.0492表1u=0.1m/s，Δx=0.2m结果对比表三种格式与解析解对比图局部放大图（2）网格数n=5；u=2.0m/s（Pe=4.0，F=2.0，D=0.5）不同格式与解析解对比图是否满足守恒性、输运性、有界性？有界性判据：①斯卡巴勒准则；≤1（至少一个离散方程满足小于1）②正系数；③源项负斜率。中心差分：节点n˃˃2D一阶上风：˃混合格式：靠近边界B的结果较解析解偏大，这是由于混合格式在处理边界节点是过高估计了上游的结果。PnbaanbaPaFnbaPado10,i=1,M!中心差分fai(i)=0.110continuedo20,j=1,1000do11,i=1,Mxfai(i)=fai(i)11continuefai(1)=((2*D+F)*fai0+(D-0.5*F)*fai(2))/(3*D+0.5*F)do30,i=2,M-1fai(i)=((D+0.5*F)*fai(i-1)+(D-0.5*F)*fai(i+1))/(2*D)30continuefai(M)=((D+0.5*F)*fai(M-1)+(2*D-F)*fai1)/(3*D-0.5*F)if(abs(fai(M)-xfai(M)).lt.0.00001)thengoto12endif20continue12write(2,*)aa+deltx/2,fai(1)do60,i=2,M-1y=aa+deltx/2+(bb-aa)*real(i-1)/real(M)write(2,*)y,fai(i)60continuewrite(2,*)bb-deltx/2,fai(M)close(unit=2)do13,i=1,M!一阶上风upfai(i)=0.113continuedo70,j=1,1000do14,i=1,Mxfai(i)=upfai(i)14continueupfai(1)=((2*D+F)*fai0+D*upfai(2))/(3*D+F)do80,i=2,M-1upfai(i)=((D+F)*upfai(i-1)+D*upfai(i+1))/(2*D+F)80continueupfai(M)=((D+F)*upfai(M-1)+(2*D-F)*fai1)/(3*D)if(abs(upfai(M)-xfai(M)).lt.0.00001)thengoto15endif70continue15write(4,*)aa+deltx/2,upfai(1)do300,i=2,M-1y=aa+deltx/2+(bb-aa)*real(i-1)/real(M)write(4,*)y,upfai(i)300continuewrite(4,*)bb-deltx/2,upfai(M)close(unit=4)（3）网格数n=20；u=2.0m/s（Pe=1.0，F=2.0，D=2.0）不同格式与解析解对比图离散方程满足守恒性、有界性和输运性三个物理特征。数据显示，混合格式的结果与中心差分的结果一致，相较于一阶上风格式误差较小，因为该情况下扩散占的比重较对流大，一阶上风会过高估计上游信息对下游的影响。2.外节点法：（1）网格数n=4；u=0.1m/s（Pe=0.2，F=0.1，D=0.5）为了获得与内节点法相同的网格间隔，外节点法取n=4，此时两种网格划分方式对应的F、D相同，可以认为是对同一物理问题的求解。不同格式与解析解对比图节点位置解析解数值解与误差率中心差分误差率一阶上风误差率混合误差率1x=0.20.87110.87110.00000.8612-0.01140.8669-0.00482x=0.40.71380.71380.00010.6980-0.02210.7078-0.00843x=0.60.52150.52170.00030.5053-0.03120.5163-0.01014x=0.80.28680.28690.00060.2756-0.03890.2840-0.0098节点位置解析解数值解与误差率中心差分误差率一阶上风误差率混合误差率1x=0.10.93880.94210.00350.9348-0.00430.94210.00352x=0.30.79640.80100.00580.7914-0.00630.80060.00533x=0.50.62250.62760.00820.6193-0.00510.62760.00824x=0.70.41000.41620.01510.41290.00710.41630.01545x=0.90.15050.15790.04920.16520.09770.15790.0492表3u=0.1m/s，Δx=0.2m外节点法表2u=0.1m/s，Δx=0.2m内节点法外节点法的中心差分格式结果误差要小于内节点法的结果，即外节点法在处理一维无源弱对流-扩散问题上要优于内节点法。（2）网格数n=4；u=2.0m/s（Pe=4.0，F=2.0，D=0.5）不同格式与解析解对比图外节点法对应的一阶上风格式在u=2.0m/s，Δx=0.2m并没有出现像内节点法那样数值超过0~1范围的情况。分析原因知，因为外节点法在处理节点n时未单独进行离散，没有出现˃的异常，并且自动满足，符合有界性准则。虽然一阶上风格式得到了收敛的数值解，但在节点n处计算结果误差较大。混合格式数值解误差率不大，但略微高估了上游的影响。nbaPanbaPa节点位置解析解数值解与误差率中心差分误差率一阶上风误差率混合误差率1x=0.21.0000-1.4E+292**0.9987-0.00131.00000.00002x=0.41.0000-5.4E+292**0.9923-0.00771.00000.00003x=0.60.9997-1.3E+293**0.9603-0.03941.00000.00034x=0.80.9817-2.0E+293**0.8002-0.18481.00000.0187（3）外节点法改进，网格数n=4；u=2.0m/s（Pe=4.0，F=2.0，D=0.5）外节点法在处理该一维无源对流-扩散问题上，避免了内节点法带来的其中一个离散方程˃的不恰当方程，解决了越界的异常。但是一阶上风格式在节点n处过低估计了上游的影响，而混合格式稍微高估了上游的影响，故对外节点法提出改进措施。式（1）、式（2）分别为一阶上风格式与混合格式在节点n处的离散方程：upfai(M)=((D+F)*upfai(M-1)+D*fai1)/(2*D+F)（1）mixfai(M)=(aw*mixfai(M-1)+ae*fai1)/ap（2）nbaPa相应增大和减小上游节点对n节点的影响程度，并且要保证，故将式（1）、式（2）改造为以下式（3）、式（4）：nbaPaupfai(M)=(1.15*(D+F)*upfai(M-1)+D*fai1)/(2*D+F)（3）mixfai(M)=(aw*mixfai(M-1)+ae*fai1)/(1.05*ap)（4）节点位置解析解数值解与误差率中心差分误差率一阶上风误差率混合误差率1x=0.21.0000-1.4E+292**0.9997-0.00031.00000.00002x=0.41.0000-5.4E+292**0.9980-0.00201.00000.00003x=0.60.9997-1.3E+293**0.9897-0.01001.00000.00034x=0.80.9817-2.0E+293**0.9485-0.03380.9524-0.0299不同格式与解析解对比图表4u=2.0m/s，Δx=0.2m改进外节点法一阶上风在节点n