Jon, thanks for catching that typo. I wish it were the problem! The solution with grad phi=0 does not satisfy the zero-flux boundary condition on species 3. In the absence of a potential gradient, that species would need to also have zero gradient (like species 2 in the current solution).
So the question becomes how to enforce that zero-flux BC. Any suggestions? I could drop the conservation equation and replace it with something that enforces zero-flux everywhere, but I’d prefer to use a boundary condition as that would be more applicable in higher dimensions. ☘ > On Jun 12, 2019, at 1:36 PM, Guyer, Jonathan E. Dr. (Fed) via fipy > <[email protected]> wrote: > > Scott - > > Eq2 appears to have a typo: > > -Eq2 = DiffusionTerm(coeff=z2*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) > +Eq2 = DiffusionTerm(coeff=z2*C2, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) > > Eq3 has (much) better coupling if written like this: > > -Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) + (C3.faceGrad).divergence > +Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) - DiffusionTerm(coeff=-z1/z3, > var=C1) - DiffusionTerm(coeff=-z2/z3, var=C2) > > [Note that this way you don't need to evaluate .faceGrad at the boundaries, > so your issue #650 doesn't come up. We don't upwind at the boundaries because > it's complicated and expensive to do right and it's unusual to care.] > > With these changes, the solution converges in one sweep, but does not agree > with the Analytical Solution you give from West. > > I don't understand a couple of things though. As posed, what governs \Phi? I > get \Phi = 0 and see no reason why it would be otherwise. Typically, I would > expect Poisson's equation to be obeyed. Since you assert charge > electroneutrality, this reduces to Laplace's equation, but in that case, the > slope of \Phi should be constant (unless the permittivity is non-uniform, but > you don't have a permittivity). In the Analytical plot, \Phi has definite > curvature. > > My guess is you don't need Eq3 (it'll get taken care of automatically) and > you need some form of Poisson's equation. > > - Jon > >> On Jun 11, 2019, at 4:48 PM, Scott Calabrese Barton <[email protected]> wrote: >> >> I’ve had some trouble using fipy for the 1D, steady-state ternary >> electrolyte problem, which involves coupled, nonlinear migration and >> diffusion. >> >> A text version of the model is listed below. A more complete development is >> at >> https://urldefense.proofpoint.com/v2/url?u=https-3A__gcc01.safelinks.protection.outlook.com_-3Furl-3Dhttp-253A-252F-252Fbit.ly-252F2wJDXob-26amp-3Bdata-3D02-257C01-257Cjonathan.guyer-2540nist.gov-257C1687dc7e080545a63f0708d6eeae58e7-257C2ab5d82fd8fa4797a93e054655c61dec-257C1-257C1-257C636958829930551626-26amp-3Bsdata-3DcSBGKKJG0-252Ft-252FMtTZbnTwenKAVsKPHa-252F1C8NgupkRI0E-253D-26amp-3Breserved-3D0&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=mskW-yZCCtkB6VJ1e2iNv6InRpjYLWWU0Xklumrtr0o&e= >> >> >> from fipy import * >> >> # Parameters >> nx=50 >> C20=4.0/3.0 >> z1, z2, z3= 2.0, -2.0, 1.0 >> >> # Mesh >> mesh = Grid1D(nx=nx, dx=1.0/nx) >> >> # Variables >> C1= CellVariable(mesh = mesh, value = 1.0) >> C2= CellVariable(mesh = mesh, value = 4.0/3.0) >> Phi = CellVariable(mesh = mesh, value = -0.5) >> >> # Electroneutrality >> C3 = -1/z3 * (z1*C1 + z2*C2) >> >> # Boundary Conditions >> C1.constrain(1.0, mesh.facesRight) >> C1.constrain(0.0, mesh.facesLeft) >> C2.constrain(C20, mesh.facesRight) >> Phi.constrain(0, mesh.facesRight) >> >> # Governing Equations >> Eq1 = DiffusionTerm(coeff=z1*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C1) >> Eq2 = DiffusionTerm(coeff=z2*C1, var=Phi) + DiffusionTerm(coeff=1.0, var=C2) >> Eq3 = DiffusionTerm(coeff=z3*C3, var=Phi) + (C3.faceGrad).divergence >> >> Eqns = Eq1 & Eq2 & Eq3 >> >> # Solution >> res = 1e+10 >> restol= 1e-3 >> >> while res > restol: >> res = Eqns.sweep() >> print(res) >> >> >> However, this system does not converge. I suspect the issue is with boundary >> conditions. I’ve tried zeroing out the diffusion coefficients at the >> boundary, but that does not seem to make a difference here. I am able to use >> this approach for the analogous binary problem: >> https://urldefense.proofpoint.com/v2/url?u=https-3A__gcc01.safelinks.protection.outlook.com_-3Furl-3Dhttp-253A-252F-252Fbit.ly-252F2wOzxfK-26amp-3Bdata-3D02-257C01-257Cjonathan.guyer-2540nist.gov-257C1687dc7e080545a63f0708d6eeae58e7-257C2ab5d82fd8fa4797a93e054655c61dec-257C1-257C1-257C636958829930551626-26amp-3Bsdata-3DCkRvDWJQl-252FjHRTwtxd6y1Krs3VuALc2bAOXKJYsdmEk-253D-26amp-3Breserved-3D0&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=IOcn8PY99cvXW8ZO26L9V7s11gHIs2_uDyFOvmq8yKg&e= >> >> >> Any insights or suggestions would be welcome. >> >> Cheers, >> Scott☘ >> >> >> _______________________________________________ >> fipy mailing list >> [email protected] >> https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=dTROrUlTOg_pDStEAloQau0rkmx3ZuaHQy-NSohak0c&e= >> >> [ NIST internal ONLY: >> https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=9qKo982gn082WTnromnDXeqFdGEvh6A8v7Li-yEEiPs&e= >> ] > > > _______________________________________________ > fipy mailing list > [email protected] > https://urldefense.proofpoint.com/v2/url?u=http-3A__www.ctcms.nist.gov_fipy&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=dTROrUlTOg_pDStEAloQau0rkmx3ZuaHQy-NSohak0c&e= > > [ NIST internal ONLY: > https://urldefense.proofpoint.com/v2/url?u=https-3A__email.nist.gov_mailman_listinfo_fipy&d=DwIGaQ&c=nE__W8dFE-shTxStwXtp0A&r=GFdsMpHIw-aduHn9jOIpyQ&m=MiNELfmZGpKXQZs2VXW1BGZb0xHBiSGDtwc09TBIN_4&s=9qKo982gn082WTnromnDXeqFdGEvh6A8v7Li-yEEiPs&e= > ] > _______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
