2D_Cyl_Piston

TITLE "Piston"

COORDINATES

Ycylinder

SELECT

regrid=off  { disable the adaptive mesh refinement }

painted

nrupdate=5  { increase the number of Newton iterations per timestep, to

   be sure the solution is converged before proceding. }

DEFINITIONS

my_ngrid = 30    { later applied to the NGRID control }

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

smooth = cspeed0*(rad/my_ngrid)

evisc = max(visc,smooth)

SELECT

ngrid= my_ngrid

INITIAL VALUES

rho = rho0

u = 0

v = 0

P = P0

EQUATIONS

{ Eulerian gas equations: (FlexPDE will add motion terms) }

rho:   dt(rho) + dr(rho*u*r)/r + dz(rho*v) = smooth*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)) = smooth*div(grad(P))

vm:  div(grad(vm))=0     { balance mesh velocities }

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