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.0000_m 0.0000_m 0.0000_m
geometry.prob_hi = 0.0005_m 0.0005_m 0.0005_m
amr.max_level = 4
geometry.is_periodic = 0 0 0

# Timestep and duration
timestep = 0.000001_s
stop_time = 0.11_s

# Phase field parameters
pf.eps = 0.000001_m
pf.lambda = 0.001_J/m^2
pf.eta.ic.type = constant
pf.eta.ic.constant.value=1.0
pf.kappa = 1.0_J/m^2
pf.w1 = 1.0_1
pf.w12 = 2.0_1
pf.w0 = 0.0_1
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.file.unit = m
phi.ic.psread.file.name = tests/SCPSpheresElastic/2d.xyzr
phi.ic.psread.eps = 1.13e-5_m

# Use the full-feedback propellant model
propellant.type = fullfeedback
# Arrhenius rate law params
propellant.fullfeedback.m_ap = 1.45e5_m/s
propellant.fullfeedback.m_htpb = 1.04e2_m/s
propellant.fullfeedback.E_ap = 11000.0_K
propellant.fullfeedback.E_htpb = 7500.0_K
propellant.fullfeedback.mob_ap = true
propellant.fullfeedback.bound = 300.0_K
# Thermal transport params
propellant.fullfeedback.rho_ap = 1950.0_kg/m^3
propellant.fullfeedback.rho_htpb = 920.0_kg/m^3
propellant.fullfeedback.k_ap = 0.4186e0_W/m/K
propellant.fullfeedback.k_htpb = 0.13_W/m/K
propellant.fullfeedback.cp_ap = 1297.9_J/kg/K
propellant.fullfeedback.cp_htpb = 2418.29_J/kg/K
# Mass-to-heat flux estimation params
propellant.fullfeedback.a1 =  1.11400_1
propellant.fullfeedback.a2 =  0.46000_1
propellant.fullfeedback.a3 =  2.79700_1
propellant.fullfeedback.b1 =  0.32300_1
propellant.fullfeedback.b2 =  0.42000_1
propellant.fullfeedback.b3 =  0.32250_1
propellant.fullfeedback.c1 = -0.09906_1
propellant.fullfeedback.mlocal_ap = 1000.0_kg/m^3/s
propellant.fullfeedback.mlocal_comb = 0.0_kg/m^3/s
propellant.fullfeedback.phi.zeta_0 = 1.13e-5_m
propellant.fullfeedback.phi.zeta = 1.13e-5 _m
propellant.fullfeedback.Pref = 1.0_MPa 


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

# Laser source term
laser.ic.type = constant
laser.ic.constant.value = 1.0e6_W/m^2
thermal.hc = 1.0e7_W/m^2

# Pressure
chamber.pressure = 4.0_MPa # [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