subroutine hyperbolic ********************************************************************** * * * This subroutine is to solve the * * hyperbolic equation of the density n(i,j). * * * * Version 0.5 - Last Modified 5/18/04 * * * ********************************************************************** INCLUDE "fpvar.inc" dimension nt(ST:NZ,ST:NX) re=0.05 eps0=1.0 epmax=1.0e18 dxz=dx*dz ********************************************************************** * * * Determine the time step 'dt' by first finding the * * maximum velocity in any given cell in either the x- * * or z-direction. * * * ********************************************************************** DT=0.0 310 EPS1=0.0 DO 312, I=3,MOZ-2 DO 312, J=1,MOX VIX1=ABS(F(I,J)/(N(I,J)*DX)) VIZ1=ABS(G(I,J)/(N(I,J)*DZ)) 312 DT=MAX(DT,VIX1,VIZ1) DT=0.35/DT IF (DT.GT.10.) DT=10. IF (DT.LT.1.) DT=1. ********************************************************************** * * * Calculate the time advanced solution for lower order flux. * * * ********************************************************************** DO 322, I=ST+2,NZ-2 AVEN=0.0 DO 320, J=1,MOX-1 320 AVEN=AVEN+N(I,J) AVEN=AVEN/REAL(MOX-1) do 322, j=ST+2,nx-2 dntd(i,j)=n(i,j)+dt*(produ(i)-recom(i)*n(i,j)) $ -(FL(I,J)-FL(I,J-1)+GL(I,J)-GL(I-1,J))/DXZ $ -dt*((cosi*xk-sini*zkr(i))**2)*(n(i,j)-aven) $ *ti(i)/(1.0+rvi(i)**2) 322 continue ********************************************************************** * * * Calculate the time advanced solution for higher order flux * * in the inner region (i=3:moz-2, j=1:mox). * * * ********************************************************************** do 330, i=3,moz-2 do 330, j=1,mox nt(i,j)=dntd(i,j)-(c1(i,j)*af(i,j)-c1(i,j-1)*af(i,j-1) $ +c2(i,j)*ag(i,j)-c2(i-1,j)*ag(i-1,j))/dxz 330 eps1=max(eps1,abs(nt(i,j)-n(i,j))) ********************************************************************** * * * Put the values of nt(i,j) to n(i,j). * * * ********************************************************************** do 340, i=3,moz-2 do 340, j=1,mox 340 n(i,j)=nt(i,j) ********************************************************************** * * * Calculate the values of n(2,j), n(1,j), n(moz-1,j), and * * n(moz,j) for j=1,mox under suitable boundary conditions. * * These values are not calculated in the calculations for * * the inner region. * * * ********************************************************************** do 350, j=1,mox n(2,j)=n(4,j)-(cosi/sini)*(n(3,j+1)-n(3,j-1))*dz/dx n(1,j)=n(3,j)-(cosi/sini)*(n(2,j+1)-n(2,j-1))*dz/dx n(moz-1,j)=n(moz-3,j) $ +(cosi/sini)*(n(moz-2,j+1)-n(moz-2,j-1))*dz/dx n(moz,j)=n(moz-2,j) $ +(cosi/sini)*(n(moz-1,j+1)-n(moz-1,j-1))*dz/dx 350 continue ********************************************************************** * * * Restrict sudden change of n(i,j) caused by numerical error. * * * ********************************************************************** do 360, i=1,moz do 360, j=1,mox if (n(i,j).lt.0.0) then n(i,j)=0.25*(n(i+1,j)+n(i-1,j)+n(i,j+1)+n(i,j-1)) end if 360 continue ********************************************************************** * * * Periodic boundary condition in the horizontal direction. * * * ********************************************************************** 370 do 372, j=-5,0 do 372, i=1,moz 372 n(i,j)=n(i,mox+j-1) do 374, j=mox+1,nx do 374, i=1,moz 374 n(i,j)=n(i,j-mox+1) open (status='old',access='append',unit=8,file='trouble.dat') write(8,*)'end of hyperbolic' close(8) return end ********************************************************************** *////////////////////////////////////////////////////////////////////* *\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\* *////////////////////////////////////////////////////////////////////* ********************************************************************** ********************************************************************** ** ** ** The functions used above for calculating the ** ** higher-order flux. ** ** ** ********************************************************************** ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * * * * ********************************************************************** REAL function c1(i,j) if (af(i,j).ge.0.0) then c1=min(r1(i,j+1),r0(i,j)) else c1=min(r1(i,j),r0(i,j+1)) end if return end ********************************************************************** REAL function c2(i,j) if (ag(i,j).ge.0.0) then c2=min(r1(i+1,j),r0(i,j)) else c2=min(r1(i,j),r0(i+1,j)) end if return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * * * * ********************************************************************** REAL function r1(i,j) INCLUDE "fpvar.inc" p(i,j)=max(0.0,ag(i-1,j))-min(0.0,ag(i,j))+ $ max(0.0,af(i,j-1))-min(0.0,af(i,j)) q(i,j)=(wmax(i,j)-dntd(i,j))*dxz if (p(i,j).gt.0.0) then r1=min(1.0,q(i,j)/p(i,j)) else r1=0.0 end if return end ********************************************************************** REAL function r0(i,j) INCLUDE "fpvar.inc" p(i,j)=max(0.0,ag(i,j))-min(0.0,ag(i-1,j))+ $ max(0.0,af(i,j))-min(0.0,af(i,j-1)) q(i,j)=(dntd(i,j)-wmin(i,j))*dxz if (p(i,j).gt.0.0) then r0=min(1.0,q(i,j)/p(i,j)) else r0=0.0 end if return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * Set of Velocity*Density Calculation routines: * * * * F is for the x-direction for a given cell, * * G is for the z-direction * * * ********************************************************************** REAL function f(i,j) INCLUDE "fpvar.inc" vix1=uebx+rix1(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ -ri2(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*sini/((1.0+rvi(i)**2)*b) $ +sini*cosi*gg/((1.0+rvi(i)*rvi(i))*vin(i)) $ -rix1(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ +ri2(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) f=vix1*n(i,j) return end ********************************************************************** REAL function g(i,j) INCLUDE "fpvar.inc" viz1=-ri2(i)*(uox+u1x(i,j)+eox/(rvi(i)*b)) $ +riz1(i)*(uoz+u1z(i,j)+eoz/(rvi(i)*b)) $ -eoy*cosi/((1.0+rvi(i)**2)*b) $ -(rvi(i)**2+sini**2)*gg/((1.0+rvi(i)**2)*vin(i)) $ +ri2(i)*(phi(i,j+1)-phi(i,j-1))/(2.0*dx*rvi(i)*b) $ -riz1(i)*(phi(i+1,j)-phi(i-1,j))/(2.0*dz*rvi(i)*b) g=viz1*n(i,j) return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * * * * ********************************************************************** REAL function FL(I,J) INCLUDE "fpvar.inc" fl=(0.5*dt*(f(i,j+1)+f(i,j))-0.25*dx*(n(i,j+1)-n(i,j)))*dz return end ********************************************************************** REAL function GL(I,J) INCLUDE "fpvar.inc" gl=(0.5*dt*(G(i+1,j)+G(i,j))-0.25*dz*(N(i+1,j)-n(i,j)))*dx return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * * * * ********************************************************************** REAL FUNCTION fh(I,J) INCLUDE "fpvar.inc" fh=(7.0*(F(i,j+1)+F(i,j))-(F(i,j+2)+F(i,j-1)))*DT*dz/12.0 RETURN END ********************************************************************** REAL function gh(i,j) INCLUDE "fpvar.inc" gh=(7.0*(g(i+1,j)+g(i,j))-(g(i+2,j)+g(i-1,j)))*dt*dx/12.0 return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * * * * ********************************************************************** REAL function ag(i,j) INCLUDE "fpvar.inc" ag=gh(i,j)-gl(i,j) t1=ag*(dntd(i+1,j)-dntd(i,j)) t2=ag*(dntd(i+2,j)-dntd(i+1,j)) t3=ag*(dntd(i,j)-dntd(i-1,j)) if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) ag=0.0 return end ********************************************************************** REAL function af(i,j) INCLUDE "fpvar.inc" AF=FH(I,J)-FL(I,J) t1=af*(dntd(i,j+1)-dntd(i,j)) t2=af*(dntd(i,j+2)-dntd(i,j+1)) t3=af*(dntd(i,j)-dntd(i,j-1)) if((t1.lt.0.0).and.((t2.lt.0.0).or.(t3.lt.0.0))) af=0.0 return end ********************************************************************** *////////////////////////////////////////////////////////////////////* ********************************************************************** * * * Set of max/min functions designed to retrieve max/min value * * of n or dntd for (i,j) and neighbouring cells. * * * * Used in functions r0 and r1 * * * ********************************************************************** REAL FUNCTION WMAX(I,J) INCLUDE "fpvar.inc" WMAX=MAX(N(I-1,J),N(I,J),N(I,J-1),N(I,J+1),N(I+1,J), $ DNTD(I-1,J),DNTD(I,J),DNTD(I,J-1),DNTD(I,J+1),DNTD(I+1,J)) RETURN END ********************************************************************** REAL FUNCTION WMIN(I,J) INCLUDE "fpvar.inc" WMIN=MIN(N(I-1,J),N(I,J),N(I,J-1),N(I,J+1),N(I+1,J), $ DNTD(I-1,J),DNTD(I,J),DNTD(I,J-1),DNTD(I,J+1),DNTD(I+1,J)) RETURN END