SCPSpheresElastic

This test simulates the deflagration of a solid composite propellant (SCP) consisting of packed AP in an HTPB matrix. It combines the phase field method for interface tracking with thermal transport and finite-deformation elasticity.

The following full-length simulation results from running the input file directly with no modification. On eight processors, it should take about one hour to run.

Notes

The domain is 0.0005 x 0.0005; however, the sphere packing file is defined over a larger domain (0.001 x 0.001).
The parameters are set to be as course as possible in both space and time without causing instability or solver divergence. This means that the results may or may not be completely accurate. It also may mean that you see some jitter or spurious fluctuation in the elastic solution (especially in the stress and strain fields.) This should go away with increased resolution, so it is recommended that you increase the AMR level and decrease the timestep in order to get properly converged results.
The parameter elastic.solver.fixed_iter replaces elastic.solver.max_iter. This allows the solution to continue even if full convergence is not attained. However, as long as it is inside the Newton solver, the NR convergence check will prevent the solver from continuing unless it is within tolerance.

References

Maycon Meier and Brandon Runnels. Finite kinematics diffuse interface mechanics coupled to solid composite propellant deflagration. Computer Methods in Applied Mechanics and Engineering, 427:117040, 2024.
Maycon Meier, Emma Schmidt, Patrick Martinez, J Matt Quinlan, and Brandon Runnels. Diffuse interface method for solid composite propellant ignition and regression. Combustion and Flame, 259:113120, 2024.

2d-serial-short

	Two-dimensional
	Serial
	Validated using check script
	./bin/alamo-2d-g++ tests/SCPSpheresElastic/input stop_time="0.001"

2d-parallel-long

	Two-dimensional
	Parallel (8 procs)
	Validated using check script
	mpiexec -np 8 ./bin/alamo-2d-g++ tests/SCPSpheresElastic/input stop_time="0.02"

Input file (../../tests/SCPSpheresElastic/input)

#@
#@ [2d-serial-short]
#@ dim=2
#@ nprocs=1
#@ args = stop_time=0.001
#@ check-file = reference/2d-serial-short.dat
#@
#@ [2d-parallel-long]
#@ dim=2
#@ nprocs=8
#@ args = stop_time=0.02
#@ check-file = reference/2d-parallel-long.dat
#@

alamo.program = flame
plot_file = tests/SCPSpheresElastic/output

# AMR Parameters
amr.plot_int = 100
amr.max_level = 5
amr.blocking_factor = 4
amr.base_regrid_int = 100
amr.grid_eff = 1
amr.refinement_criterion = 0.1
amr.refinement_criterion_temp = 10.0
amr.phi_refinement_criterion = 0.5
amr.cell.all = 1
amr.node.all = 1


# Geometry and AMR level specification
amr.n_cell = 8 8 8
geometry.prob_lo = 0.0 0.0 0.0
geometry.prob_hi = 0.0005 0.0005 0.0005
amr.max_level = 4
geometry.is_periodic = 0 0 0

# Timestep and duration
timestep = 0.000001
stop_time = 0.11

# Phase field parameters
pf.eps = 0.000001 
pf.lambda = 0.001
pf.eta.ic.type = constant
pf.eta.ic.constant.value=1.0
pf.kappa = 1.0
pf.w1 = 1.0
pf.w12 = 2.0
pf.w0 = 0.0
small = 1E-4

# Eta boundary condition
pf.eta.bc.type = constant
pf.eta.bc.constant.type.xlo = dirichlet
pf.eta.bc.constant.type.xhi = dirichlet
pf.eta.bc.constant.type.ylo = neumann
pf.eta.bc.constant.type.yhi = neumann
pf.eta.bc.constant.type.zlo = neumann
pf.eta.bc.constant.type.zhi = neumann
pf.eta.bc.constant.val.xlo = 0.0
pf.eta.bc.constant.val.xhi = 1.0
pf.eta.bc.constant.val.ylo = 0.0
pf.eta.bc.constant.val.yhi = 0.0
pf.eta.bc.constant.val.zlo = 0.0
pf.eta.bc.constant.val.zhi = 0.0

# Phi initial condition (spheres)
phi.ic.type = psread
phi.ic.psread.filename = tests/SCPSpheresElastic/2d.xyzr
phi.ic.psread.eps = 1.13e-5 

# Use the full-feedback propellant model
propellant.type = fullfeedback
# Arrhenius rate law params
propellant.fullfeedback.m_ap = 1.45e5
propellant.fullfeedback.m_htpb = 1.04e2
propellant.fullfeedback.E_ap = 11000.0
propellant.fullfeedback.E_htpb = 7500.0
propellant.fullfeedback.mob_ap = 1
propellant.fullfeedback.bound = 300.0
# Thermal transport params
propellant.fullfeedback.rho_ap = 1950.0
propellant.fullfeedback.rho_htpb = 920.0
propellant.fullfeedback.k_ap = 0.4186e0
propellant.fullfeedback.k_htpb = 0.13
propellant.fullfeedback.cp_ap = 1297.9
propellant.fullfeedback.cp_htpb = 2418.29
# Mass-to-heat flux estimation params
propellant.fullfeedback.a1 = 1.114
propellant.fullfeedback.a2 = 0.46
propellant.fullfeedback.a3 = 2.797
propellant.fullfeedback.b1 = 0.323
propellant.fullfeedback.b2 = 0.42
propellant.fullfeedback.b3 = 0.3225
propellant.fullfeedback.c1 = -0.09906
propellant.fullfeedback.mlocal_ap = 1000.0
propellant.fullfeedback.mlocal_comb = 0.0
propellant.fullfeedback.phi.zeta_0 = 1.13e-5
propellant.fullfeedback.phi.zeta = 1.13e-5 


# Activate thermal transport
thermal.on = 1

# Temperature boundary conditions
thermal.temp.bc.type = constant
thermal.temp.bc.constant.type.xlo = neumann
thermal.temp.bc.constant.type.xhi = neumann
thermal.temp.bc.constant.type.ylo = neumann
thermal.temp.bc.constant.type.yhi = neumann
thermal.temp.bc.constant.type.zlo = neumann
thermal.temp.bc.constant.type.zhi = neumann
thermal.temp.bc.constant.val.xlo = 0.0
thermal.temp.bc.constant.val.xhi = 0.0
thermal.temp.bc.constant.val.ylo = 0.0
thermal.temp.bc.constant.val.yhi = 0.0
thermal.temp.bc.constant.val.zlo = 0.0
thermal.temp.bc.constant.val.zhi = 0.0

# Temperature initial conditions
temp.ic.type="constant"
temp.ic.constant.value=500

# Laser source term
laser.ic.type = constant
laser.ic.constant.value = 1.0e6
thermal.hc = 1.0e7

# Pressure
chamber.pressure = 4.0 # [MPa]


# Thermoelasticity
elastic.on = 1
elastic.interval = 100
# Elastic model params
model_ap.mu = 140
model_ap.kappa = 150
model_ap.eps0 = 2.217e-5 0.0 0.0 0.0 2.217e-5 0.0 0.0 0.0 2.217e-5
model_htpb.mu = 8
model_htpb.kappa = 210
model_htpb.eps0 = 51e-7 0.0 0.0 0.0 51e-7 0.0 0.0 0.0 51e-7

# Elastic solver information
elastic.solver.verbose = 4
elastic.solver.fixed_iter = 500
elastic.solver.nriters = 200
elastic.solver.nrtolerance = 1E-5
elastic.print_model = 1
elastic.print_residual = 1

# Mechanics boundary conditions
elastic.bc.type = constant
elastic.bc.constant.type.xlo = trac trac 
elastic.bc.constant.type.ylo = trac disp 
elastic.bc.constant.type.yhi = trac disp 
elastic.bc.constant.type.xloylo = trac disp 
elastic.bc.constant.type.xloyhi = trac disp 

# No traction on free surface
elastic.traction = 0.0