﻿ Sample Problems > Usage > Moving_Mesh > 2d_lagrangian_shock

# 2d_lagrangian_shock

Navigation:  Sample Problems > Usage > Moving_Mesh >

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

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