2d_piston_movingmesh

<< Click to Display Table of Contents >>

Navigation:  Sample Problems > Applications > Fluids >

2d_piston_movingmesh

Previous pageReturn to chapter overviewNext page

 

{  2D_PISTON_MOVINGMESH.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 .PG7 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,2600) 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.5)

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