#cls                                                                            
TEXT(Piston moved by crank & conn. rod                                          
title                                                                           
libref=783                                                                      
                                                                                
  DISPLAY                                                                       
                                    valve vent                                  
              /////////////////////›  ›///                                      
              --------------------------/                                       
                     ------            ›/                                       
         crank/ .      pis›            ›/                                       
             / conn.   ton›            ›/                                       
            /      rod.   ›            ›/                                       
           *_ . _ . _ . _*. _ . _ . _ .›_ . _  . Symmetry axis                  
                          ›            ›/                                       
                          ›            ›/                                       
                                                                                
                     z1  zwn                                                    
                                                                                
   The compression of gas by a piston-cylinder combination is                   
   simulated by the use of the In-Form BOX function.                            
                                                                                
   The "fluid" within the piston is given a velocity which                      
   corresponds the rotation of the crank.                                       
                                                                                
                                                                                
   At a precribed instant a valve provides a fixed-pressurse outlet             
   in the cylinder wall. This is achieved by use of a boundary                  
   condition embodying the In-Form IF condition.                                
  ENDDIS                                                                        
                                                                                
    GROUP 1. Run title and other preliminaries                                  
TEXT(3D Piston by In-Form)                                                      

   INFORM1BEGIN
           Initial data of problem                                              
           -----------------------                                              
REAL(CRAD); CRAD=.05            ! Cylinder radius, m.                           
REAL(CHIG); CHIG=.1111          ! Cylinder height, m.                           
REAL(CR); CR=.050               ! Crank radius, m.                              
REAL(DR); DR=4.                 ! Connecting-rod centre distance                
                                ! /crank radius                                 
REAL(RPS); RPS=40.              ! Revolutions per second.                       
REAL(VOLF); VOLF=0.8E-3*RPS     ! Mean vol flowrate                             
REAL(ANGFRST); ANGFRST=0.       ! Crank angle at start, deg.                    
REAL(ANGLAST); ANGLAST=180.     ! CrANGLASTk engle at end, deg.                 
INTEGER(KNZ1)                   ! Z cell number where outlet vent               
                                ! starts                                        
                                                                                
           Auxiliary data                                                       
           --------------                                                       
REAL(PI,PI2); PI=3.14159; PI2=2.*PI  ! Pythagoras's constant                    
REAL(ANGLE); ANGLE=ANGLAST - ANGFRST ! Crank-angle change, deg.                 
REAL(DEGPS); DEGPS=RPS*360.     ! Revolutions per second, radians.              
REAL(DET)                       ! Time step, sec.                               
REAL(AV); AV=RPS*PI*2.          ! Angular velocity                              
REAL(ANGL)                      ! Angle                                         
REAL(Z1)                        ! Z coordinate of piston at start               
CHAR(FORM)                      ! Formula                                       
   INFORM1END
                                                                                
    GROUP 2. Transience; time-step specification                                
STEADY=F                                                                        
TFIRST=ANGFRST/DEGPS                                                            
TLAST=ANGLAST/DEGPS                                                             
LSTEP=10                                                                        
DET=(TLAST-TFIRST)/LSTEP                                                        
GRDPWR(T,LSTEP,TLAST-TFIRST,1.0)     ! uniform time steps                       
                                                                                
    GROUP 3. X-direction grid specification                                     
CARTES=F;NX=4                                                                   
XFRAC(1)=0.3;XFRAC(2)=0.5;XFRAC(3)=0.7;XFRAC(4)=1.                              
XULAST=PI                                                                       
                                                                                
    GROUP 4. Y-direction grid specification                                     
NY=4                                                                            
YFRAC(1)=0.5; YFRAC(2)=.8; YFRAC(3)=.95; YFRAC(4)=1.                            
YVLAST=CRAD                                                                     
                                                                                
    GROUP 5. Z-direction grid specification                                     
NZ=100                                                                          
GRDPWR(Z,NZ,CHIG,1.0)                                                           
KNZ1=NZ-5                                                                       
    GROUP 7. Variables stored, solved & named                                   
SOLVE(P1,U1,V1,W1,H1);SOLUTN(P1,Y,Y,Y,N,N,N)                                    
STORE(RHO1,MARK)                                                                
                                                                                
    GROUP 9. Properties of the medium (or media)                                
ENUL=1.E-5;ENUT=1.E-3                                                           
   ** Temperature from enthalpy field.                                          
TMP1=LINH                                                                       
   **  the specific heat...                                    
CP1=1000.
TMP1B=1./cp1                                                                   
RHO1=IDEALGAS        ! density from ideal gas law                               
DRH1DP=IDEALGAS      ! compressibility likewise                                 
RHO1B=29./8305.6     ! molecular weight of air/Gas_Constant                     
                                                                                
    GROUP 11. Initialization of variable or porosity fields                     
FIINIT(P1)=1.0E5;FIINIT(H1)=9.0E5                                               
                                                                                
    GROUP 13. Boundary conditions and special sources                           
WALL (HEAD,LOW,1,NX,1,NY,1,1,1,LSTEP)    ! Cylinder Head                        
WALL (CYL,NORTH,1,NX,NY,NY,1,NZ,1,LSTEP) ! Cylinder Wall                        
WALL (PCR,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP)  ! Piston Crown                         
                                                                                
  INFORM13BEGIN                                                                 
           In-Form statements                                                   
           ------------------                                                   
                                                                                
   ** Piston geometry at first time step                                                          
(MAKE of ZWO)                                                                   
(MAKE of ZWN)                                                                   
(MAKE of VENT)                                                                  
                                                                                
   **** ZWO is Z coordinate of piston at previous time step                   
PATCH(ONECELL1,CELL,1,1,1,1,1,1,1,1)                                            
(STORE1 of ZWO at ONECELL1 is 0)                                                
ANGL=ANGFRST*PI/180.+AV*DET                                                     
                                                                                
   ***** Z1=CR*(SQRT(DR**2-SIN(ANGL)**2)-COS(ANGL)-DR+1)                        
Z1=DR**2-SIN(ANGL)**2                                                           
Z1=CR*(SQRT(Z1)-COS(ANGL)-DR+1.)                                                
                                                                                                                                                
   **** ZWN is Z coordinate of piston at current time step                      
(STORE1 of ZWN at ONECELL1 is Z1)                                             
PATCH(PISTON1,CELL,1,NX,1,NY,1,NZ,1,1)                                          
                                                                                                                                               
   **** Using BOX function for formation of piston geometry                       
(INFOB at PISTON1 is BOX(0.,0.,0.,.1,.1,Z1,0.,0.,0.$                          
) with INFOB_1)                                                                 
   
  !!!!! note the use of the IF condition !!!!!                                                                             
   *** Second and next time steps                                               
PATCH(ONECELL2,CELL,1,1,1,1,1,1,2,LSTEP)                                        
(STORE1 of ZWO at ONECELL2 is ZWN with IF(ISWEEP.EQ.1))                         
                                                                                
FORM=:CR:*(SQRT(:DR:^2-SIN(AV*TIM)^2)-COS(AV*TIM)-:DR:+1.)                      
(STORE1 of ZWN at ONECELL2 is :FORM: with$                                      
 IF(ISWEEP.EQ.1))                                                               
                                                                                
   **** Using BOX function for formation of piston geometry                     
PATCH(PISTON2,CELL,1,NX,1,NY,1,NZ,2,LSTEP)                                      
(INFOB at PISTON2 is BOX(0.,0.,0.,.1,.1,ZWN,0.,0.,0.)$                          
 with INFOB_1)                                                                  
                                                                                
   ** Velocities components of piston                                           
                                                                                
   *** First time step                                                          
(SOURCE of W1 at PISTON1 is Z1 with INFOB_1!FIXVAL)                           
                                                                                
   *** Second and next time steps                                               
(SOURCE of W1 at PISTON2 is (ZWN-ZWO)/DT with$                                  
 INFOB_1!FIXVAL)                                                                
                                                                                
   ** Set the switch for opening of the outlet valve                            
      VENT=0 - the output vent is closed                                        
      VENT=1 - the output vent is opened                                        
                                                                                
   *** Condition of opening of the outlet vent                                  
FORM=IF(ISWEEP.EQ.LSWEEP.AND.P1[,,NZ-1].GE.3.E5)                                
   *** First time step                                                          
(STORE1 of VENT at ONECELL1 is 0.)                                              
   *** Second and next time steps                                               
       Set VENT is 1 when FORM is TRUE                                          
(STORE1 of VENT at ONECELL2 is 1. with :FORM:)                                  
                                                                                
   ** Boundary condition at the output vent                                     
                                                                                
PATCH(EXHAUST,NORTH,1,1,NY,NY,KNZ1,NZ,2,LSTEP)                                  
(SOURCE of P1 at EXHAUST is -10. with IF(VENT.NE.0.))                           
                                                                                
   ** MARK variable is used for the image of the piston position                
                                                                                
   *** Initialization of MARK variable                                          
(STORED of MARK is 0 with IF(ISWEEP.EQ.1))                                      
   *** Image of a location of the piston                                        
(STORED of MARK is 1 with INFOB_1!IF(ISWEEP.EQ.1))                              
  INFORM13END                                                                   
                                                                                
    GROUP 15. Termination of sweeps                                             
LSWEEP=30                                                                       
RESREF(P1)=2.0E-4*VOLF                                                          
RESREF(W1)=0.2*RPS*RESREF(P1)                                                   
RESREF(U1)=RESREF(W1);RESREF(V1)=RESREF(W1)                                     
RESREF(H1)=RESREF(P1)*300000.                                                   
                                                                                
    GROUP 16. Termination of iterations                                         
LITER(P1)=10;LITER(U1)=1;LITER(V1)=1;LITER(W1)=1;LITER(H1)=5                    
                                                                                
    GROUP 17. Under-relaxation devices                                          
RELAX(P1,LINRLX,.05)                                                            
RELAX(V1,FALSDT,5.E-5)                                                          
RELAX(W1,FALSDT,5.E-5)                                                          
RELAX(U1,FALSDT,5.E-5)                                                          
                                                                                
    GROUP 18. Limits on variables or increments to them                         
VARMIN(P1)=1.E4                                                                 
                                                                                
    GROUP 22. Spot-value print-out                                              
IXMON=1;IYMON=1;IZMON=NZ-1                                                      
TSTSWP=-1                                                                       
                                                                                
    GROUP 23. Field print-out and plot control                                  
NYPRIN=NY/5;NZPRIN=NZ/5;NPRINT=LSWEEP;NTPRIN=5                                  
XZPR=T;IYPRF=1;IYPRL=1                                                          
IPLTL=LSWEEP                                                                    
                                                                                
  Group 24. Dumps For Restarts                                                  
IDISPA=1; CSG1='P'                                                              
DISTIL=T                                                                        
EX(P1)=1.912E+05                                                                
EX(U1)=1.406E-01; EX(V1)=8.936E-02; EX(W1)=2.305E+00                            
EX(H1)=1.001E+06; EX(MARK)=8.700E-01; EX(RHO1)=4.367E-01                        
libref=783                                                                      
  PHOTON USE                                                                    
  p                                                                             
  p1                                                                            
                                                                                
                                                                                
  up z; gr z 1                                                                  
  gr out x 1; gr out x 4                                                        
  con w1 x 1 fi;0.001;  con w1 x 3 fi;0.001                                     
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Velocity contours at time step=1                                          
  msg Press RETURN to continue                                                  
  pause; vec off; red                                                           
  cont p1 x 1 fil;0.001                                                         
  cont p1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Pressure contours at time step=1                                          
  msg Press RETURN to continue                                                  
  pause                                                                         
  cont h1 x 1 fil;0.001                                                         
  cont h1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Temperature contours at time step=1                                       
  msg Press RETURN to continue                                                  
  pause                                                                         
  p                                                                             
  p5                                                                            
                                                                                
                                                                                
  up z; gr z 1                                                                  
  gr out x 1; gr out x 4                                                        
  con w1 x 1 fi;0.001;  con w1 x 3 fi;0.001                                     
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Velocity contours at time step=5                                          
  msg Press RETURN to continue                                                  
  pause; vec off; red                                                           
  cont p1 x 1 fil;0.001                                                         
  cont p1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Pressure contours at time step=5                                          
  msg Press RETURN to continue                                                  
  pause                                                                         
  cont h1 x 1 fil;0.001                                                         
  cont h1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Temperature contours at time step=5                                       
  msg Press RETURN to continue                                                  
  pause                                                                         
  p                                                                             
  p10                                                                           
                                                                                
                                                                                
  up z; gr z 1                                                                  
  gr out x 1; gr out x 4                                                        
  con w1 x 1 fi;0.001;  con w1 x 3 fi;0.001                                     
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Velocity contours at time step=10                                         
  msg Press RETURN to continue                                                  
  pause; vec off; red                                                           
  cont p1 x 1 fil;0.001                                                         
  cont p1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Pressure contours at time step=10                                         
  msg Press RETURN to continue                                                  
  pause                                                                         
  cont h1 x 1 fil;0.001                                                         
  cont h1 x 3 fil;0.001                                                         
  surf mark x .99 col 4                                                         
  surf mark y .99 col 4                                                         
  msg Temperature contours at time step=10                                      
  msg Press e to END                                                            
  ENDUSE