Tips on solving coupled system

[icojb25]

Firstly, thanks for t great tool, this is truly invaluable for those who work in such an area, and for making it available.

Second, I would like to solve the following system:

$$ \frac{\partial \rho}{\partial t} + \frac{\partial G}{\partial z} = 0 $$

$$ \frac{\partial G}{\partial t} + \frac{\partial}{\partial z}(\frac{G^2}{\rho}) + \frac{\partial P}{\partial z} + \frac{f}{D_H}\frac{|G|G}{2 \rho} = 0 $$

$$ \frac{\partial \rho h}{\partial t} + \frac{\partial G h}{\partial z} = \frac{D P}{D t} $$

From the paper, the linearised system is, where the * quantities are from the previous timestep:

$$ \frac{\partial \rho^*}{\partial t} + \frac{\partial G}{\partial z} = 0 $$

$$ \frac{\partial G}{\partial t}+\frac{\partial}{\partial z}(\frac{G^{* 2}}{\rho^*})+\frac{\partial P}{\partial z}+\frac{f}{D_H}\frac{G^2 *}{2 \rho *} = 0 $$

$$ \frac{\partial \rho^* h}{\partial t} + \frac{\partial G^* h}{\partial z} = \frac{D P}{D t} $$

From the papers i have seen (solved with another method, Least squares spectral method (LSSM)), the solution vector is $[G, P, h]$ so i am guessing $\rho$ is calculated using a gas law or an equation of state ... I was wondering whether this could be solved using Fipy - I guess it could be, and whether one needs to solve in a coupled or segregated fashion (sweeps), and whether anyone had any tips on how to get started.

I guess the transient terms (TransientTerm) are trivial, and the gradient terms can be modelled using the .grad or ImplicitSourceTerm or similar ...

EDIT. Fix up wrong derivative term (should $\frac{\partial G}{\partial z}$ should have been $\frac{\partial G}{\partial t}$ and $\frac{\partial}{\partial z}(\frac{G^2}{\partial \rho})$ should have been $\frac{\partial}{\partial z}(\frac{G^2}{\rho})$ and also in the linearised equations

[guyer]

So $\rho$ is a density or something, and not a radial coordinate? From some of the terms, I initially though this was in cylindrical coordinates. Can you check the expression:
$$\frac{G^2}{\partial \rho}$$
Should this be
$$\frac{\partial G^2}{\partial \rho}$$
or is $\partial$ a variable?

As to the forms of the terms, I try to avoid 1D expressions because the role of the derivatives is ambiguous:

• Are $G$ and $P$ scalar or vector?
• Is $\partial G/\partial z$ equal to $\nabla G \equiv \partial_i G$ or is it equal to $\nabla \cdot G \equiv \partial_i G_i$? Even if $G$ is known to be a vector field, $\nabla \vec{G} \equiv \partial_j G_i$ is not the same operation as $\nabla \cdot \vec{G} \equiv \partial_i G_i$.

It may technically not matter for the 1D analytical expressions, but the discretizations are very different and they're not the same in higher dimensions at all.

If I had to guess, $G$ is vector and $\partial G/\partial z \equiv \nabla\cdot \vec{G}$. This would be represented with a ConvectionTerm by FiPy.

If $P$ is pressure, then it should be scalar and $\partial P/\partial z \equiv \nabla P$. This would probably use .grad.

The problem with this interpretation of $G$ and $P$ is that

$$ \begin{align} \frac{\partial G}{\partial z} + \frac{\partial}{\partial z}(\frac{G^2}{\partial \rho}) + \frac{\partial P}{\partial z} + \frac{f}{D_H}\frac{|G|G}{2 \rho} &= 0 \\ \nabla\cdot G + \nabla\stackrel{?}{\cdot}(\frac{G^2}{\partial \rho}) + \nabla P + \frac{f}{D_H}\frac{|G|G}{2 \rho} &= 0 \\ \text{scalar} + \text{?} + \text{vector} + \text{vector} &= 0 \end{align} $$

is dimensionally inconsistent (I really wish that \underbrace{} worked properly in GFM math).

If $G$ and $P$ are both scalar, then

$$ \begin{align} \frac{\partial G}{\partial z} + \frac{\partial}{\partial z}(\frac{G^2}{\partial \rho}) + \frac{\partial P}{\partial z} + \frac{f}{D_H}\frac{|G|G}{2 \rho} &= 0 \\ \nabla G + \nabla(\frac{G^2}{\partial \rho}) + \nabla P + \frac{f}{D_H}\frac{|G|G}{2 \rho} &= 0 \\ \text{scalar} + \text{scalar} + \text{scalar} + \text{scalar} &= 0 \end{align} $$

is OK, but

$$ \begin{align} \frac{\partial \rho}{\partial t} + \frac{\partial G}{\partial z} &= 0 \\ \frac{\partial \rho}{\partial t} + \nabla G &= 0 \\ \text{scalar} + \text{vector} &= 0 \end{align} $$

is not.

So, I think $G$ and $P$ need to both be vector. The equations would all be dimensionally consistent, but then $\vec{P} \equiv p \mathsf{I}$ is the volumetric stress contribution to the Cauchy stress tensor, rather than the pressure, which is rather confusing notation.

No good ever comes from writing PDEs in 1D.

Another point of confusion is the interpretation of $D P/ Dt$. This looks like a material derivative, but what is its flow velocity?

[icojb25]

Hi there,
Thank you for the reply, its much appreciated. Apologies for the typos:

So ρ is a density or something, and not a radial coordinate? From some of the terms, I initially though this was in cylindrical coordinates. Can you check the expression: G2∂ρ Should this be ∂G2∂ρ or is ∂ a variable?

You are correct, this was a typo, it should have been $\frac{\partial}{\partial z}(\frac{G^2}{\rho})$. I should also have indicated the meanings of the terms.

• $\rho$ is density
• $G$ is mass flux
• $P$ is static pressure and
• $h$, $f$, $D_h$ are enthalpy, friction factor, and hydraulic diameter
• $z$ and $t$ are the spatial and time coordinates

As to the forms of the terms, I try to avoid 1D expressions because the role of the derivatives is ambiguous:

• Are G and P scalar or vector?

I guess $G$ is effectively scalar since its a 1D flux, $P$ is scalar (pressure)

• Is ∂G/∂z equal to ∇G≡∂iG or is it equal to ∇⋅G≡∂iGi? Even if G is known to be a vector field, ∇G→≡∂jGi is not the same operation as ∇⋅G→≡∂iGi.

It may technically not matter for the 1D analytical expressions, but the discretizations are very different and they're not the same in higher dimensions at all.
If I had to guess, G is vector and ∂G/∂z≡∇⋅G→. This would be represented with a ConvectionTerm by FiPy.

Ok, thanks a lot

If P is pressure, then it should be scalar and ∂P/∂z≡∇P. This would probably use .grad.

Ok, again, thanks

The problem with this interpretation of G and P is that

∂G∂z+∂∂z(G2∂ρ)+∂P∂z+fDH|G|G2ρ=0∇⋅G+∇⋅?(G2∂ρ)+∇P+fDH|G|G2ρ=0scalar+?+vector+vector=0

is dimensionally inconsistent (I really wish that \underbrace{} worked properly in GFM math).

If G and P are both scalar, then

∂G∂z+∂∂z(G2∂ρ)+∂P∂z+fDH|G|G2ρ=0∇G+∇(G2∂ρ)+∇P+fDH|G|G2ρ=0scalar+scalar+scalar+scalar=0

is OK, but

∂ρ∂t+∂G∂z=0∂ρ∂t+∇G=0scalar+vector=0

is not.

So, I think G and P need to both be vector. The equations would all be dimensionally consistent, but then P→≡pI is the volumetric stress contribution to the Cauchy stress tensor, rather than the pressure, which is rather confusing notation.

No good ever comes from writing PDEs in 1D.

:)

Another point of confusion is the interpretation of DP/Dt. This looks like a material derivative, but what is its flow velocity?

Yes, the paper doesnt actually say, but i assumed it was. These equations are 1D models of mass, momentum and energy convservation (modelling pressure and density wave instabilities in thermohydraulic systems)