Coupled PDEs with dependent sources formulation

[dgallegos47]

Please forgive what may be a basic question but I am looking for some guidance on the correct approach of applying FiPy to an existing model. An existing model in the literature for describing polymer decomposition is described by:

$$\frac{\partial T}{\partial t} = \frac{1}{\rho c_p}\left(\frac{\partial}{\partial y}\left(\kappa\frac{\partial T}{\partial y}\right) - \lambda\rho\Delta h_{v}k-\rho c_p v \frac{\partial T}{\partial y}\right)\tag{1}$$

$$\frac{\partial \lambda}{\partial t} = -k\lambda + \frac{\lambda}{\rho}\frac{\dot{m}^{''}_{v}}{\partial y} - v \frac{\partial \lambda}{\partial y}\tag{2}$$

$$\dot{m}_{v}^{''} = \beta(1-\lambda)\rho\frac{\partial T}{\partial y}\tag{3}$$

The material properties (coefficients) are dependent on the polymer mass fraction $\lambda(y, t)$ and need to be updated at each time step. Similarly, the rate constant $k$ is dependent on temperature according to an Arrhenius style equation and is therefore a function of cell temperature. The mass flux $\dot{m}^{''}_v$ is evaluated at each time step using the computed mass fraction and temperature field (eq. 1).

Generally, I was hoping for some guidance on the implementation. Considering Eq. 1, and trying to separate the terms in the transient-diffusion-convection-source format:

$$\rho c_p \frac{\partial T}{\partial t} + \rho c_p v \frac{\partial T}{\partial y} = \frac{\partial }{\partial y}\left(\kappa\frac{\partial T}{\partial y}\right) - \lambda \rho \Delta h_{v} k$$

How would I specify the non-unity coefficient of the convective term? Would I need to re-formulate such that the $\rho c_p$ term leading the convective term winds up in the source term such that the general form is respected? My primary question is regarding the computation of the non-constant material properties and other variables that are necessary for the computation of the above equations. Should I define them all as CellVariables() such they are handled (somehow) in the eqns.solve() or should I handle the update of the material properties and assignment independently?

There is some documentation for coeffcients that are dependent on the solution variable (i.e. $\rho = \rho(T)$), I'm unsure how to incorporate coeffcients that are dependent on the result of a second, coupled equation. Further, how can I incorporate the results from both equations into the computation of the mass flux (eq. 3) such that the next timestep incorporates the updated values for all fields.

Related, there was a post by @wd15 on the mailing list archive that discussed a problem that is relevant to the phase change - moving interface processes, does anyone know where I can find such an implementation for reference? Particularily interested in how the moving interfaces are treated.

[guyer]

How would I specify the non-unity coefficient of the convective term? Would I need to re-formulate such that the ρcp term leading the convective term winds up in the source term such that the general form is respected?

You'll need to work the chain rule to get all of the coefficients inside the derivatives where FiPy expects them, e.g., ConvectionTerm(coeff=rho * cp * v, var=T) represents $\frac{\partial}{\partial y}\left(\rho c_p v T\right)$. Applying the chain rule gives $\frac{\partial}{\partial y}\left(\rho c_p v T\right) \equiv \rho c_p v \frac{\partial T}{\partial y} + T \frac{\partial}{\partial y}\left(\rho c_p v\right)$. That second term on the right needs to be accounted for, possibly via further application of the chain rule.

Alternatively, dig into the derivation of your equations. I'm not a fluids person, but I've found that most presentations assume that $\rho$, $c_p$, or $\kappa$ do not vary in time or space, and then they are applied to situations where that's not true. I don't know what the correct formulation is, but it needs to account for all of the sources of variation of the coefficients.

My primary question is regarding the computation of the non-constant material properties and other variables that are necessary for the computation of the above equations. Should I define them all as CellVariables() such they are handled (somehow) in the eqns.solve() or should I handle the update of the material properties and assignment independently?

If you have a functional expression for them, then you should write that, e.g., if $k = k_0 \exp(-Q/k_B T)$, then write

k = k0 * fp.numerix.exp(-Q/(kB * T))

where T is the CellVariable you are solving for in (1) and k0, Q, and kB are defined as whatever they are, either constants or functions of variables. k will then be updated whenever T changes.

There is some documentation for coeffcients that are dependent on the solution variable (i.e. $\rho = \rho(T)$), I'm unsure how to incorporate coeffcients that are dependent on the result of a second, coupled equation. Further, how can I incorporate the results from both equations into the computation of the mass flux (eq. 3) such that the next timestep incorporates the updated values for all fields.

You can, but I would not, write an expression for mass flux, which will update whenever there are changes in the variables it depends on. I wouldn't do this, though, because, as written, Eq. (2) and Eq. (3) will have poor implicit properties. Instead, substitute (3) into (2) to get

$$ \frac{\partial \lambda}{\partial t} = -k \lambda + \frac{\lambda}{\rho}\frac{\partial}{\partial y}\left[\beta\left(1-\lambda\right)\rho\frac{\partial T}{\partial y}\right] - v\frac{\partial \lambda}{\partial y} \tag{2'} $$

This puts the second term in the form of a diffusion term in T, which will provide better coupling with Eq. (1).

Related, there was a post by @wd15 on the mailing list archive that discussed a problem that is relevant to the phase change - moving interface processes, does anyone know where I can find such an implementation for reference? Particularily interested in how the moving interfaces are treated.

See examples.reactiveWetting.liquidVapor1D and examples.reactiveWetting.liquidVapor2D.