Library case Y303:

A THIN-SLAB CASTER WITH ELECTROMAGNETIC BRAKE



TALK=f;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 on the free surface plane
   0.1E+04 0.025E+04 CR
   pause

   ENDUSE

    GROUP 1. Run title and other preliminaries

TEXT(CSM - BFC CASTER with electromagnetic brake
TITLE

  DISPLAY

  Library case Y303:

  A THIN-SLAB CASTER WITH ELECTROMAGNETIC BRAKE

  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.
  A magnetic field is applied in the region of the mould. All its
  attributes (magnitude, frequency, and position) can be modified
  from the Q1.
  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


    **************************************************************
    *                PARAMETERS  OF  THE CASE                    *
    **************************************************************


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 DATA                    *
    ************************************************************

CARTES=T

    ***************** DIMENSIONS *******************************

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)

    *************** CELL DISTRIBUTION  **************************

     * 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,FX,FY,FZ,BX,BY,BZ,BS)
STORE(BX2,BY2,BZ2,B2,BXY,BXZ,BYZ,YP,ZP)
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.)

  ****************** MAGNETIC BRACKING  SETTING **********************

        Set constants
REAL(MU,BMAX,SIGMA)
SIGMA=7.22E5    ; MU=1.26E-6    ; BMAX=0.25

    **  General form of the magnetic field implemented,
        depending on the real variables A1, A2, A3, B, C, D, E.
        and on the integer variables    N1, M1, P1, K1
                                        N2, M2, P2, K2
                                        N3, M3, P3, K3
    Bx = A1 * SIN(B*y*N1)) * COS(C*y*M1) * SIN(D*z*P1) * COS(E*z*K1)
    By = A2 * SIN(B*y*N2)) * COS(C*y*M2) * SIN(D*z*P2) * COS(E*z*K2)
    Bz = A3 * SIN(B*y*N3)) * COS(C*y*M3) * SIN(D*z*P3) * COS(E*z*K3)

    where you can adjust    A1 changing RG(36)
                            A2 changing RG(40)
                            A3 changing RG(41)
                            B  changing RG(35)
                            C  changing RG(39)
                            D  changing RG(33)
                            E  changing RG(38)
                            N1 changing IG(21)
                            M1 changing IG(22)
                            P1 changing IG(23)
                            K1 changing IG(24)
                            N2 changing IG(25)
                            M2 changing IG(26)
                            P2 changing IG(27)
                            K2 changing IG(28)
                            N3 changing IG(29)
                            M3 changing IG(30)
                            P3 changing IG(31)
                            K3 changing IG(32)
     To get rid of the effect of a cosine, set the corresponding M
     or K constant to zero: then cos(...)=1.
     To get rid of the effect of a sine, set the corresponding
     constant N or P to zero: this will replace sin(...) by 1.

    ** Settings corresponding to B such as:
                Bx = 0
                By = 0
                Bz = Bmax * sin(2*pi/largh*y)

RG(33)=1.
RG(35)=2*PI/LARGH
RG(36)= 0.
RG(37)=SIGMA
RG(38)=1.
RG(39)=1.
RG(40)= 0.
RG(41)= BMAX


IG(21)=0  ; IG(22)=0  ; IG(23)=0  ; IG(24)=0
IG(25)=0  ; IG(26)=0  ; IG(27)=0  ; IG(28)=0
IG(29)=1  ; IG(30)=0  ; IG(31)=0  ; IG(32)=0
     ** Settings corresponding to the first Magnetic field B you sent us

    RG(33)=PI/SPES
    RG(35)=2*PI/LARGH
    RG(36)=BMAX/1.732
    RG(40)=BMAX/1.732
    RG(41)=BMAX/1.732
    RG(37)=SIGMA
    RG(38)=PI/2
    RG(39)=1

    IG(21)=1  ; IG(22)=0  ; IG(23)=0  ; IG(24)=0
    IG(25)=1  ; IG(26)=0  ; IG(27)=0  ; IG(28)=0
    IG(29)=1  ; IG(30)=0  ; IG(31)=1  ; IG(32)=0

    ** Set the limits in x direction of the patch:
            first cell in x=IXMAGF
             last cell in x=IXMAGL
INTEGER(IXMAGF,IXMAGL)
IXMAGF=1
IXMAGL=10

    * Send these values to the ground coding
IG(11)=IXMAGF
IG(12)=IXMAGL

    ** Patch to set the electromagnetic force
PATCH(FMAG,VOLUME,IXMAGF,IXMAGL,1,NY,1,NZ,1,1)
COVAL(FMAG,U1,FIXFLU,GRND)
COVAL(FMAG,V1,FIXFLU,GRND)
COVAL(FMAG,W1,FIXFLU,GRND)

  ***** 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=50 
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
STOP