|
<< Click to Display Table of Contents >> 2d_lagrangian_shock |
  
|
|
<< Click to Display Table of Contents >> 2d_lagrangian_shock |
|
{ 2D_LAGRANGIAN_SHOCK.PDE
This example demonstrates moving meshes in 2D by solving Sod's shock tube problem
(a 1D problem) on a 2D moving mesh.
Mesh nodes are given the local fluid velocity, so the model is fully Lagrangian.
Ref: G.A. Sod, "A Survey of Several Finite Difference Methods for Systems of Nonlinear
Hyperbolic Conservation Laws", J. Comp. Phys. 27, 1-31 (1978)
See also Kershaw, Prasad and Shaw, "3D Unstructured ALE Hydrodynamics with the
Upwind Discontinuous Finite Element Method", UCRL-JC-122104, Sept 1995.
Other versions of this problem can be found in the "Applications | Fluids" folder.
}
TITLE "Sod's Shock Tube Problem - 2D Lagrangian"
SELECT
ngrid= 100
regrid=off
VARIABLES
rho(1)
u(1)
P(1)
xm=move(x)
DEFINITIONS
len = 1
wid = 0.01
gamma = 1.4
{ define a damping term to kill unwanted oscillations }
eps = 0.001
v = 0
rho0 = 1.0 - 0.875*uramp(x-0.49, x-0.51)
p0 = 1.0 - 0.9*uramp(x-0.49, x-0.51)
INITIAL VALUES
rho = rho0
u = 0
P = p0
EULERIAN EQUATIONS
{ equations are stated as appropriate to the Eulerian (lab) frame.
FlexPDE will add motion terms to convert to Lagrangian form for moving mesh }
{ since the equation is really in x only, we add dyy(.) terms with natural(.)=0
on the sidewalls to impose uniformity across the fictitious y coordinate }
rho: dt(rho) + u*dx(rho) + rho*dx(u) = dyy(rho) + eps*dxx(rho)
u: dt(u) + u*dx(u) + dx(P)/rho = dyy(u) + eps*dxx(u)
P: dt(P) + u*dx(P) + gamma*P*dx(u) = dyy(P) + eps*dxx(P)
xm: dt(xm) = u
BOUNDARIES
REGION 1
{ we must impose the same equivalence dt(xm)=u on the side boundaries
as in the body equations: }
START(0,0)
natural(u)=0 dt(xm)=u line to (len,0)
value(xm)=len value(u)=0 line to (len,wid)
dt(xm)=u natural(u)=0 line to (0,wid)
value(xm)=0 value(u)=0 line to close
TIME 0 TO 0.375
MONITORS
for cycle=5
grid(x,10*y)
elevation(rho) from(0,wid/2) to (len,wid/2) range (0,1)
elevation(u) from(0,wid/2) to (len,wid/2) range (0,1)
elevation(P) from(0,wid/2) to (len,wid/2) range (0,1)
PLOTS
for t=0 by 0.02 to 0.143, 0.16 by 0.02 to 0.375
grid(x,10*y)
elevation(rho) from(0,wid/2) to (len,wid/2) range (0,1)
elevation(u) from(0,wid/2) to (len,wid/2) range (0,1)
elevation(P) from(0,wid/2) to (len,wid/2) range (0,1)
END