《五类黎曼间断问题数值模拟》由会员分享，可在线阅读，更多相关《五类黎曼间断问题数值模拟（9页珍藏版）》请在金锄头文库上搜索。
1、五类黎曼间断问题数值模拟题目：1对于一维流 Euler 方程：UF0, （当地线性化后为UA U0）txtx试用通量的矩阵分裂格式的思想，构造此方程的数值解法。左右边界不考虑均使用 Neumann 边界条件：U0xx12编写 MatLab 程序计算本章第六节中所述的5 类黎曼间断问题中的任意三类问题。3模拟初场如下的非定常流， 并尝试分析常数 k 取何数值可使流场最大马赫数大于或小于 0.3：p(0, x)kmax( ax x0 ,0) 1u(0, x)0(0, x)Const可设 x1,1, a0.1, x0 0计算过程如下：1.通量的矩阵分裂格式：将方程组：UFt0x离散后得：U i n
2、1U i nt ( Fi 1/2 Fi 1/2 )x方程空间导数项的差分计算转化为网格单元界面处的通量计算。而双曲型问题的主要特点是存在特征线，并且扰动以特征速度沿特征线传播。不同的特征分量，传播方向不同。为体现迎风的思想，需要考虑这种方向性。 以一维情况为例，由于影响来自上游，对于向右传播的分量，应该采用左边的值计算，对于向左传播的分量，采用右边的值计算。按照这种思想，需要进行特征分裂。对于方程组UF0 ，当地线性化后为：UAU0 。可txtx将A写成A TT 1 的形式，其中diag( 1,2 ,n ) 。方程组左乘T 1，有：UUT 1T 1TT 10tx定义 U ch ，使得： T 1
3、U ch，则有：U ChU Ch0Utx将 分解成正部和负部，即，其中：1 () ，1 ()22F 的表对这两种方向不同的特征波，考虑其影响域，可以导出界面处通量达式。由于：FA UFRFLA(U RU L )F 也可以写成：F T T 1 U T (T 1 U) T UChT FCh对于特征形式的通量FChU Ch ，同样有：FCh =(FCh )R(FCh ) L( UCh )R (UCh ) L U Ch UCh UCh U ChFChFCh也即有： (FCh ) R(FCh )LFChFCh界面处的通量为：(FCh )界面 =(FCh )L +(FCh )R也可以直接由左侧或右侧得到：
4、(FCh ) 界面 =(FCh ) RFCh(FCh ) R+( FCh ) R( FCh ) R -(FCh )L (FCh ) 界面 =(FCh ) LFCh(FCh ) L +( FCh ) L( FCh ) R-(FCh )L 上式左乘 T ，得：F界面 =FRTFChF界面 =FLTFCh对从左右两侧得到的通量进行平均，则得：F1 (FRFL )1 T U Ch22“”表示代入 Roe 平均值进行计算。Roe平均值：uLDuR， HH LDH R， DRuD1D1LR ， a(1)(H12Lu)2前次作业已经导出矩阵 A的形式，并求出其特征值：1ua, 2ua, 3u ，且相应的矩阵
5、T 及向量 U ch 的形式为：111TMPuauauHuaHuau 2212 (pau)2aU ch T 1 U12 (pau)2ap2a因此，通量的矩阵分裂格式：U i n 1U i nt (Fi 1/2 Fi 1/2 )xF1 (FRFL )1 T UCh222.黎曼间断问题：( L ,uL , pL ),0x1一维流动，初始条件：(,u, p)(x,0)，( R ,uR, pR ),1x0即所有物理量在t=0 时刻的值在x=0 处有间断，左右两侧常数分布。这种初始条件的间断一般是不满足间断关系的，因而是不稳定的。在t0 后，这种初始间断立即分解为若干满足间断关系的间断及中心稀疏波或激波
6、。由于初始情况不同，一般该问题可以分解为下图五种类型中的一种。（虚线表示接触间断，单条线表示激波，一束线表示稀疏波）Riemann 问题解的类型的判别规则：因此对于该问题，用 roe 格式离散，初始条件按上面的方式给定，边界条件采用 Neumann 边界条件：U0x x 1即可编写程序进行计算，现计算出三类结果，罗列如下：（ 1）左、右都是激波L1, uL3, pL1.2 ，R1, uR0, pR0.01（ 2）左稀疏波，右激波L1,uL1, pL3 ，R1,uR0, pR0.01（ 3）中间间断面，两边稀疏波L2,uL0, pL1.5，R1, uR0, pR1上面三幅图中，速度、压力、密度同
7、时出现突变，对应于激波间断；速度、压力、密度同时出现连续变化，对应于稀疏波；而速度、压力不变，密度突变，对应接触间断。3.带压力脉冲的初场计算：按条件：p(0, x)k max( ax x0,0)1u(0, x)0(0, x)Const可设 x 1,1,a0.1, x00在程序中设定初场，即： pk max(0.1x ,0)1 ， u0 ，const ，在每个时间步计算各节点的马赫数分布，实际计算中发现，在时间推进的过程中，马赫数分布逐渐在两侧出现尖峰，之后尖峰增大，直到达到某个最大值，然后再减小，并且尖峰的位置向两侧传播。因此在每个时间步求出各节点马赫数的极大值，看其是否超过0.3。计算结果如下：CFL=0.5 ， N=100 时， k17.3 ，最大马赫数大于0.3。 k17.3 时，某一时刻的马赫数分布如下，最大马赫数0.30039经过计算发现，不同的CFL 数和节点个数N 的设置对计算结果有影响。节点个数增多，局部马赫数大于0.3 需要的 k 值减小； CFL 数减小，需要的 k 值则增大。具体如下：（ 1） CFL=0.5 ， N=200 时， k14.8 时即出现了马赫数大于0.3（ 2） CFL=0.3 ， N=100 时， k18.3 时出现了马赫数大于0.3
