GROUP 1. Run title and other preliminaries
   PHOTON USE
   AUTOPLOT
   file
   phi 5
 
   cl
   msg MHD PLANE CHANNEL FLOW
   msg Hartmann number = 10 Reynolds number =100
   msg Velocity (W1) profile
   msg Blue line --- PHOENICS solution
   msg crosses ---   analytical solution
   da 1 w1;da 1 w1a
   col3 1;blb4 2
   msg press  to continue
   pause
 
   msg press  to end
   pause
   end
   END_USE
  DISPLAY
 
TEXT(1D Laminar MHD Channel Flow
TITLE
 
  DISPLAY
   This problem concerns the steady fully-developed laminar
   flow of an incompressible electrically-conducting fluid
   in the positive z-direction of a plane channel. A uniform
   magnetic field By is imposed normal to the walls and a
   current jx is induced in the fluid in the x-direction,
   together with a magnetic field Bz in the z-direction. The
   problem neglects end effects, secondary flows, Hall effect
   and ion-slip phenomena.
  ENDDIS
 
   The dimensionless momentum equation to be solved is:
 
   dp/dz + d/dy(dw/dy)/Re + Ha*Ha(K-w)/Re = 0
 
   where y =y/yin
         z =z/yin
         w =w/win
         p=p/(rho*win**2)
         Re=win*yin/enul
         Ha=sig*(By*yin)**2/(rho1*enul)
         K =Ex/(win*By)
 
   Here, yin is the channel half width (m), By is
   the imposed magnetic flux density in the +ve y
   direction (volt.s/m**2), sig is the electric
   conductivity (ohm/m) and Ex is electric field
   intensity in the +ve x-direction (volt/m).
 
   The Hartmann number Ha represents the ratio of
   the electromagnetic forces to the viscous forces.
 
   The voltage ratio K is the ratio of the voltage
   to the open-circuit voltage. If K=0, this corresponds
   to a short-circuit condition. When K=1 the net current
   flow is zero, which is known as the open-circuit
   condition. This is the classical Hartmann problem.
   Further, if K < 1 the channel will act as a MHD pump,
   whereas if K > 1 the channel will act as a flowmeter.
 
   The analytical solution to this problem has been
   presented in 'Engineering Magnetohydrodynamics',
   Chapter 10, G.W.Sutton and A.Sherman, McGraw Hill,
   (1965).
 
REAL(YIN,HIN,WIN);YIN=1.0;WIN=1.0
    GROUP 2. Transience; time-step specification
CARTES=T
    GROUP 4. Y-direction grid specification
NY=40;GRDPWR(Y,NY,YIN,1.)
    GROUP 5. Z-direction grid specification
ZWLAST=0.1
    GROUP 7. Variables stored, solved & named
SOLVE(W1);STORE(W1A)
    GROUP 8. Terms (in differential equations) & devices
TERMS(W1,N,N,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
RHO1=1.;ENUT=0.
REAL(REY,FLOWIN);REY=100.;ENUL=1./REY
REAL(HA,KVOLT);HA=10.0;KVOLT=1.0
HA
KVOLT
REY
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN
  ** compute analytical solutions
REAL(WA,GR,HAY);INTEGER(JJM1)
DO JJ=1,NY
+PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1)
+GR=0.5*YFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1
+GR=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1))
ENDIF
+GR=GR*YVLAST;HAY=GR*HA
+WA=HA*(COSH(HA)-COSH(HAY))/(HA*COSH(HA)-SINH(HA))
+INIT(IN:JJ:,W1A,ZERO,WA)
ENDDO
    GROUP 13. Boundary conditions and special sources
PATCH(WALL,NWALL,1,NX,NY,NY,1,NZ,1,1)
COVAL(WALL,W1,1.0,0.0)
  ** activate fully-developed pressure-drop calculation
FDFSOL=T;USOURC=T
FLOWIN=RHO1*WIN*YIN
PATCH(FDFW1DP,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(FDFW1DP,W1,FLOWIN,GRND1)
 
PATCH(MHDFOR,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(MHDFOR,W1,HA*HA/REY,KVOLT)
 
    GROUP 15. Termination of sweeps
LSWEEP=5;LITHYD=8
    GROUP 16. Termination of iterations
RESREF(W1)=1.E-12*WIN*YIN
    GROUP 17. Under-relaxation devices
REAL(DTF);DTF=10.0*ZWLAST/WIN
RELAX(W1,FALSDT,DTF)
    GROUP 19. Data communicated by satellite to GROUND
FLOWIN=RHO1*WIN*YIN
    GROUP 22. Spot-value print-out
IYMON=NY;TSTSWP=-1
    GROUP 23. Field print-out and plot control
NPLT=1;NYPRIN=1;NZPRIN=1;NTPRIN=2
    GROUP 24. Dumps for restarts