Arbitrary Lagrange/Eulerian Moving Mesh

FlexPDE allows the definition of mesh-moving equations, and applies an Arbitrary Lagrange/Eulerian (ALE) formulation. This allows several optional behaviors:

By locking the mesh to the fluid velocity, you can create a fully Lagrangian model.
Or, you can define a relaxive mesh within moving boundaries to maintain mesh integrity.
Or, by omitting the mesh moving equations, you can perform a fully Eulerian computation.

In any case, FlexPDE automatically corrects the PDE system to represent the motion of the mesh.

The problem shown here computes the motion of a gas in a compressor cylinder.

Animations

Mesh

Gas Pressure

(All animations are FlexPDE native graphical output, assembled into a movie by APNG Maker)

VIEW PROBLEM SCRIPT

{ 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(threshold = rho0/10) { gas density }
u(threshold = vpeak/10) { radial velocity }
v(threshold = vpeak/10) { axial velocity }
P(threshold = P0/10) { pressure }
vm(threshold = 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,610)
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