This setup illustrates a small example on the setup of a moving geometry. The fluid is initially in rest. A piston (geometry) is located in the middle area of the domain that has a length of 1.0. The piston moves with Mach 0.4 towards the right boundaty, due to the sudden movement, a shock is formed ahead of the piston. Behind the piston rarefaction can be observed. The piston is modelled as a porous medium. An exact solution can be found in literature for this test case e.g. in Toro. This test case is used to ensure conservation of mass, momentum and energy. If a shock is not formed ahead of the piston or it has a wronge velocity, conservation is not maintained.
The complete configuration is provided in ateles.lua:
logging = {level = 10}
--debug = {logging = {level=10}}
-------------------------------------------
order = 16
degree = order - 1
nElems = 128
------------------------------------------
tmax= 5 * 1e-06
-- global simulation options
simulation_name='ateles'
--wallclock = 2*60*60 - 4*60
sim_control = {
time_control = {
min = 0,
max = tmax,
--max = {sim=tmax, clock = wallclock},
interval = {iter = 1},
}
}
check = {interval = 1}
---- Restart settings
NOrestart = {
-- read = './restart/ateles_lastHeader.lua',
write = './restart/',
time_control = {
min = 0,
max = tmax,
interval = tmax/20,
align_trigger = {sim = true},
},
}
-- Segments for tracking--
segments = math.ceil(1.0/(1/nElems)) * (degree +1) * 4
-- tracking --
tracking = {
{ label = 'line',
variable = {'density', 'pressure', 'velocity'},
shape = {kind='canoND',
object = { origin = {0.0, 0.0, 0.0},
vec = {1.0, 0.0, 0.0},
segments = {segments, segments, segments}},
},
time_control = {
min = 0,
max = tmax,
interval = tmax
},
folder = './',
output = {format = 'asciiSpatial', use_get_point = true}
},
}
-- Mesh definitions --
mesh = {
predefined = 'line_bounded',
origin = {0.0},
length = 1.0,
element_count = nElems
}
--physical data
gamma = 1.4
velocityX = 150.0
dens = 1.0
press = 1.0e5
-- Define a piston inside the domain --
xmin = 0.4
xmax = 0.44
temp = press/(dens*(1.0/1.4))
vel_init = 0.0
c = math.sqrt((press*gamma)/dens)
u = velocityX
a = - (gamma/(2*c^2))*((4*c^2/gamma)+(gamma + 1) * u^2)
b = (gamma/(2*c^2))* ((2*c^2/gamma) - u^2*(gamma - 1))
p_ratio1 = -a/2 + math.sqrt((a/2)^2 - b)
p_ratio2 = -a/2 - math.sqrt((a/2)^2 - b)
density_ratio = (1 + ((gamma + 1)/(gamma - 1 ))*
math.max(p_ratio1,p_ratio2))/(((gamma + 1)/(gamma -1))+
math.max(p_ratio1,p_ratio2))
densL = density_ratio * dens
pressL = math.max(p_ratio1,p_ratio2) * press
temperature_ratio = (1 +(pressL - press)/press ) *
((2*gamma + (gamma-1)*((pressL-press)/press))/
(2*gamma + (gamma+1)*((pressL - press)/press)))
tempL = temperature_ratio * temp
function inside_piston(x,y,z,t)
xa = xmin + velocityX * t
xb = xmax + velocityX * t
if (xa <= x and x <= xb ) then
return true
else
return false
end
end
function velocityRelax(x,y,z,t)
return velocityX
end
function temperature(x,y,z,t)
xa = xmin + velocityX * t
xb = xmax + velocityX * t
diff_half = (xb -xa)/2
if (x >= xa and x <=(xa + diff_half)) then
return temp
elseif (x >(xa + diff_half) and x<=xb) then
return tempL
else
return temp
end
end
function characteristic(x,y,z,t)
if inside_piston(x,y,z,t) then
return 1.0
else
return 0.0
end
end
dx = 0.00000001
variable = {
{
name = 'Xi',
ncomponents = 1,
vartype = "st_fun",
st_fun = {
{ const = 0.0 },
{
fun = characteristic,
shape = {
kind = 'canoND',
object = {
origin = {0.3 ,0.0, 0.0},
vec = {
{ 0.6, 0.0, 0.0 },
{ 0.0, dx, 0.0 },
{ 0.0, 0.0, dx },
},
}
}
}
}
},
{
name = 'relax_velocity',
ncomponents = 1,
vartype = "st_fun",
st_fun = {
{ const = 0.0 },
{
fun = velocityRelax,
shape = {
kind = 'canoND',
object = {
origin = {0.3 ,0.0,0.0},
vec = {
{ 0.6, 0.0, 0.0 },
{ 0.0, dx, 0.0 },
{ 0.0, 0.0, dx },
},
}
}
}
}
},
{
name = 'relax_temperatur',
ncomponents = 1,
vartype = "st_fun",
st_fun = {
{ const = 0.0 },
{
fun = temperature,
shape = {
kind = 'canoND',
object = {
origin = {0.3 ,0.0,0.0},
vec = {
{ 0.6, 0.0, 0.0 },
{ 0.0, dx, 0.0 },
{ 0.0, 0.0, dx },
},
}
}
}
}
},
}
-- timing settings (i.e. output for performance measurements, this table is otional)
timing_file = 'timing.res' -- the filename of the timing results
phi = 1.0
beta = 1e-6
eta_v = phi^2 * beta^2
eta_t = 0.4 * phi * beta
-- Equation definitions --
equation = {
name = 'euler_1d',
numflux = 'hll',
isen_coef = 1.4,
r = 1.0/1.4,
porosity = phi,
viscous_permeability = eta_v,
thermal_permeability = eta_t,
material = {
characteristic = 'Xi',
relax_velocity = 'relax_velocity',
relax_temperature = 'relax_temperatur',
modereduction = false,
}
}
-- (cv) heat capacity and (r) ideal gas constant
equation["cv"] = equation["r"] / (equation["isen_coef"] - 1.0)
-- Scheme definitions --
scheme = {
-- the spatial discretization scheme
spatial = {
name = 'modg_1d',
modg_space = 'Q',
m = degree,
},
stabilization = {
{
name = 'spectral_viscosity',
alpha = 36,
order = 24,
isAdaptive = true,
--recovery_order = 1.0e-4,
recovery_density = 1.0e-1,
recovery_pressure = 1.0e-1
},
--{
-- name = 'covolume',
-- alpha = 36,
-- order = 24,
-- beta = 1.0 - 2.0/degree,
-- },
},
-- the temporal discretization scheme
temporal = {
name = 'imexRungeKutta',
steps = 4,
control = {
name = 'cfl',
cfl = 0.2,
},
},
}
-- ...the general projection table
projection = {
kind = 'l2p',
material = {
factor = 3.0
},
}
vel_init = 0.0
c = math.sqrt((press*gamma)/dens)
u = velocityX
a = - (gamma/(2*c^2))*((4*c^2/gamma)+(gamma + 1) * u^2)
b = (gamma/(2*c^2))* ((2*c^2/gamma) - u^2*(gamma - 1))
p_ratio1 = -a/2 + math.sqrt((a/2)^2 - b)
p_ratio2 = -a/2 - math.sqrt((a/2)^2 - b)
density_ratio = (1 + ((gamma + 1)/(gamma - 1 ))* math.max(p_ratio1,p_ratio2))/
(((gamma + 1)/(gamma -1))+ math.max(p_ratio1,p_ratio2))
densL = density_ratio * dens
pressL = math.max(p_ratio1,p_ratio2) * press
ushock= (u/c)/(1 - density_ratio^(-1))* c
velocity_2 = ushock - u
velocity_1 = ushock
x_position = 0.44 + ushock*tmax
-- Initial condition
function iniVel(x,y,z,t)
if inside_piston(x,y,z,0.0) then
return velocityX
else
return 0.0
end
end
function inipress(x,y,z,t)
if (x > 0.42 and x < 0.44) then
return pressL
else
return press
end
end
function inidens(x,y,z,t)
if (x > 0.42 and x < 0.44) then
return densL
else
return dens
end
end
initial_condition = {
density = inidens,
pressure = inipress,
velocity = iniVel,
}
-- Boundary definitions
boundary_condition = {
{
label = 'west',
kind = 'outflow',
pressure = press,
}
,
{
label = 'east',
kind = 'outflow',
pressure = press,
}
,
}
Projection: l2p
Polynomial representation: Q
Filtering: spectral_viscosity
Timestepping: imexRungeKutta, 4 steps
Boundary conditions: outflow
Others: porous material (geometry), over-integration for the geometry (Piston) lua function for geometry definition