TALK=T;RUN( 1, 1) s124

  photon use
  p
 
 
 
  gr x 1
  msg Press RETURN to continue
  pause
  msg temperature contours in solid and fluid
  con tem1 x 1 fi;0.05
  dump tem1.gif
  msg Press RETURN to continue
  pause
  msg contours of horizontal displacement
  con w1 x 1 fi;0.05
  dump w1.gif
  msg Press RETURN to continue
  pause
  msg outline of solid
  set property off
  con prps x 1;val 1; 100
  msg displacements in solid
  vec x 1 y 1 10 z 6 15 co 2
  dump disp.gif
  msg Press RETURN to continue
  pause
  msg velocity vectors in the fluid
  vec x 1 y 1 m z 1 5  co 1
  vec x 1 y 11 m z 6 15 co 1
  vec x 1 y 1 m z 16 20 co 1
  dump vel.gif
  msg Press RETURN to continue
  pause
  con off
  gr off
  red
  gr ou x 1
  msg  radial stress contours
  con stry x 1 y 1 10 z 6 15 fi;0.001
  msg Press RETURN to continue
  msg  radial strain contours
  con epsy x 1 y 1 10 z 6 15 fi;0.001
  msg Press RETURN to continue
  pause
  con off
  gr off
  red
  msg vectors only
  con prps x 1;val 1; 2
  msg Press RETURN to continue
  pause
  msg radial strains
  con epsy x 1 y 1 10 z 6 15 fi;0.001
  dump epsy.gif
  msg
  
  msg Press e to END
  ENDUSE
#cls
TEXT(Centre-Heated, Edge-Cooled Block
TITLE
  DISPLAY
  STRESS ANALYSIS IN SOLIDS - Centre-heated, edge-cooled block
 
  2-dimensional (y-z), cartesian, steady, elliptic simulation
 
 
  The problem simulated is sketched below. Metal block is heated at
  the centre and externally cooled. The centre is maintained at
  1500 K. Cool external gas flows through a porous media. Thermal
  stresses which arise from temperature difference cause metal
  dilatation.
 
                      ext. load
             +----------|---------+
             |          v         -->
             |       +-----+           outflow
    inflow   |     ->| bar |<-     -->               y ^
         -->         |     |                          |
             +--------------------+                   +---->
                                                           z
  ENDDIS
#pause
  STOREd variables are as follows:
  STRX Stress distribution in axial direction
  STRY Stress distribution in radial direction
 
 
REAL(WIN,PI)
REAL(RESCO,TIN,heatsor)
INTEGER(IYNORT,IZLOW,IZHI,UNIT)
UNIT=5
  unit=2
  UNIT=1
heatsor=1.0e5
  heatsor=0.0
  ** switch for stress & strain post-processing
CALSTR=T
  ** porous media resistance coeff
RESCO=1.E5
 
  ** Thermal expansion coeff (linear)
  ** for constant alpha set EXCOLI=expansion coefficient
  ** if client-specified temperature dependent alpha is
     required, set EXCOLI=-1.0
  ** if linear variation of alpha with T is required, i.e.
     ALPHA = EXCOC1 + EXCOC2*T then set EXCOLI=0, EXCOC1 & EXCOC2
  EXCOLI=1.0E-05
  ** poisson ratio
  POISSN=0.3
  POISSN=0.
  
WIN=1.0
TIN=300.0
tsurr=300
 
  ** grid settings
IYNORT=unit*2
IZLOW=unit+1;IZHI=unit*3
  CARTES=F
 
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
XULAST=0.01
GRDPWR(X,1,XULAST,1.0)
 
    GROUP 4. Y-direction grid specification
NREGY=2;YVLAST=1.0
IREGY=1;GRDPWR(Y,unit*2,0.3,1.0)
IREGY=2;GRDPWR(Y,unit*2,0.3,1.0)
 
    GROUP 5. Z-direction grid specification
NREGZ=3;ZWLAST=1.0
IREGZ=1;GRDPWR(Z,unit,0.1,1.0)  ! diminish z-direction sizes 10-fold
IREGZ=2;GRDPWR(Z,unit*2,0.01,1.0)
IREGZ=3;GRDPWR(Z,unit,0.1,1.0)
 
    GROUP 7. Variables stored, solved & named
   *  Solve for P1, V1, W1 and TEM1 by whole-field method
SOLVE(P1,V1,W1,TEM1)
SOLUTN(P1,Y,Y,Y,N,N,N)
SOLUTN(V1,Y,Y,Y,P,P,P)
SOLUTN(W1,Y,Y,Y,P,P,P)
SOLUTN(TEM1,Y,Y,Y,N,N,Y)
   *  Store other variables
STORE(PRPS,DILA,DVO1,DRH1)
STORE(EPSY,STRY,EPSZ,STRZ,EPST)
STRA=T    ! activate calculation of stress and strain in solid
   *
    GROUP 8. Terms (in differential equations) & devices
TERMS(TEM1,N,Y,Y,Y,Y,Y)
CONVAC=T  ! use the vorticity method as convergence accelerator
 
    GROUP 9.  Properties of the medium (or media)
  ** set via prps values
  TEXT(Choose Fluid Materials
    71 start of ....fluidmat
store(prps)
integer(air20 , airisent,   airideal,   water20,    mercury, freon)
integer(3gasideal, stm100, stmisent, stmideal)
air20 = 0;      airisent=1; airideal=2; 3gasideal=30; stm100=23
stmisent=24;  stmideal=25
water20=67 ;mercury=66; freon=64
    71 end of ....fluidmat
 ** LOAD( 71) from the PHOENICS Input Library
 ** LOAD( 71) from the PHOENICS Input Library
  TEXT(Choose Solid Materials
    70 start of ....solidmat
store(prps)
  The following settings correspond to the IMAT (ie PRPS) values.
  Note that only the first 6 characters of the names of the
  integers are significant
integer(alumin,copper,epoxy,fibregl,steel,glass,phase1,phase2)
alumin= 100;    copper=103; epoxy=104; fibregl=105; steel=111
glass=  106
    70 end of ...solidmat
 ** LOAD( 70) from the PHOENICS Input Library
 ** LOAD( 70) from the PHOENICS Input Library
    GROUP 11. Initialization of fields of variables,
              porosities, etc.
  ** working fluid is air
FIINIT(PRPS)=air20
  ** Initialize Temperature and density (to air density) Field
FIINIT(TEM1)=TIN
  ** Body properties are those of steel
PATCH(BODY,INIVAL,1,NX,1,IYNORT,IZLOW,IZHI,1,1)
INIT(BODY,PRPS,0.0,steel);INIT(BODY,TEM1,0.0,TIN)
 
 
    GROUP 13. Boundary conditions and special sources
PATCH(INLET,LOW,1,NX,1,IYNORT,1,1,1,LSTEP)
COVAL(INLET,P1,FIXFLU,1.189*WIN)
  COVAL(INLET,W1,ONLYMS,WIN)
COVAL(INLET,TEM1,ONLYMS,TIN)
 
  ** outlet boundary condition, name EXIT (at NORTH or HIGH)
PATCH(EXIT,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP)
COVAL(EXIT,P1,1.0,0.0);COVAL(EXIT,TEM1,ONLYMS,SAME)
  ** porous-medium resistances in parts of domain accessible to
     fluid
  PATCH(PORMED1,PHASEM,1,1,1,NY-1,1,IZLOW-1,1,LSTEP)
  COVAL(PORMED1,V1,RESCO,0.0)
  PATCH(PORMED15,PHASEM,1,1,IYNORT+1,NY-1,IZLOW,IZHI,1,LSTEP)
  COVAL(PORMED15,V1,RESCO,0.0)
  PATCH(PORMED2,PHASEM,1,1,1,NY-1,IZHI+1,NZ,1,LSTEP)
  COVAL(PORMED2,V1,RESCO,0.0)
  PATCH(PORMED3,PHASEM,1,1,1,NY,1,IZLOW-2,1,LSTEP)
  COVAL(PORMED3,W1,RESCO,0.0)
  PATCH(PORMED35,PHASEM,1,1,IYNORT+1,NY,IZLOW-1,NZ,1,LSTEP)
  COVAL(PORMED35,W1,RESCO,0.0)
  PATCH(PORMED4,PHASEM,1,1,1,IYNORT,IZHI+1,NZ,1,LSTEP)
  COVAL(PORMED4,W1,RESCO,0.0)
  ** HEAT-SOURCE boundary condition, name HOT
PATCH(HOT,VOLUME,1,NX,1,1,IZLOW,IZHI,1,LSTEP)
COVAL(HOT,TEM1,FIXFLU,heatsor)
 
  ** fix displacement to zero at iy=IYNORT, along larger-z half
  PATCH(FIXV1,NORTH,1,NX,IYNORT,IYNORT,(IZHI+IZLOW)/2+1,IZHI,1,LSTE$
  COVAL(FIXV1,V1,FIXVAL,0.0)
 
  ** hold w1 to zero at base of beam
PATCH(FIXW1,HIGH,1,NX,1,1,IZLOW-1,IZHI,1,LSTEP)
COVAL(FIXW1,W1,FIXVAL,0.0)
  ** hold v1 to zero at south boundary by wall patch
 
PATCH(FIXV2,SWALL,1,NX,1,1,IZLOW,IZHI,1,LSTEP)
COVAL(FIXV2,V1,1.0,0.0)


   
  ** bending the beam
 
  
  
PATCH(BEAM,VOLUME,1,NX,1,IYNORT,IZLOW,IZHI,1,1)
PATCH(TIP,VOLUME,1,1,IYNORT,IYNORT,IZLOW,IZHI,1,1)
  ** provide torque at beam end in form of vorticity gradients
REAL(TORQUE)
TORQUE=-1.0
(source of w1 at beam is :torque: with fixflux)
 
(source of w1 at tip is -2*:unit:*(:torque:) with fixflux)
  ** end of sources specification



LSWEEP=100
 
 
  ** GROUP 16. Termination criteria for inner iterations.
LITER(P1)=20; LITER(V1)=20; LITER(W1)=20; LITER(TEM1)=20
RESREF(P1)=1.E-20;RESREF(V1)=1.E-20
RESREF(W1)=1.E-20;RESREF(W1)=1.E-20
ENDIT(P1)=1.E-20;ENDIT(V1)=1.E-20;
ENDIT(W1)=1.E-20 ;ENDIT(W1)=1.E-20
SELREF=F
  ** GROUP 19. Special data
  SPEDAT(SET,STRAIN,POISSN,R,0.3) !  set Poisson's ratio
relax(w1,linrlx,0.25)
relax(w1,linrlx,0.25)
  ** GROUP 21. Frequency and extent of field printout.
IYPRL=IYNORT
if(unit.eq.5) then
 IZPRF=6 ;IZPRL=15
endif
if(unit.eq.1) then
 IZPRF=1 ;IZPRL=nz
 iyprl=ny
endif
izprf=izlow
izprl=izhi
NPRINT=LSWEEP ; NZPRIN=1 ; NYPRIN=1
 
    GROUP 20. Preliminary print-out
  ** Assign cell-indices of spot-point monitoring location
IXMON=1;IYMON=IYNORT/2;IZMON=(IZLOW+IZHI)/2
    GROUP 23. Variable-by-variable field printout and plot
              and/or tabulation of spot-values and residuals.
  ** GROUP 24. Preparation for continuation runs.
TSTSWP=-1
 
(stored var v1an is anco(v1))
(stored var v1as is asco(v1))
(stored var v1ah is ahco(v1))
(stored var v1al is alco(v1))
(stored var v1ap is apco(v1))
(stored var v1rs is resi(v1))
 
(stored var w1an is anco(w1))
(stored var w1as is asco(w1))
(stored var w1ah is ahco(w1))
(stored var w1al is alco(w1))
(stored var w1ap is apco(w1))
(stored var w1rs is resi(w1))
store(pdcy,pdcz)
  lsweep=5
  nprint=1
dbsoda=t
debug=t
  dbgphi(v1)=t;dbcomp=t;dbcmph=t;dbcmpn=t
iswdb1=1;iswdb2=lsweep
izdb1=3;izdb2=3;dbindx=t
fiinit(w1)=0.0
fiinit(v1)=0.0
dbsol2=t
isolz=0;isoly=0
  izprf=12
  lsweep=1
STOP