双曲型守恒律组的高阶Godunov格式

双曲型守恒律组的高阶Godunov格式丁岩1.Godunov格式简述模型方程为下述双曲型守恒律组:dUdF(U)八——+———=0

dt dxU(x,0)二U

o其中U=(p,pu,pv,E)t，F(U)=(pu,pu2+p,u(E+p))t。假设所采用的是均匀网格，Ax是网格宽度。给定tn时刻的单元均值分布｛Un｝，构造分片常数分布函数:j问题:Un（x）=Un，当x<x<x

j j j—1/2 j+1/2问题:在每个边界x处，（近似）求解如下初始条件的j+1/2U=Un，U=UnLjR j+1对充分小的时间步长At（满足CFL<1条件），记tn+1=tn+At，可得到具有相似特性的近似解x—XU(x,t) U( j^2)，x<x<x，tn<t<tn+iTOC\o"1-5"\h\zj jt—tn j—1/2 j+1/2从而有整个计算区域内的近似解：U(x,t)=U(x,t)，x<x<x，tn<t<tn+Atj j—1/2 j+1/2将它代入 式，并在x<x<x， tn<t<tn+At上积分有：j—1/2 j+1/2Jxj+1/2U(x,"+i)dx=Jxj+1/2U(x,"+i)dx=xjXj-1/2xjxj-1/2tnF[U(xjj-1/2/)]dt-J=F[U(xtn jj+1/2注意到〃（X）是是 问题沿射线x=x，tn<t<tn+At的解，因而TOC\o"1-5"\h\zjj+1/2 j+1/2是一常数，于是得到如下格式：— — At- 。Un+— — At- 。Un+1=Un- (F—F)，jj Ax j+1/2 j—1/2中给出了利用问题的解计算通量时需要考虑的十种j+1/2 j+1/2中给出了利用问题的解计算通量时需要考虑的十种情形之一。2.二阶Godunov格式(MUSCL)通过以下步骤计算数值通量F ，以实现j+1/2格式的二阶推广。II=[x,x]中，由j j—1/2 j+1/2U(x)=Un+ j-A，xe[x,x]j j Ax j j—1/2 j+1/2不同于 格式采用分片常数插值，这里在每个tn时刻的均值｛Un｝构造分片线性函数，即jx一xA是适当选择的U(x)在I上的斜率j j j端点x,x 处的值为：j-1/2j+1/2八八1，UL=Un——A，jj2j

矢量(差分形式)，U(X)在I的两个j j称之为是具有三个分量的矢量，因而中实际有六个称之为是具有三个分量的矢量，因而中实际有六个注意到对于在每个l中，将由给出的 按如下格j式在时间上推进1At：2Ul-UL+1竺[F(UL)—F(UR)]TOC\o"1-5"\h\zj j 2Ax j jUr-Ur+1A[F(Ul)—F(Ur)]

j j 2Ax j j在每个 边界点处现在存在有两个通量F(Ul)和F(Ur)，而且他们一般来说是j j问题：格式中完全问题：格式中完全为计算间的通量F，可以通过求解如下初始条件的j+1/2U-Ur，U-UlTOC\o"1-5"\h\zLjR j+1来得到相似性解u (x-xj+1/2)，则数值通量F可以按照j+1/2 t—tn j+1/2相同的方式来得到，即人F-F(U (0))j+1/2 j+1/2— —At」 。Un+1-Un— (F—F)j j Ax j+1/2 j—1/2注：关于斜率A的选择j

中，A可以简单的取作

j1+-1+-(1-w)A

j-1/2 2 j+1/2A=-(1+w)A其中A=Un—Unj-1/2 j j-1A其中A=Un—Unj-1/2 j j-1A三j+1/2Un-Un且WG[-1,1]j+13=0时，上述格式正是二阶的方法。j根据定理，在大梯度变化区域，上述格式将产生数值振荡。所以中A的取法不能保证格式的j性质。可以采用如下的A替代上面的Aj其中的0其中的0(r)是一斜率限制子j足以下条件的g(r)是 的：，采用如此A的格式对于满

j0<己(r)<min{1(r),己(r)}，r>0LR其中(r)(r)=自(r)=RA

r=-j-1-3+(1+3)r21-3+(1+3)rj+1/2基于此，可以构造以下满足上述条件是斜率限制子:TOC\o"1-5"\h\zQ r<0g(r)=J2r， 0<r<1/2sb1, 1/2<r<1min{r,g(r),2}, r>1L R

0,《 2rmin{7—1 1+rm(r)m(r)=<r,

mbmin{1,m(r)},

R各种限制子与 区域间的关系（3.三阶Godunov格式（PPM）方法（）方法（）是所研1的一种三阶本插值函数，来代替所采用的常数函数和所采用的线性函数。种三阶本插值函数，来代替所采用的常数函数和所采用的线性函数。格式。它的主要特点是采用二次多项式函数作为网格内部的基先讨论单个线性方程情况:du du „—+a——0，a—const.dt dx

网格j的空间步长Ax=x-x 可以是不等距的。j j+1/2j-1/2已知tn时刻的离散分布｛Un｝，如解在网格j中的平均值。假设根据｛Un｝构造了分j j段抛物线函数：Un(x)=U +0[(u —u)+u(1—0)]j L,j R,j L,j6,j=u—^[(u—u)—u (1—^)]R,j R,j L,j6,jTOC\o"1-5"\h\zxx—x x—xq= j-i/2— n=——j+i/2 x—x x—xj+1/2 j-1/2 j+1/2 j-1/2利用q+n=1，匕(1—q)=n(1—n)容易验证 式给出的二次多项式的两种表达式是等价的。要求所构造的函数un(x)满足:1邛…

un= un(x)dxjAxjjxj—1/2相当于成立以下关系式:i/ 、iun=u

j+ (u —uun=u

jL，j 2凡jL，j6 6,j由此确定ci/u=6[un-—(u+u由此确定ci/u=6[un-—(u+u)]6,jj2R，jL,jun(X)中的另外两个参数u,u与｛un｝的关系按以下方式确定:j L,jR,j j首先根据｛un｝计算u在X 处的近似值u ，j j+i/2 j+i/2X处的值可以由下式得到：j+i/2因为u的不定积分4X)在网格边界A(x )=A =ZunAxj+i/2 j+i/2 kkk<j为计算u，j+i/2构造通过五个点(A,x)，k=0,±i,±2的四次插值多项式，并j+k+i/2 j+k+i/2对其微分得到u=竺Ij+i/2dxX山2求uj+i/2的公式可以用un,un,un和AX,M,M,Mj-i j j+ij-i j j+i j+2 X2 Ax X2 Axk=-i j+k2AxAx

j+i__j-

Ax +Axj+ijAx +Ax——j^i j-2Ax+Axj j+iAx+Ax

―j+i :Ax+2Ax

j .j2(un-un): j+i jj+1」表示如下:Ax

u=un+ j (un —un)+j+i/2 j Ax+Ax j+i j乙j j+iAx+Ax Ax+AxTOC\o"1-5"\h\z-Ax——j- j-8u+Ax—j+ 48uj2Ax+Ax j+i j+iAx+2Ax jj j+i j j+i这里的8u是所构造的二次插值函数在网格j中的近似斜率，由下式给出:jAx8u= j j Ax+AxAx8u= j j Ax+Ax+Axj-i j j+i j- j-(un—un)+ j j+i(un—un)Ax+Ax j+i j Ax +Ax jj-ij j+i j-i j在实际的计算中，还须对8u作单调性处理，即采用以下8u代替8u：sgn(8u)min(8sgn(8u)min(8u,2j0,其它un-unj j-i,2un—un),

j+i j如果3〃-un)(un-un)>0

j+i jj j-i这种修正不仅可以得到更加锐利的捕捉到解的间断，而且能够保证u 落在unj+i/2 j和un的范围之内。j+i上述计算u 的公式对非均匀的网格具有三阶精度，如果网格是均匀的，则:j+1/2ujuj+1/2由和得到u后，j+1/27/ 、1/ 、二—(un+un)—-(un-un)12jj+112j+2j-1即令：u—uLu—uL,j+1R,j在大多数情况下，由上式得到的u,uL,jR,j=u， Vjj+1/2就是二次多项式un(x)表达式 中的两端值。但在以下两种情况，还需要对u,uL,jR,jj重新赋值进行调整:当(u-un)(un-u)<0时，则在此网格中有局部极值。这时重新对u,u赋值:R,jjjL,j L,jR,ju—u—unR，jL,jj使得在该网格中插值多项式un(X)三un。j jj作法是：将u,uR,jL,j当unj作法是：将u,uR,jL,j值；当U落在靠近U的一段中时，修改值；当U落在靠近U的一段中时，修改U值。R,j用数学表达式来描述这种作R,j法时，有:、/ 1/ 、、、/ 1/ 、、1/)(un——(u +u))>—(uj2 R,j L,j 6-u)2时，R,jL,ju保持不变，重新赋值R，ju=3u=3un-2u 。L,j j R,j当(u -uR,j L,j、/ 1/ 、、)(un—(u +u))<-j2R，j一1(u-u)2时，u保持不变，重新赋值6R，jL，j L,ju =3u =3un-2u 。R,j j L,j一旦得到了u和u的值之后R,jL,jun(X)的平均：就可以写出计算un+1的格式，首先定义插值函数jfuj+1/2,L,、1h,、，(y)= j+1/2un(x)dX其中的y表示一正数，fuj+其中的y表示一正数，fuj+1/2,R容易验证:j+1/2X.,j+1/2j,+1/2+yUn(X)dxfu (y)=u ——y—(u —u—(1-—y^—)u),TOC\o"1-5"\h\z1+1/2,^ r,j 2Ax R，j l,j' 3Ax 6,/j jfu (y)=u ————(u —u+(1— ——)u )j+1/2,R L,j+12Ax R，j+1 l,j+1 3Ax 6,j+1j+1 j+1则un+1可以按以下格式得到：jun+1=un+a (u —u)jj Ax j-1/2 j+1/2j其中：fu (aAt), 若a>0u=1j+1/2,Lj+1/2 fu (—aAt), 若a<0lj+1/2,R下面将讨论，如果将单个方程的 方法推广应用到 坐标的流体力学方程组:StSt 3(rau)St 3m0,TOC\o"1-5"\h\zSu Sp 八—ra=0,

St SmSES(rapu)—— =0St Sm其中的t=1/p为比容，p为密度，u为速度，p为压力，E=e+2u2为总能，e为内能，状态方程为p=(y-1)pe,y>1,r为空间坐标，m为 质量坐标，它们之间有以下关系式相联系:m=Jrpradrt为时间，对t求导为导数a=0,1,2分别对应于平面对称，柱对称，球对称的一维情形。Am―m—m，j j+1/2 j—1/2令Am为Am―m—m，j j+1/2 j—1/2j j—1/2 j+1/2这些量不随时间变化。现设给定了tn时刻网格均值分布｛Tn,un,En｝以及网格端点jjj

的空间坐标｛r的空间坐标｛rn ｝，j+1/2(m,tn)dmj-1/2要求用方法计算tn+1=tn+At时刻的网格平均值｛Tn+1,Un+1,En+1｝及端点位置g+1｝要求用方程组的 算法的计算步骤如下:j+1/2方程组的 算法的计算步骤如下:由网格平均值Tn,Un,En用以下公式计算压力的网格平均值小：/ 、 L1/ 、pn=p(Tn,en)，en=En一一Qn)2

jjjjj2j其计算精度为二阶精度。应用线性方程情况所描述的构造分段抛物函数Un(X)的方法，只需要以质量坐标m代替前述的x，分别对所给的网格均值分布｛Tn｝,｛un｝,｛坐标m代替前述的x，jjj构造相应的分段抛物线函数Tn(X),Un(X)pX，)确定其中有关的参数{T ,T ,T }u{ u, u, }p{p,p,。}L,jRj6,jLjRj6,jLjRj6j,p±：这两个区j+1/2按下式定义T,u,p,p±：这两个区j+1/2j+1/2 j+1/2m之间的区域。j+1/2域分别是经过点(m ,tn+1)的土特征线与｛tm之间的区域。j+1/2j+1/2e+ =fe (AtCnAn)j+1/2 j+1/2,L jje- =fe (AtCnAn)j+1/2 j+1/2,R j+1j+1这里e=t,u,p，c2=ypp，TOC\o"1-5"\h\z(rn )a+1-(rn )a+1An= j+1/2 M/2 j (a+1)(rn-rn)j+1/2 j-1/2分别以e=e+和e=e- (e=T,u,p)为左右状态求解黎曼问题，得到j+1/2,L j+1/2 j+1/2,R j+1/2速度和压力在网格边界j+速度和压力在网格边界j+2处的时间平均值Uj+1/2,P+1/2。所以在这儿作了与一阶方法类似的工作，只是这里的左右