piston.pde
{ 
  This problem models the flow of a perfect gas in a compressor cylinder. 
  The initial gas pressure is chosen as 1e-4 Atm, to show interesting swirling. 
  The boundaries of the domain are moved according to the oscillation of the piston, 
  while the interior mesh is tensioned within the moving boundaries. 
  This results in a mixed Lagrange/Eulerian model, in that the mesh is moving, 
  but with different velocity than the fluid. 
}  
TITLE "Piston"  
   
COORDINATES  
  Ycylinder  
   
SELECT  
 regrid=off  { disable the adaptive mesh refinement }  
 painted     { paint all contours }
   
DEFINITIONS  
  my_ngrid = 30    { later applied to the NGRID control and smoother }  
  stroke = 8       { cylinder stroke cm }  
  rad = 4            { cylinder bore radius cm }  
  zraise = 1           { raised height of piston center }  
  rraise = 3*rad/4     { radius or raised piston center }  
  gap = 2              { piston/head clearance }  
  gamma = 1.4  
   
  rho0 = 0.001  
  P0 = 100  { initial pressure (dyne/cm2) = 1e-4 Atm }  
  visc = 0.15         { kinematic viscosity, cm^2/sec }  
   
  rpm = 1000  { compressor speed }  
  period = 60/rpm  {seconds}  
  vpeak = (pi*stroke/period)  
  { velocity profile: }  
  vprofile =vpeak*sin(2*pi*t/period)  
  { the piston shape: }  
  zpiston = if r<rraise then zraise else zraise*(rad-r)/(rad-rraise)  
  { the time-dependent piston profile: }  
  zprofile = zpiston+0.5*stroke*(1-cos(2*pi*t/period))  
  ztop = stroke+gap+zraise  { maximum z postion }  
   
VARIABLES  
  rho(rho0/10)      { gas density }  
  u(vpeak/10)       { radial velocity }  
  v(vpeak/10)       { axial velocity }  
  P(P0/10)          { pressure }  
  vm(vpeak/10)      { mesh velocity }  
  zm=move(z)        { mesh position }  
   
DEFINITIONS  
  { sound speed }  
  cspeed = sqrt(gamma*P/rho)  
  cspeed0 = sqrt(gamma*P0/rho0)  
  { a smoothing coefficient: }  
  smoother = cspeed0*(rad/my_ngrid)  
  evisc = max(visc,smoother)  
   
SELECT  
  ngrid= my_ngrid  
   
INITIAL VALUES  
  rho = rho0  
  u = 0  
  v = 0  
  P = P0  
   
EULERIAN EQUATIONS  
  { Eulerian gas equations: (FlexPDE will add motion terms) }  
  rho:   dt(rho) + dr(rho*u*r)/r + dz(rho*v) = smoother*div(grad(rho))  
  u:   dt(u) + u*dr(u)+v*dz(u) + dr(P)/rho  = evisc*div(grad(u))-evisc*u/r^2  
  v:   dt(v) + u*dr(v)+v*dz(v) + dz(P)/rho  = evisc*div(grad(v))  
  P:   dt(P) + u*dr(P)+v*dz(P) + gamma*P*(dr(r*u)/r+dz(v)) = smoother*div(grad(P))  
  vm:  dzz(vm)=0     { balance mesh velocities in z only }  
  zm:  dt(zm)=vm          { node positions - move only in z }  
   
BOUNDARIES  
  { use a piston and compression chamber with beveled edge, to create a swirl }  
  REGION 1  
    START(0,zraise)  
    value(u)=0 value(v)=vprofile value(vm)=vprofile dt(zm)=vprofile  
          line to (rraise,zraise) to (rad,0)  
    value(u)=0 nobc(v) nobc(vm) nobc(zm)  
          line to (rad,stroke+gap)  
    value(u)=0 value(v)=0 value(vm)=0 dt(zm)=0  
          line to (rraise,ztop) to (0,ztop)  
    value(u)=0 nobc(v) nobc(vm) nobc(zm)  
          line to close  
   
    { add a diagonal feature to help control cell shapes at upper corner }  
    FEATURE start(rraise,zraise) line to (rad,stroke+gap)  
   
TIME 0 TO 2*period by 1e-6  
   
PLOTS  
  for t=0 by period/120 to endtime  
  { control the frame size and data scaling to create a useable movie 
    ( the movie can be created by replaying the .PG5 file and selecting 
      EXPORT MOVIE, or we could add PNG() commands here to create 
      it directly) }  
  grid(r,z)  frame(0,0, rad,ztop)  
  contour(rho)  frame(0,0, rad,ztop) fixed range(0,0.01) contours=50  nominmax  
  contour(u)  frame(0,0, rad,ztop) fixed range(-500,500) contours=50  
  contour(v)  frame(0,0, rad,ztop) fixed range(-550,550) contours=50  
  contour(P)  frame(0,0, rad,ztop) fixed range(0,2300) contours=50  nominmax  
  vector(u,v)  frame(0,0, rad,ztop) fixed range(0,550)  
  contour(cspeed) frame(0,0, rad,ztop) fixed range(0,600)  
  contour(magnitude(u,v)/cspeed)  frame(0,0, rad,ztop) fixed range(0,1.1)  
   
  history(vprofile/vpeak,zprofile/stroke)   range(-1,1)  
          report(vpeak)  report(stroke)  
  history(globalmax(P), globalmin(P))  
  history(integral(P))  
  history(globalmax(rho), globalmin(rho))  
  history(integral(rho))  
  history(deltat)  
END