Capacitance in TransientReactiveTransport

[shillcoat]

Hello,

I'm an MME student at UW doing research at the MPI lab. I've been working on modifying OpenPNM for simulations of the porous flow in transpiration cooling (the repository is here if you're interested). The work is specifically looking at the impact of surface pressure fluctuations, so I've been building on the TransientReactiveTransport algorithm to implement what is effectively a Stokes flow algorithm with changing boundary conditions.

I noticed that both the implicit and Crank-Nicholson time schemes introduce a capacitance to the problem due to the way the A and b matrices are updated:

def _t_update_A(self):
        r"""
        A method to update 'A' matrix at each time step according to 't_scheme'
        """
        network = self.project.network
        Vi = network[self.settings['pore_volume']]
        dt = self.settings['t_step']
        A = self._f1 * self._A_steady
        A.setdiag(A.diagonal() + self._f2/dt * Vi)
        self._A = A

def _t_update_b(self):
        r"""
        A method to update 'b' array at each time step according to
        't_scheme' and the source term value
        """
        quantity = self.settings['quantity']
        network = self.project.network
        phase = self.project.phases()[self.settings['phase']]
        Vi = network[self.settings['pore_volume']]
        dt = self.settings['t_step']
        X = self[quantity]
        b = self._f2 * (1-self._f1) * -self._A_steady * X + self._f2 * Vi/dt * X
        self._update_iterative_props()
        for item in self.settings['sources']:
            Ps = self.pores(item)
            b[Ps] -= self._f2 * (1-self._f1) * phase[f"{item}.rate"][Ps]
        self._b = b

I was hoping someone would be able to explain or point me to a paper that explains what equations the modifications to A and b correspond to, as it's clear the transient algorithm builds further on steady state algorithms than just implementing basic time stepping. Any insight you could provide would be greatly appreciated.

Thank you,
Sophie

[jgostick]

Hi Sophie, I'm sorry I can't answer this right away. Also, since your at UW maybe we should just meet in person? I can say that we have update the transient algorithms quite a bit in the upcoming V3 release, including a explicit solvers and the use of scipy's ivp_solve module.

[shillcoat]

Hi Dr. Gostick,

No problem, I understand you're probably quite busy. Meeting in person may be helpful and more efficient if it's not too much trouble for you. My supervisor, Dr. Jean-Pierre Hickey, would also like to join, which I assume is not an issue. I will send an email so we can coordinate this.

I did notice that the transient solver has been updated quite a bit on the dev branch, I'll take a closer look. For the moment I've been working from the latest release (v2.8.2) because of the high backwards incompatibility of v3, but perhaps I am better migrating over.

Thank you,
Sophie