[johnuthin]

Hello Foamers,

I hope this is the right place to post.

I am trying to solve for a liquid molten-salt nuclear fuel with a custom solver I have. The difficulty arises in modelling delayed neutron precursor drift in the solver which uses an advection term.

I am solving the following using SIMPLE.



When adding this last term (), the solution seems to diverge immediately, even when working on a uniform 2D cartesian mesh. To implement this equation, I have written a series of 8 different delayed neutron equations that mimic the following:

	fvScalarMatrix delayedThree
	(
		fvm::ddt((1.0-keffNew)*rho, concThree)
		- Beta_3/keffNew*rho*fissionTerm
		+ fvm::Sp(lamb_3*rho, concThree)
		+ fvm::div(phi,concThree)
	);

I have also multiplied all terms by rho in the above to account for the unit difference between phi and U in my implementation. Its possible that I have made a very bad and silly mistake and a second set of eyes might be helpful. I am also ignoring heat transfer in my implementation so rho is constant throughout the entire problem.

The error I get once I implement this last term is below. I am still learning openfoam through this implementation and I am asking if there is something obvious that I am missing in the advection term that is causing the FPE. All other parts of the solver work fine so far up to this point. I have also tried preconditioning my velocity fields and that leads to the same error. I am currently using Gauss upwind in my divSchemes and have had no other luck on any divSchemes that I have tested.

Thank you for any help!!!

:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib/x86_64-linux-gnu/libpthread.so.0
#3  Foam::scalarProduct<double, double>::type Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::PBiCGStab::scalarSolve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
#6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#8  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#9  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
#10  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#11  Foam::fvMatrix<double>::solve() at ??:?
#12  ? at ??:?
#13  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#14  ? at ??:?
Floating point exception