解欧拉方程为什么会得到激波？

robinlxw · 发表于 2006-9-16 08:58:32

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？注册

x

去逛了一下流体网站，有一个topic是正本清源我觉得挺好的，搞清概念是很重要的。这里把我以前考虑过的问题也拿来大家讨论一下，我当时是觉得很有意思的：
问题1：欧拉方程组的能量方程实际上就是一个等熵方程，所以，解欧拉方程应当得到一个等熵的流场，但是，事实上，我解过欧拉方程，用的是有限差方法，TVD格式，可以得到激波，而激波前后是不等熵的，也就是说，我解得的流场不是一个等熵流场，这不是互相矛盾吗？
为什么？？？？，
忙去了，以后再加一些问题进来，呵呵

流星客 · 发表于 2006-9-16 19:20:29

你看看等熵流动的基本假设,定常是一个必要条件.求解EULER方程是引入了非定常计算.这只是问题个一个方面.

robinlxw · 发表于 2006-9-16 22:19:07

欧拉方程是等熵方程并不需要定常假设，实际上等熵方程（能量方程）就是熵的质点导数为0，而质点导数可分解为对时间的偏导+...，所以说并没有引入定常假设

流星客 · 发表于 2006-9-17 00:49:47

建议看看空气动力学学基础吧

robinlxw · 发表于 2006-9-17 08:31:00

我错了？不用看空气动力学的书，流体力学基础，清华大学，潘文全，里面写得很清楚

鲁秦 · 发表于 2006-9-17 22:24:23

数值黏性问题

robinlxw · 发表于 2006-9-18 08:47:42

总版主你真强啊，十分佩服，和我的想法是一样的，我只是想印证一下我的想法，这边都没有你总版主这么强的高人，谢谢总版主。可不可以加你的QQ啊？或者电子邮件？我的电子邮件Robinlxw@163.com，[br][br][以下内容由 robinlxw 在 2006年09月18日 08:53am 时添加] [br]
纠正一下，我说这边没总版主这么强的高人，是指我们这边，不是是论坛上面

鲁秦 · 发表于 2006-9-18 13:01:49

我不怎么上QQ
给我发邮件吧,我抽时间看，工作有点忙
chuanchao_liu@163.com

spacewarship · 发表于 2006-9-20 21:53:11

好象前提不对吧
我们通常求解的欧拉方程只是忽略了粘性力和热传导的作用,并没有引入等熵假设,
欧拉方程组中的能量方程并不是等熵方程.
而欧拉方程在超音速流时无论是定常还是非定常,都是双曲型的,仍然满足激波产生的条件,所以欧拉方程的理论解也可以有激波的产生,只不过这时的解应是理论上没有厚度的无粘激波.
但是实际的数值解由于引入了差分格式固有的数值粘性,才使得数值解出的激波有被抹平或者振荡的特性,又与实际流动中激波现象一致

robinlxw · 发表于 2006-9-21 09:27:54

回复9楼：
我们在数值解的时候是没有引入等熵假设，因为根本就不需要，欧拉方程组中的能量方程本来就是等熵方程，参考流体力学基础，清华大学，潘文全，144页
恰恰相反，数值求解时，对于有激波并且用激波捕捉法的格式，还一般都要求格式是变熵的，这一点也是构造构式的要点之一，参考张院士的书差分方法的原理及应用
欧拉方程在超音速时是双曲型的说法是错误的，这是众所周知的。对于钝头体，正激波后的亚音速区是扫椭圆型的，而后体的超音速区才是双同曲型的，正因这个，所以才要用非定常的方法解钝头体的欧拉方程，参考候天相老师的《非定常方法求解...》。
您提到的最后一点倒是事实，不过我已经不搞这个很多年了，也不能再对它进行深入的研究了，不过并不能充分解释我的问题，我还是认为总版主和我自己的观点是对的。
大家多多讨论，或许我是错的

yangeryanzi · 发表于 2006-11-9 23:45:40

真强呀！都是牛人！

lwd1981 · 发表于 2007-3-3 18:50:15

我觉得这里面实际上有一个至关重要的问题,被上面的讨论所忽略了.那就是等熵与绝热的区别.等熵是绝热+可逆.等熵关系和激波前后的绝热关系在积分形式下在数学形式上是一样的,在超音速情况下,绝热关系是不能够全程微分的,但是等熵关系可以,这是二者的区别.另外我们谈欧拉方程计算激波的时候,提法实际上并不准确.特别是对于熟悉CFD的人士而言.我们知道再进行计算时要采用相应的方程形式,如微分形式、积分形式;守恒形式、非守恒形式等.我们用于计算激波的形式一般是积分而且是守恒形式.此时的欧拉方程可以说是绝热方程,也可以说是等熵方程.但本质上讲应该是绝热方程才更为准确.因为前面讲了等熵关系和激波前后的绝热关系在积分形式下在数学形式上是一样的.从这个角度上讲我们不能够得到欧拉方程算不出激波的结论,因为激波的存在是个局部间断,而此时的欧拉方程是个积分(整体,准确的讲是激波前后整体所满足的)方程.这是从物理理论上的分析.另外,在从CFD计算的角度来分析一下,CFD的计算格式(主要讲激波捕捉格式)的构造一般都满足了熵增差分关系,来避免激波前后的震荡,也就是人为加入了数值粘性.同时,还有各种格式计算过程中所存在的色散和耗散,也是一种数值耗散,他们提高了激波的分辨率,这也使得CFD中欧拉方程可以算出激波.因此,无论从物理还是从数值方面我们都不能够得到欧拉方程算不出激波的结论!
另外,顺便说一句.在超音速时,无论是定常还是非定常情况,欧拉方程都是双曲型的!在定常跨音速时是混合型的.在定常亚音速时是椭圆型的!在非定常情况,欧拉方程都是双曲型的!![br][br][以下内容由 lwd1981 在 2007年03月20日 03:24pm 时添加] [br]

nozzle · 发表于 2007-4-2 14:05:01

我觉得楼主有两个方面需要再仔细考虑一下：
1、欧拉方程忽略了黏性，但黏性与激波的产生无关，故得不出欧拉方程算不出激波的结论；
2、欧拉方程的微分形式是等墒的，这是做了流场所有参数连续的假设，而积分形式的方程无这样的假设，仅仅绝热；计算存在激波这样的间断面、存在墒增的过程只能使用积分形式；在CFD中为捕获激波从而引入数值黏性，这正是描述此过程的欧拉方程在CFD中的数值描述！

panjiaying · 发表于 2011-4-27 09:14:53

   !二维欧拉方程
program NND2d
implicit real(a-h,o-z)
parameter(nx=400,ny=400,ighost=2)
real  u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real minmod

gama=1.4
dx=1.0/nx
dy=1.0/ny

call init(u,nx,ny,ighost,gama)

dt=CFL(u,dx,dy,nx,ny,ighost,gama)
dt=dx*dy
   supt=0.08
   t=0
   do while(t.lt.supt)
     if(t+dt.gt.supt) dt=supt-t
t=t+dt
   write(*,*) 't=',t,'dt=',dt

   call solveall(u,nx,ny,dt,dx,dy,ighost,gama)

call bound(u,nx,ny,ighost)

enddo

   call output(u,nx,ny,dx,dy,ighost,gama)

   end

!子程序

subroutine init(u,nx,ny,ighost,gama)!初始化
implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)

do i=-ighost,nx/2
do j=-ighost,ny/2
u(i,j,0)=0.5
u(i,j,1)=0.5*3.0
u(i,j,2)=0.5*1.0
u(i,j,3)=1.0/(gama-1)+0.5*0.5*10
enddo
enddo

do i=nx/2+1,nx+ighost
do j=-ighost,ny/2
u(i,j,0)=1.5
u(i,j,1)=1.5*3.0
u(i,j,2)=-1.5*1.0
u(i,j,3)=1.0/(gama-1)+0.5*1.5*10
enddo
enddo

do i=nx/2+1,nx+ighost
do j=ny/2+1,ny+ighost
u(i,j,0)=0.5
u(i,j,1)=-0.5*3.0
u(i,j,2)=-0.5*1.0
u(i,j,3)=1.0/(gama-1)+0.5*0.5*10
enddo
enddo

do i=-ighost,nx/2
do j=ny/2+1,ny+ighost
u(i,j,0)=1.5
u(i,j,1)=-1.5*3.0
u(i,j,2)=1.5*1.0
u(i,j,3)=1.0/(gama-1)+0.5*1.5*10
enddo
enddo

end

subroutine solveall(u,nx,ny,dt,dx,dy,ighost,gama)
implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real gt(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ft(-ighost:nx+ighost,-ighost:ny+ighost,0:3)

   call Solvex(u,ft,nx,ny,ighost,gama)
   call Solvey(u,gt,nx,ny,ighost,gama)

r1=dt/dx
r2=dt/dy

do i=-ighost+2,nx+ighost-2
do j=-ighost+2,ny+ighost-2
do k=0,3
u(i,j,k)=u(i,j,k)-r1*(ft(i,j,k)-ft(i-1,j,k))
   & -r2*(gt(i,j,k)-gt(i,j-1,k))
enddo
enddo
enddo

end

subroutine Solvex(u,ft,nx,ny,ighost,gama)
implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
   real fl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real fr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real dfl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real dfr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ffl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ffr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ft(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real minmod

eps=1e-6

   do j=-ighost,ny+ighost
do i=-ighost,nx+ighost
p=(gama-1.0)*(u(i,j,3)-0.5*u(i,j,1)**2/u(i,j,0)
   & -0.5*u(i,j,2)**2/u(i,j,0))
      uu=u(i,j,1)/u(i,j,0)
c=sqrt(gama*p/u(i,j,0))
vv=u(i,j,2)/u(i,j,0)
H=c*c/(gama-1)+0.5*(uu**2+vv**2)

a1=uu
a2=uu-c
a3=uu+c

a1l=0.5*(a1+sqrt(a1*a1+eps*eps))
a2l=0.5*(a2+sqrt(a2*a2+eps*eps))
a3l=0.5*(a3+sqrt(a3*a3+eps*eps))

a1r=0.5*(a1-sqrt(a1*a1+eps*eps))
a2r=0.5*(a2-sqrt(a2*a2+eps*eps))
a3r=0.5*(a3-sqrt(a3*a3+eps*eps))

fr(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1r+a2r+a3r)
fr(i,j,1)=u(i,j,0)/2/gama*(2*uu*(gama-1)*a1r+(uu-c)*a2r
   &       +(uu+c)*a3r)
fr(i,j,2)=u(i,j,0)*vv/2/gama*(2*(gama-1)*a1r
   & +a2r+a3r)
fr(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1r
   & +(H-c*uu)*a2r+(H+uu*c)*a3r)

fl(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1l+a2l+a3l)
fl(i,j,1)=u(i,j,0)/2/gama*(2*uu*(gama-1)*a1l+(uu-c)*a2l
   &       +(uu+c)*a3l)
fl(i,j,2)=u(i,j,0)*vv/2/gama*(2*(gama-1)*a1l
   & +a2l+a3l)
fl(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1l
   & +(H-c*uu)*a2l+(H+uu*c)*a3l)

enddo

do i=-ighost+1,nx+ighost
do k=0,3
dfl(i,j,k)=fl(i,j,k)-fl(i-1,j,k)
dfr(i,j,k)=fr(i,j,k)-fr(i-1,j,k)
enddo
enddo

do i=-ighost+1,nx+ighost-1
   do k=0,3
y=minmod(dfl(i,j,k),dfl(i+1,j,k))
ffl(i,j,k)=fl(i,j,k)+0.5*y
enddo
enddo

do i=-ighost+1,nx+ighost-2
do k=0,3
y=minmod(dfr(i+1,j,k),dfr(i+2,j,k))
ffr(i,j,k)=fr(i+1,j,k)-0.5*y
enddo
enddo

do i=-ighost+1,nx+ighost-2
do k=0,3
ft(i,j,k)=ffl(i,j,k)+ffr(i,j,k)
enddo
enddo
enddo
end

subroutine Solvey(u,gt,nx,ny,ighost,gama)
implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
   real gl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real gr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real dgl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real dgr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ggl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real ggr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real gt(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
real minmod

eps=1e-6
   do i=-ighost,nx+ighost
do j=-ighost,ny+ighost
p=(gama-1.0)*(u(i,j,3)-0.5*u(i,j,1)**2/u(i,j,0)
   & -0.5*u(i,j,2)**2/u(i,j,0))
      uu=u(i,j,1)/u(i,j,0)
c=sqrt(gama*p/u(i,j,0))
vv=u(i,j,2)/u(i,j,0)
H=c*c/(gama-1)+0.5*(uu**2+vv**2)

a1=vv
a2=vv-c
a3=vv+c

a1l=0.5*(a1+sqrt(a1*a1+eps*eps))
a2l=0.5*(a2+sqrt(a2*a2+eps*eps))
a3l=0.5*(a3+sqrt(a3*a3+eps*eps))

a1r=0.5*(a1-sqrt(a1*a1+eps*eps))
a2r=0.5*(a2-sqrt(a2*a2+eps*eps))
a3r=0.5*(a3-sqrt(a3*a3+eps*eps))

gl(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1l+a2l+a3l)
gl(i,j,1)=u(i,j,0)/2/gama*(2*vv*(gama-1)*a1l
   & +(vv-c)*a2l+(vv+c)*a3l)
gl(i,j,2)=u(i,j,0)*uu/2/gama*(2*(gama-1)*a1l
   & +a2l+a3l)
gl(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1l
   & +(H-c*vv)*a2l+(H+vv*c)*a3l)

gr(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1r+a2r+a3r)
gr(i,j,2)=u(i,j,0)/2/gama*(2*vv*(gama-1)*a1r
   & +(vv-c)*a2r+(vv+c)*a3r)
gr(i,j,1)=u(i,j,0)*uu/2/gama*(2*(gama-1)*a1r
   & +a2r+a3r)
gr(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1r
   & +(H-c*vv)*a2r+(H+vv*c)*a3r)

enddo



do j=-ighost+1,ny+ighost
do k=0,3
dgl(i,j,k)=gl(i,j,k)-gl(i,j-1,k)
dgr(i,j,k)=gr(i,j,k)-gr(i,j-1,k)
enddo
enddo

do j=-ighost+1,ny+ighost-1
do k=0,3
y=minmod(dgl(i,j,k),dgl(i,j+1,k))
ggl(i,j,k)=gl(i,j,k)+0.5*y
enddo
enddo

do j=-ighost+1,ny+ighost-2
do k=0,3
y=minmod(dgr(i,j+1,k),dgr(i,j+2,k))
ggr(i,j,k)=gr(i,j+1,k)-0.5*y
enddo
enddo

do j=-ighost+1,ny+ighost-2
do k=0,3
gt(i,j,k)=ggl(i,j,k)+ggr(i,j,k)
enddo
enddo
   enddo
end

CFL=0.2*min(dx,dy)/dmaxvel
end

subroutine bound(u,nx,ny,ighost)
implicit real (a-h,o-z)
   real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)

   c=1.0/3.0
   do k=0,3
do j=-ighost,ny+ighost
do i=-ighost+1,-ighost,-1
u(i,j,k)=c*(4*u(i+1,j,k)-u(i+2,j,k))
enddo
do i=nx+1,nx+ighost
u(i,j,k)=c*(4*u(i-1,j,k)-u(i-2,j,k))
enddo
enddo
do i=-ighost,nx+ighost
do j=-ighost+1,-ighost,-1
u(i,j,k)=c*(4*u(i,j+1,k)-u(i,j+2,k))
enddo
do j=ny+1,ny+ighost
u(i,j,k)=c*(4*u(i,j-1,k)-u(i,j-2,k))
enddo
enddo
end do

do k=0,3
u(-1,nx+1,k)=u(0,nx,k)
u(-2,nx+ighost,k)=u(-1,nx+1,k)
u(-1,-1,k)=u(0,0,k)
u(-2,-2,k)=u(-1,-1,k)
u(nx+1,-1,k)=u(nx,0,k)
u(nx+ighost,-2,k)=u(nx+1,-1,k)
u(nx+1,nx+1,k)=u(nx,nx,k)
u(nx+ighost,nx+ighost,k)=u(nx+1,nx+1,k)
enddo

end

real function minmod(x,y)
   implicit real(a-h,o-z)
   if(x*y.le.0.) then
minmod = 0.
   elseif(abs(x).ge.abs(y)) then
minmod = y
   else
minmod = x
   endif
   return
   end

real function CFL(u,dx,dy,nx,ny,ighost,gama)
   implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)

   dmaxvel=1e-10
do i=1,nx
do j=1,ny
uu=u(i,j,1)/u(i,j,0)
vv=u(i,j,2)/u(i,j,0)
p=(gama-1)*(u(i,j,3)-0.5*u(i,j,0)*(uu**2+vv**2))
vel=sqrt(gama*p/u(i,j,0))+sqrt(uu**2+vv**2)
if(vel.gt.dmaxvel) dmaxvel=vel
enddo
enddo

subroutine output(u,nx,ny,dx,dy,ighost,gama)
implicit real(a-h,o-z)
real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)

   open(unit=1,file='result.plt')
   write(1,*)'TITLE = "EXAMPLE: 3D GEOMETRIES" '
   write(1,*)'VARIABLES = "X ", "Y  ", "rou" "u" "v" "p" "E"'
   write(1,*)'ZONE T="Floor", I=',nx,' J=',ny,' F=POINT '
do i=1,nx
do j=1,ny
rou=u(i,j,0)
uu=u(i,j,1)/rou
vv=u(i,j,2)/rou
p=(gama-1)*(u(i,j,3)-0.5*rou*(uu*uu+vv*vv))
write(1,*) i*dx,j*dy,rou,uu,vv,p,u(i,j,3)
enddo
enddo
   close(1)

end

可以和我交流，QQ121493418

通流 · 发表于 2011-4-27 09:22:31

我们在数值解的时候是没有引入等熵假设，因为根本就不需要，欧拉方程组中的能量方程本来就是等熵方程，参考流体力学基础，清华大学，潘文全，144页
------------------------------------------------------------------------------------------
是看错了还是书写错了？

		自动登录	找回密码
密码			注册

解欧拉方程为什么会得到激波？

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。