Library case Y302:

STEEL SOLIDUFICATION IN A THIN-SLAB CASTER




TALK=T;RUN( 1, 1)

 sppnam=tundish  

   PHOTON USE
   p







   vi z
   up -x
   vec z 1 sh
   set vec ref
   0.41
   mag
   8
   0.17559E+04 0.29376E+04 CR
   text
   4

   Velocity vectors at symmetry plane
   0.14987E+04 0.27620E+04 CR
   pause
   cl
   vec z 6 sh
   text
   4

   Velocity vectors after the nozzle
   0.14987E+04 0.27620E+04 CR
   pause
   cl
   vec z 7 sh
   text
   4

   Velocity vectors on the free surface plane
   0.14987E+04 0.27620E+04 CR
   pause
   cl
   mag
   1
   con solf z 1 fi;0.001
   text
   4

   Solid fraction contours at symmetry plane
   0.1+04 0.025E+04 CR
   pause
   cl
   con solf z 6 fi;0.001
   text
   4

   Solid fraction contour after the nozzle
   0.1E+04 0.025E+04 CR
   pause
   cl
   con solf z 7 fi;0.001
   text
   4

   Solid fraction contours on the free surface plane
   0.1E+04 0.025E+04 CR
   pause
   cl
   con tem1 z 1 fi;0.001
   text
   4

   Temperature contours at symmetry plane
   0.1+04 0.025E+04 CR
   pause
   cl
   con tem1 z 6 fi;0.001
   text
   4

   Temperature contour after the nozzle
   0.1E+04 0.025E+04 CR
   pause
   cl
   con tem1 z 7 fi;0.001
   text
   4

   Temperature contours in the free surface plane
   0.1E+04 0.025E+04 CR
   pause
   ENDUSE

    GROUP 1. Run title and other preliminaries

TEXT(CSM BFC THIN-SLAB FUNNEL-SHAPED CASTER
TITLE

  DISPLAY

  Library case Y302:

  STEEL SOLIDIFICATION IN A THIN-SLAB CASTER

  The case considered is 3D steady, turbulent flow and heat transfer
  with solidification in the mould and sub-mould regions of a
  vertical steel slab caster with a funnel-shaped mould. Liquid
  steel with 30oC superheat is introduced into the mould through a
  submerged entry nozzle, and the strand is continuously withdrawn
  at a fixed speed of 0.6m/min. Turbulence is accounted for by means
  of the 2-equation k-e model. The fixed-grid enthalpy-porosity
  technique is used to represent the solidification.

  ENDDIS

    GROUP 1. Run title and other preliminaries

BOOLEAN(SOLID,CONCP,CONTK,KEMOD,RECTCAS)
       ** SOLID     = T (solidification activated)
                    = F (no solidification)
          CONCP     = T (constant specific heat)
                    = F (variable Cp)
          CONTK     = T (constant thermal conductivity )
                    = F (variable k)
          KEMOD     = T (k-e turbulence model)
                    = F (uniform turbulent viscosity)
          SHORTCAS  = T (lenth of the casting: 7m)
                    = F (length of the casting: 12 m)
          RECTCAS   = T (cartesian rectangular casting, with a
                         thickness defined by LZCAST)
                    = F (funnel)
SOLID=T;CONCP=F;CONTK=F;KEMOD=T;RECTCAS=F
REAL(TMPIN,RHOL,PI,HIN,RADIUS,ROLVEL,TLIQ,TSOL);PI=3.14159

       ** RHOL  = Melt density (kg/m^3)
          TLIQ  = Liquidus temperature (K)
          TSOL  = Solidus temperature (K)
          TMPIN = Melt inlet temperature (K)
RHOL=7300.0
ROLVEL=0.
TLIQ=1803.
TSOL=1753.
TMPIN=1833.
RADIUS=0.

   ** inlet enthalpy calculation, HIN
REAL(ACP,BCP,CCP,CPIN);ACP=0.136;BCP=273.15;CCP=523.0
IF(CONCP) THEN
+ HIN=700.0*TMPIN
ELSE
+ CPIN=0.5*ACP*TMPIN-ACP*BCP+CCP;HIN=CPIN*TMPIN
CPIN
ENDIF
HIN

    ** Mass flow rate through 1/4 of the nozzle area
REAL(GMASIN,VELCOL,SPES,LARGH,CASTLEN,DHYD,REY)
REAL(ANOZ,DNOZ,CLNOZ,VELAL)

    ** VELCOL  is the casting speed (m/s)
     * SPES    is the z-thickness of the mould (m)
     * LARGH   is the y-width of the mould (m)
     * CASTLEN is the x-length of the casting (m)
     * DNOZ    is the inlet nozzle diameter (m)
VELCOL=5./60.
SPES=.06
LARGH=1.6
CASTLEN=8.211
DNOZ=116.E-03

    ** Mass outflow rate
GMASIN=VELCOL*SPES*LARGH*RHOL
    ** Exploit 1/4 symmetry
GMASIN=0.25*GMASIN

    ** Nozzle Geometry: Area available to flow = 1/8 of nozzle area
ANOZ=0.25*PI*DNOZ*DNOZ
    ** Exploit 1/4 symmetry
ANOZ =0.25*ANOZ
    ** Equivalent Cartesian length of nozzle
CLNOZ=ANOZ**0.5

    ** Melt inlet velocity (m/s)
VELAL=GMASIN/(RHOL*ANOZ)

    ************************************************************
    *           Geometrical Settings of the Casting            *
    ************************************************************

CARTES=T
    *** Grid specifications

INTEGER(NX_NO1,NX_NO2,NX_MLD,NX_ZN1,NX_ZN2,NX_ZN3,NX_ZN4)
INTEGER(NY_NOZ,NY_F2,NY3,NY4,NZ_NOZ,NZ_2)
REAL(LMOLD,LFUNL,LZCART,Y_NOZ,Z_NOZ)
REAL(X_NO1,X_NO2,LZN1,LZN2,LZN3,LZN4,LMLD,LF2,LY3,LY4,LZ2)
INTEGER(NX_ZN5,NX_ZN6)
REAL(LZN5,LZN6)
   Geometrical data
    *  *_NO*  Dimensions of the nozzle
    *  LZN*   Length of zone *
    *  LMLD   Length of the mould after the nozzle
    *  LMOLD  Mould length (m)
    *  LMLD   Length of the mould after the nozzle
    *  LFUNL  Funnel length in Y direction.
              divided into two zones: Nozzle (Region 1)
              & LF2 (Region 2)
    *  LY3    Length of the region after the funnel in
              Y direction (Region 3)
    *  LY4    Last thin cells in Y direction (region 4)
    *  LZCART is the dimension in Z of a rectangular cartesian
              casting LZCART must be greater than Z_NOZ

LMOLD=1.050; LFUNL=0.6; LZCART= 0.07
X_NO1= 0.1818 ; X_NO2= 0.0909 ;Y_NOZ = CLNOZ ;Z_NOZ = CLNOZ
LMLD = LMOLD-(X_NO1+X_NO2)
LZN1=0.0494;LZN2=0.9191;LZN3=1.872;LZN4=1.2155;LZN5=1.235;LZN6=1.87

LF2 = LFUNL-Y_NOZ; LY3 =  0.5*LARGH- LFUNL- 0.002 ;LY4 = 0.002
LZ2 =  LZCART - (Z_NOZ)

     * Number of cells in X direction
         Nozzle area      : nx_NO1, nx_NO2
         Mould  (after the nozzle,up to end of mould)  : nx_mld
         Zone 1 (LZN1)    : nx_zn1
         Zone 2 (LZN2)    : nx_zn2
         Zone 3 (LZN3)    : nx_zn3
         Zone 4 (LZN4)    : nx_zn4
         Zone 5 (LZN5)    : nx_zn5
         Zone 6 (LZN6)    : nx_zn6

NX_NO1= 5; NX_NO2= 3; NX_MLD= 10
NX_ZN1=1;NX_ZN2=5;NX_ZN3=8;NX_ZN4=5;NX_ZN5=4;NX_ZN6=5
NX=NX_NO1+NX_NO2+NX_MLD+NX_ZN1+NX_ZN2+NX_ZN3+NX_ZN4
NX=NX+NX_ZN5+NX_ZN6

     * Number of cells in Y direction
         Nozzle area                     : nY_NOZ
         Funnel region, after the nozzle : NY_F2
         Region 3, after the funnel      : NY3
         Region 4, thin cells            : NY4

NY_NOZ= 5; NY_F2= 15; NY3=12; NY4=2
NY    = NY_NOZ+ NY_F2+ NY3+NY4

     * Number of cells in Z direction
         Nozzle area        : nZ_NOZ
         Rest of the domain : nz_2
NZ_NOZ= 5; NZ_2= 2; NZ= NZ_NOZ+ NZ_2

RSET(M,NX,NY,NZ)
XSI= CASTLEN; YSI= 0.5*LARGH; ZSI= LZCART; RSET(D,CHAM)

   GROUP 2. Transience; time-step specification
STEADY=T

   GROUP 3  X-direction grid specification
    * 9 Regions are set in X-direction:
          - The 2 first to set the nozzle
          - The 3rd for the Mould
          - The 4th to 9th to set the 6 zones for the temperature
            boundary conditions
       For each region, the length and number of cells is defined
       above and can be modified above.
REAL(GRPX3)
GRPX3=1.
NREGX=9
IREGX= 1; GRDPWR(X,NX_NO1,X_NO1,1.)
IREGX= 2; GRDPWR(X,NX_NO2,X_NO2,1.)
IREGX= 3; GRDPWR(X,NX_MLD,LMLD,:GRPX3:)
IREGX= 4; GRDPWR(X,NX_ZN1,LZN1,1.)
IREGX= 5; GRDPWR(X,NX_ZN2,-LZN2,1.3)
IREGX= 6; GRDPWR(X,NX_ZN3,LZN3,1.)
IREGX= 7; GRDPWR(X,NX_ZN4,LZN4,1.)
IREGX= 8; GRDPWR(X,NX_ZN5,LZN5,1.)
IREGX= 9; GRDPWR(X,NX_ZN6,LZN6,1.)

   GROUP 4  Y-direction grid specification
    * 4 Regions are set in Y-direction
          - The 1st, to set the nozzle
          - The 2nd, to define the end of the funnel
          - The 3rd, region after the funnel
          - The 4th, to refine the grid at the boundary of the
            domain
      As for X, the number of cells and lenth have to be
      modified above (if needed).

REAL(GRPY2)
GRPY2=1.4
NREGY= 4
IREGY= 1; GRDPWR(Y, NY_NOZ, Y_NOZ,  1.)
IREGY= 2; GRDPWR(Y,  NY_F2,   LF2,  :GRPY2:)
IREGY= 3; GRDPWR(Y,    NY3,   -LY3, .8)
IREGY= 4; GRDPWR(Y,    NY4,   LY4, 1.)

   GROUP 5  Z-direction grid specification
    * 2 Regions are set in Z-direction
        - The 1st, to set the nozzle
        - The 2nd,for the rest of the domain
NREGZ= 2
IREGZ= 1; GRDPWR(Z, NZ_NOZ, Z_NOZ, 1.)
IREGZ= 2; GRDPWR(Z,   NZ_2,   LZ2, 1.)

   GROUP 6. Body-fitted coordinates or grid distortion)

IF(.NOT.RECTCAS) THEN
+ BFC=T
+ INTEGER(II12,II13,JJ12,II12_,II13_,JJ12_,NZ_)
   * Imn is the number of cells in the region M to N, in I direction
   * Jmn is the number of cells in the region M to N, in J direction
   * Imn_ =Imn+1
   * Jmn_ =Jmn+1

+ II12= NX_NO1+ NX_NO2; II13= II12+ NX_MLD; JJ12= NY_NOZ+ NY_F2
+ II12_=II12+1; II13_=II13+1; JJ12_=JJ12+1;NZ_=NZ+1

    *** Thin zone of the casting where nz = 0.03 m
                  - Y = region 3 & 4
                  - X = region 4, 5, 6, 7, 8 & 9
      Rq: the grid power of the two following transfers have
          to stay identical!
+ GSET(C, K:NZ+1:, F, K1, II13_,  NX,    1, NY, +,0,0,0.03)
+ GSET(T, K:NZ+1:, F, K1, II13_,  NX,    1, NY,1.)
+ GSET(C, K:NZ+1:, F, K1,    1, II13, JJ12_, NY, +,0,0,0.03)
+ GSET(T, K:NZ+1:, F, K1,    1, II13, JJ12_, NY,1.)

   *** Funnel
+ GSET(T, J:JJ12_:, F, J:NY_NOZ+1:, 1, II12, 1, NZ_NOZ, :GRPY2:)
+ GSET(T, I:II13_:, F,    I:II12_:, 1, JJ12, 1, NZ_NOZ,:GRPX3:)
    *  Funnel curve
+ REAL(DX,DXB,D0,XMEN,X20,RX,ALFAX,BETA,ZMIN,REA1)
+ INTEGER(IMAX)
+ XMEN=0.08; D0=0.055; X20=0.75; DX=D0*(1-XMEN/X20); DXB=DX/LFUNL
+ ALFAX=2.*ATAN(DXB); RX=0.5*LFUNL/SIN(ALFAX);ZMIN=SPES/2.;IMAX=12
+ DO ii=1,IMAX+1
+  XPO=0.0; YPO=(II-1)/IMAX*LFUNL
+  IF(II.LE.(IMAX/2))THEN
+  REA1=YPO/RX; BETA=ASIN(REA1); ZPO=ZMIN+DX-RX*(1-COS(BETA))
+  ELSE
+  REA1=(LFUNL-YPO)/RX; BETA=ASIN(REA1); ZPO=ZMIN+RX*(1-COS(BETA))
+  ENDIF
+  GSET(P,M:II:)
+ ENDDO
+ XPO=XC(II13_,1,NZ_); YPO=YC(II13_,1,NZ_); ZPO=ZMIN; GSET(P,P1)
+ XPO=XC(II13_,JJ12_,NZ_);YPO=YC(II13_,JJ12_,NZ_)
+ ZPO=ZMIN;GSET(P,P2)
+ XPO=XC(II13_,:NY_NOZ+1:,NZ_);YPO=YC(II13_,:NY_NOZ+1:,NZ_)
+ ZPO=ZMIN;GSET(P,P3)
+ XPO=X_NO1+X_NO2; YPO=0.0; ZPO=ZMIN+DX*(1-XPO/LMOLD); GSET(P,Q1)
+ XPO=X_NO1+X_NO2; YPO=LFUNL; ZPO=ZMIN ; GSET(P,Q2)
    * lines & curves
+ GSET(V, CU1,M2,M3,M4,M5,M6,M7,M8,M9,M10,M11,M12,M13)
+ GSET(L, C1, M2, M13,  NY_F2, :GRPY2:, CRV, CU1)
+ GSET(L, L1, P1,  P3, NY_NOZ, 1.)
+ GSET(L, L2, P3,  P2,  NY_F2,:GRPY2:)
+ GSET(L, L3, Q1,  P1, NX_MLD, :GRPX3:)
+ GSET(L, L4, Q1,  M1,   II12, 1.)
+ GSET(L, L5, Q2,  P2, NX_MLD, :GRPX3:)
+ GSET(L, L6, Q2, M13,   II12, 1.)
+ GSET(L, L7, M1,  M2, NY_NOZ,1.)


   * Frame
+ GSET(F, F1, M1, Q1, P1, P3, P2, Q2, M13, M2)
+ GSET(M, F1, +I+J, 1, 1, NZ_,TRANS)
   * Transfer
+ GSET(T, k:NZ_NOZ+1:, F1, K:NZ_:, 1, :II13:, 1, :JJ12:, 1.0)
ENDIF

    GROUP 7. Variables stored, solved & named
    ** In addition to the known PHOENICS variables the following
       non-standard variables are defined:

    ** XDAR : VALue of the DARCY source term in the x-direction
    ** YDAR : VALue of the DARCY source term in the y-direction
    ** ZDAR : VALue of the DARCY source term in the z-direction
    ** FDAR : Variable used to determine whether DARCY source term
              is active or not (1 = active ; 0 = inactive)
    ** SOLF : Solid fraction of melt
    ** LATN : Latent heat content

SOLVE(P1);SOLUTN(P1,Y,Y,Y,N,N,N)
SOLVE(U1,V1,W1)
SOLVE(TEM1)

STORE(RHO1,LATN,SOLF,HSCE,XDAR,YDAR,ZDAR,SPH1,KOND)
STORE(PRPS,ENUT,FDAR,IMB1,EPOR,NPOR)
STORE(FX,FY,FZ,BX,BY,BZ,BXY,BYZ,BXZ,B2,BX2,BY2,BZ2)
STORE(HPOR,VPOR)

    GROUP 8. Terms (in differential equations) & devices
DIFCUT=0.0

   GROUP 9. Properties of the medium (or media)
    ** Constant density (= liquidus density)
RHO1=RHOL

    ** Set the thermal conductivity
IF(CONTK) THEN
+ PRNDTL(TEM1)=-25.0
ELSE
+ PRNDTL(TEM1)=-GRND
ENDIF
    ** Set specific heat
IF(CONCP) THEN
+ CP1=700.0
ELSE
+ CP1=GRND
ENDIF
    ** Laminar viscosity
REAL(TKEAL,EPSAL)
TKEAL=0.2**2*VELAL**2
EPSAL=TKEAL**1.5*0.1643/4.5E-3
ENUL=0.9E-06
REY=VELAL*DNOZ/ENUL
REY

    ** Activate k-e turbulence model
IF(KEMOD) THEN
+ TURMOD(KEMODL)
ELSE
+ DHYD=2.*(SPES*LARGH)/(SPES+LARGH)
+ ENUT=0.01*VELCOL*DHYD
ENUT
ENDIF
    ** LAMINAR PRANDTL NUMBER FOR H1
RG(1)=0.7
    ** LIQUID DENSITY
RG(2)=RHOL
    ** SOLID DENSITY
RG(3)=7300.0
    ** LIQUID SPECIFIC HEAT
RG(4)=700.0
    ** SOLID SPECIFIC HEAT
RG(5)=600.0
    ** LIQUID CONDUCTIVITY
RG(6)=30.0
    ** SOLID CONDUCTIVITY
RG(7)=25.0
    ** SOLIDUS TEMPERATURE
RG(8)=TSOL
    ** LIQUIDUS TEMPERATURE
RG(9)=TLIQ
    ** INDEX IN SOLID FRACTION EQUATION.
RG(10)=1.0
    ** LATENT HEAT OF SOLIDIFICATION
RG(11)=2.64E+05
    ** Limit over which DARCY source terms are active (from centre o
RG(12)=2.
    ** LINEAR RELAXATION IN ENTHALPY SOURCES.
RG(13)=0.1
    ** CONSTANT IN DARCY LAW (BIG A)
RG(14)=5.E+06
    ** SMALL NUMBER IN DARCY LAW TO PREVENT DIVISION BY ZERO
RG(15)=5.0E-6
    ** VALUE IN DARCY SOURCE TERMS (X,Y,Z)
RG(17)=VELCOL;RG(18)=0.0;RG(19)=0.0

   ????????????????????????????????????????????????????????
    ** ANGULAR VELOCITY OF ROLLERS
          RG(20) =  ROLVEL/RADIUS
   ????????????????????????????????????????????????????????

    ** FRACTION OF (TLIQ-TSOL) BY WHICH TEMPERATURE IS ALLOWED TO VA
RG(21)=1.0
    ** ROTATIONAL UNIT VECTOR COORDINATES  (X,Y,Z)
          RG(22) = 0.0 ; RG(23) = -1.0 ; RG(24) = 0.0
    ** LOWER LIMIT OF SOLID FRACTION FOR DARCY SOURCES
RG(25)=0.
    ** MINIMUM ALLOWABLE TEMPERATURE IN THE SYSTEM
RG(26)=300.0
    ** SPECIFIC HEAT COEFFICIENTS
RG(27)=ACP
RG(28)=BCP
RG(29)=CCP

    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
    ** Initial fields for KE and EP

FIINIT(KE)=TKEAL;FIINIT(EP)=EPSAL
FIINIT(ENUT)=0.09*TKEAL**2/EPSAL
FIINIT(W1)=0.0;FIINIT(U1)=0.0;FIINIT(V1)=0.0
FIINIT(FDAR)=0.0
FIINIT(TEM1)=TMPIN;FIINIT(RHO1)=RHOL
FIINIT(LATN)=RG(11);FIINIT(SOLF)=0.0;FIINIT(P1)=-9957.

IF(CONCP) THEN
+ FIINIT(SPH1)=700.
ELSE
+ FIINIT(SPH1)=CPIN
ENDIF

    ** Nozzle Blockages  #=whole region $=1st cell %=last cell
REAL(IZFA,IZFO,IZS,IXFI,IXFF,IXS,IYFA,IYFO,IYS)
CONPOR(0.0,HIGH,#1,#2,#1,#1,%1,%1)
CONPOR(0.0,NORTH,#1,#1,%1,%1,#1,#1)
CONPOR(0.0,WEST,%2,%2,#1,#1,#1,#1)
INIADD=F

    GROUP 12. Convection and diffusion adjustments
    GROUP 13. Boundary conditions and special sources
   ***** WALL FUNCTIONS *******
   ** HIGH WALL
PATCH(WALLH,HWALL,1,NX,1,NY,NZ,NZ,1,LSTEP)
COVAL(WALLH,V1,GRND2,0.0);COVAL(WALLH,U1,GRND2,0.0)
COVAL(WALLH,KE,GRND2,GRND2);COVAL(WALLH,EP,GRND2,GRND2)

   ** NORTH WALL
PATCH(WALLN,NWALL,1,NX,NY,NY,1,NZ,1,LSTEP)
COVAL(WALLN,U1,GRND2,0.0);COVAL(WALLN,W1,GRND2,0.0)
COVAL(WALLN,KE,GRND2,GRND2);COVAL(WALLN,EP,GRND2,GRND2)

  ***** INLET BOUNDARY *******
INTEGER(KK1)
REAL(CELLIN)
CELLIN=NZ_NOZ*NY_NOZ
KK1= NZ_NOZ-1
PATCH(NOCPFFL1,CELL,1,1,#1,#1,1,KK1,1,LSTEP)
COVAL(NOCPFFL1,P1,FIXFLU,GMASIN/CELLIN)
COVAL(NOCPFFL1,V1,ONLYMS,0.)
COVAL(NOCPFFL1,U1,ONLYMS,VELAL);COVAL(NOCPFFL1,TEM1,ONLYMS,HIN)
COVAL(NOCPFFL1,KE,ONLYMS,TKEAL);COVAL(NOCPFFL1,EP,ONLYMS,EPSAL)

PATCH(NOCPIN2,EAST,1,1,#1,#1,%1,%1,1,LSTEP)
COVAL(NOCPIN2,P1,FIXP,0.0)
COVAL(NOCPIN2,V1,ONLYMS,0.)
COVAL(NOCPIN2,U1,ONLYMS,VELAL);COVAL(NOCPIN2,TEM1,ONLYMS,HIN)
COVAL(NOCPIN2,KE,ONLYMS,TKEAL);COVAL(NOCPIN2,EP,ONLYMS,EPSAL)

  ***** OUTLET BOUNDARY *******

PATCH(OUTLFFL,WEST,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(OUTLFFL,P1,FIXFLU,-RHOL*VELCOL)
COVAL(OUTLFFL,TEM1,ONLYMS,SAME)
COVAL(OUTLFFL,KE,ONLYMS,SAME);COVAL(OUTLFFL,EP,ONLYMS,SAME)
COVAL(OUTLFFL,U1,ONLYMS,VELCOL)
COVAL(OUTLFFL,V1,ONLYMS,0.);COVAL(OUTLFFL,W1,ONLYMS,0.)

  ***** LATENT HEAT *******
IF(SOLID) THEN
+ PATCH(LATLIN,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(LATLIN,TEM1,GRND3,GRND3)

+ PATCH(LATFIX,CELL,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(LATFIX,TEM1,FIXFLU,GRND3)

  ***** DARCY SOURCES *****

+ PATCH(DARCY,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(DARCY,U1,GRND,VELCOL)
+ COVAL(DARCY,V1,GRND,0.0);COVAL(DARCY,W1,GRND,0.0)
+ COVAL(DARCY,KE,GRND,0.0);COVAL(DARCY,EP,GRND,0.0)


    ***** HEAT FLUX PATCHES ******

        * NMLD is the number of X-cells in the mould section
        * HLING is the length of the mould
        * XPS is the position in x-direction of the centre
        * HEATFL is the heat flux

INTEGER(NMLD)
REAL(HEATFL,XPS,HLING)
NMLD=NX_NO1+NX_NO2+NX_MLD
HLING=LMOLD
   ** Heat loss through mould walls
DO II=1,NMLD
IF (II.EQ.1) THEN
XPS=XULAST*XFRAC(1)/2
ELSE
XPS=XULAST*(XFRAC(II)+XFRAC(:II-1:))/2
ENDIF
HEATFL=-1.404E6*(HLING/XPS)**0.65

PATCH(MOULDH:ii: , HIGH,  ii, ii,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(MOULDH:ii: , TEM1,  FIXFLU,  HEATFL  )

PATCH(MOULDN:ii: , NORTH,  ii, ii,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(MOULDN:ii: , TEM1,  FIXFLU,  HEATFL  )

ENDDO
   ** Heat loss through high wall of zone 1
PATCH(SLABH1 , HIGH,  #4, #4,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH1 , TEM1, 4000., 400.)

   ** Heat loss through high wall of zone 2
PATCH(SLABH2 , HIGH,  #5, #5,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH2 , TEM1, 1700., 400.)

   ** Heat loss through high wall of zone 3
PATCH(SLABH3 , HIGH,  #6, #6,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH3 , TEM1, 1060., 300.)

   ** Heat loss through high wall of zone 4
PATCH(SLABH4 , HIGH,  #7, #7,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH4 , TEM1, 855., 300.)

   ** Heat loss through high wall of zone 5
PATCH(SLABH5 , HIGH,  #8, #8,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH5 , TEM1, 650., 300.)

   ** Heat loss through high wall of zone 6
PATCH(SLABH6 , HIGH,  #9, #9,  1, NY,  NZ, NZ,  1, LSTEP)
COVAL(SLABH6 , TEM1, 330., 300.)

   ** Heat loss through north wall of zone 1
PATCH(SLABN1, NORTH,  #4, #4,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN1,  TEM1, 4000., 400.)


   ** Heat loss through north wall of zone 2
PATCH(SLABN2, NORTH,  #5, #5,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN2,  TEM1, 1700., 400.)

   ** Heat loss through north wall of zone 3
PATCH(SLABN3, NORTH,  #6, #6,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN3,  TEM1, 1060., 300.)

   ** Heat loss through north wall of zone 4
PATCH(SLABN4, NORTH,  #7, #7,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN4,  TEM1, 855., 300.)

   ** Heat loss through north wall of zone 5
PATCH(SLABN5, NORTH,  #8, #8,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN5,  TEM1, 650., 300.)

   ** Heat loss through north wall of zone 6
PATCH(SLABN6, NORTH,  #9, #9,  NY, NY,  1, NZ,  1, LSTEP)
COVAL(SLABN6,  TEM1, 330., 300.)


   ** Free Surface
+ PATCH(HEATSUR1,WEST,1,1,1,NY,#2,#2,1,LSTEP)
+ COVAL(HEATSUR1,TEM1,FIXFLU,-7.000E+03)

+ PATCH(HEATSUR2,WEST,1,1,#2,#4,#1,#1,1,LSTEP)
+ COVAL(HEATSUR2,TEM1,FIXFLU,-7.000E+03)
ENDIF

    GROUP 14. Downstream pressure for PARAB=.TRUE.
    GROUP 15. Termination of sweeps
LSWEEP= 100
REAL(MASREF);MASREF=1.E-12*GMASIN
  ** Deactivate built-in residual calculation
SELREF=F
RESREF(P1)=MASREF/RHOL
RESREF(U1)=MASREF*VELAL
RESREF(V1)=RESREF(U1);RESREF(W1)=RESREF(U1)
RESREF(KE)=MASREF*TKEAL;RESREF(EP)=MASREF*EPSAL
RESREF(TEM1)=MASREF*HIN

    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices

REAL(DTF)
IF(CONCP) THEN
+ DTF=0.1*CASTLEN/VELAL/NX
ELSE
+ DTF=0.05*CASTLEN/VELAL/NX
ENDIF
DTF
RELAX(P1,LINRLX,1.0)
RELAX(V1,FALSDT,1.e-3);RELAX(U1,FALSDT,1.e-3);RELAX(W1,FALSDT,1.e-3)
IF(KEMOD) THEN
+ RELAX(KE,FALSDT,0.2*DTF);RELAX(EP,FALSDT,0.2*DTF);KELIN=3
ENDIF
RELAX(TEM1,FALSDT,50.*DTF)

    GROUP 18. Limits on variables or increments to them
    GROUP 19. Data communicated by satellite to GROUND

    ** Set LG(1)=T to update latent heats & solid fractions at IZ
    ** Set LG(1)=F to update latent heats & solid fractions at sweep
LG(1)=T
    Set LG(3)=T to impose limits on temperature increments
LG(3)=T
    ** Set LG(4)=T for no solidification
                =F for solidification
IF(SOLID) THEN
+ LG(4)=F
ELSE
+ LG(4)=T
ENDIF

    GROUP 20. Preliminary print-out
    GROUP 21. Print-out of variables
IXMON=2;IYMON=NY/2;IZMON=NZ-2
INIFLD=F
OUTPUT(P1,Y,N,Y,Y,Y,Y);OUTPUT(SOLF,Y,N,Y,Y,Y,Y)
   ** Built-in facility to dump PHI files every IDISPA sweeps
CSG1=SW;IDISPA=10000

    GROUP 22. Spot-value print-out
    GROUP 23. Field print-out and plot control

NPRINT=LSWEEP

ITABL=3;NPLT=5;TSTSWP=-1
    GROUP 24. Dumps for restarts

   RESTRT(ALL)
DISTIL=T
+ EX(P1  )=1.794E+04;EX(U1  )=9.022E-02;EX(V1  )=1.963E-02
+ EX(EPKE)=1.770E+00;EX(VPOR)=1.000E+00;EX(HPOR)=9.964E-01
+ EX(NPOR)=9.976E-01;EX(EPOR)=9.982E-01;EX(FDAR)=1.000E+00
+ EX(KOND)=2.981E+01;EX(SPH1)=5.917E+02;EX(SOLF)=6.407E-01
+ EX(LATN)=9.486E+04;EX(RHO1)=7.300E+03
+ EX(TEM1)=1.574E+03;EX(VCRT)=1.880E-02;EX(UCRT)=9.007E-02
nx
ny
nz
STOP