求J.D.Aderson 计算流体力学入门7.3节采用预报-校正法求解拟一维喷管数值解源程序，最好是C！！求大神解

最新问答

• 心菲殿下

  你给个邮箱啊。我这有。但Aderson的源程序是fortran写的。贴这里得了
  Program nozzle
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  C A Program of quasi-one dimensional nozzle flows
  C Written by Dong Gang, May, 2003
  C
  C Based on section 7.3&7.4, chapter 7, Computational Fluid
  C Dynamics, The Basics with Applications,
  C John D. Anderson, JR. McGraw-Hill, 2002, 4
  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  parameter (n=31)
  dimension x(n),a(n),r(n),t(n),u(n),p(n),
  + am(n),fm(n),rm(n),um(n),tm(n)
  dimension delt(n),drdt(n),dudt(n),dtdt(n),drdtm(n),dudtm(n),
  + dtdtm(n),drdtav(n),dudtav(n),dtdtav(n)
  dimension rt(10000),tt(10000),pt(10000),amt(10000)
  gamma=1.4
  open(11,file='output\output.dat',status='unknown')
  open(12,file='output\outputa.dat',status='unknown')
  write(*,*)'please input the nstop'
  read (*,*)nstop
  5 write(*,*)'Subsonic-Supersonic/Purely subsonic(0/1)?'
  read (*,*)icase
  C
  C Initial conditions specified
  C
  if (icase.eq.0) then
  do i=1,n
  x(i)=(i-1)/((n-1)/3.)
  a(i)=1.+2.2*(x(i)-1.5)*(x(i)-1.5)
  r(i)=1.-0.3146*x(i)
  t(i)=1.-0.2314*x(i)
  u(i)=(0.1+1.09*x(i))*(t(i)**0.5)
  end do
  elseif (icase.eq.1) then
  write(*,*)'In purely subsonic case, input the exit pressure,'
  write(*,*)'the range is from 0.528P0 to P0'
  read(*,*) pe
  do i=1,n
  x(i)=(i-1)/((n-1)/3.)
  a(i)=1.+0.2223*(x(i)-1.5)*(x(i)-1.5)
  if (x(i).le.1.5) a(i)=1.+2.2*(x(i)-1.5)*(x(i)-1.5)
  r(i)=1.-0.023*x(i)
  t(i)=1.-0.009333*x(i)
  u(i)=0.05+0.11*x(i)
  end do
  else
  write(*,*)'There is no such case, repeat now!'
  goto 5
  end if
  C
  C Time stepping and looping
  C
  deltx=3./real(n-1)
  cfl=0.5
  do ii=1,nstop
  dt=1.
  do i=1,n
  as=sqrt(t(i))
  delt(i)=cfl*deltx/(as+u(i))
  if (dt.gt.delt(i)) dt=delt(i)
  end do
  C
  C Get solutions by MacCormack's explicit technique
  C
  do i=2,n-1
  drdt(i)=(-r(i)*(u(i+1)-u(i))-r(i)*u(i)*(log(a(i+1))-log(a(i)))
  + -u(i)*(r(i+1)-r(i)))/deltx
  dudt(i)=(-u(i)*(u(i+1)-u(i))-1./gamma*(t(i+1)-t(i)+
  + t(i)/r(i)*(r(i+1)-r(i))))/deltx
  dtdt(i)=(-u(i)*(t(i+1)-t(i))-(gamma-1.)*t(i)*(u(i+1)-u(i)+
  + u(i)*(log(a(i+1))-log(a(i)))))/deltx
  C
  rm(1)=r(1)
  um(1)=u(1)
  tm(1)=t(1)
  rm(i)=r(i)+drdt(i)*dt
  um(i)=u(i)+dudt(i)*dt
  tm(i)=t(i)+dtdt(i)*dt
  C
  drdtm(i)=(-rm(i)*(um(i)-um(i-1))-rm(i)*um(i)*(log(a(i))-
  + log(a(i-1)))-um(i)*(rm(i)-rm(i-1)))/deltx
  dudtm(i)=(-um(i)*(um(i)-um(i-1))-1./gamma*(tm(i)-tm(i-1)
  + +tm(i)/rm(i)*(rm(i)-rm(i-1))))/deltx
  dtdtm(i)=(-um(i)*(tm(i)-tm(i-1))-(gamma-1.)*tm(i)*
  + (um(i)-um(i-1)+um(i)*(log(a(i))-log(a(i-1)))))/deltx
  C
  drdtav(i)=(drdt(i)+drdtm(i))/2.
  dudtav(i)=(dudt(i)+dudtm(i))/2.
  dtdtav(i)=(dtdt(i)+dtdtm(i))/2.
  C
  r(i)=r(i)+drdtav(i)*dt
  u(i)=u(i)+dudtav(i)*dt
  t(i)=t(i)+dtdtav(i)*dt
  p(i)=r(i)*t(i)
  am(i)=u(i)/sqrt(t(i))
  fm(i)=r(i)*a(i)*u(i)
  if (i.eq.15) then
  rt(ii)=r(i)
  tt(ii)=t(i)
  pt(ii)=p(i)
  amt(ii)=am(i)
  end if
  end do
  C
  C Boundary conditions at inflow (n=1)
  C
  u(1)=2.*u(2)-u(3)
  r(1)=1.
  t(1)=1.
  p(1)=r(1)*t(1)
  am(1)=u(1)/sqrt(t(1))
  fm(1)=r(1)*a(1)*u(1)
  C
  C Boundary conditions at outflow (n=n)
  C
  if (icase.eq.0) then
  u(n)=2.*u(n-1)-u(n-2)
  r(n)=2.*r(n-1)-r(n-2)
  t(n)=2.*t(n-1)-t(n-2)
  p(n)=r(n)*t(n)
  am(n)=u(n)/sqrt(t(n))
  fm(n)=r(n)*a(n)*u(n)
  else
  p(n)=pe
  u(n)=2.*u(n-1)-u(n-2)
  t(n)=2.*t(n-1)-t(n-2)
  r(n)=p(n)/t(n)
  am(n)=u(n)/sqrt(t(n))
  fm(n)=r(n)*a(n)*u(n)
  end if
  C
  end do
  C
  C Print the results
  C
  do i=1,n
  write(11,100)x(i),a(i),r(i),u(i),t(i),p(i),am(i),fm(i)
  end do
  do ii=1,nstop
  write(12,101)ii,rt(ii),tt(ii),pt(ii),amt(ii)
  end do
  100 format(f5.2,7f8.3)
  101 format(i3,4f8.3)
  end

  浏览 233赞 92时间 2023-02-26