  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 }
zraise = 1           { raised height of 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: }
{ 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 }
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: }
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
value(u)=0 nobc(v) nobc(vm) nobc(zm)
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 }

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) }
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