Hi Daniel,
I updated your code a bit and it seems to do better. See
https://pastebin.com/51xfqWH3 for the full version.
I changed the "eqX" equation to
param = k_c * (P - Peq)
maxX = 0.9
largeValue = 1e9
constraint = largeValue * (X > maxX)
big_value = 1e9
eqX = TransientTerm(var = X) == param - ImplicitSourceTerm(param,
var=X) + maxX * constraint - ImplicitSourceTerm(constraint, var=X)
to have a constraint on the maximum value of X. I'm not sure whether
this is physical, but it does constrain the max value of X.
I changed the "eqY" equation to be,
eqT = TransientTerm(rho_b_MH * Cp_b, var=T) == DiffusionTerm(coeff
= L_eff, var=T) + rho_b_MH*(dH/M)*rc
to have the transient coefficient inside the term. Again, not sure if
that matters, but did it anyway.
I changed the constraint to be
constraint_value = FaceVariable(mesh=mesh)
T.faceGrad.constrain(constraint_value,mesh.facesLeft)
and then updated the constraint in the loop with
for sweep in range(sweeps):
constraint_value[:] = h * ((60. + 273.15) - T.faceValue)
If FiPy worked correctly that wouldn't be necessary, but evidently it doesn't.
Anyway that all seems to work at least as far as I can tell.
Cheers,
Daniel
--
Daniel Wheeler
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]