How to get discretization coefficients of matrix A in Finite Volume Method (FVM)?

Question

First we have Discretization of the Transport Equation $$ \frac{\partial \rho \phi}{\partial t} + \nabla(\rho U \phi) - \nabla (\rho \Gamma_\phi \nabla \phi) = S_\phi (\phi) $$

In Finite Volume Method it looks like: $$ \int_t^{t+\Delta t} \left[ \frac{\partial}{\partial t} \int_{V_p} \rho \phi dV + \int_{V_p} \nabla \cdot (\rho U\phi)dV - \int_{V_p} \nabla \cdot (\rho \Gamma_\phi \nabla \phi)dV \right]dt = \int_t^{t+\Delta t} \left(\int_{V_p}S_\phi(\phi)dV \right)dt $$

and after transformations:

$$ \int_t^{t+\Delta t} \left[ \left( \frac{\partial \rho \phi}{\partial t} \right)_p V_p + \sum_f F\phi_f - \sum_f(\rho \Gamma_\phi)_f S.(\nabla\phi)_f \right]dt=\int_t^{t+\Delta t}(SuV_p+S_pV_p\phi_p)dt $$

and time discretization: $$ \frac{\rho_p\phi_p^n-\rho_p\phi_p^o}{\Delta t}V_p+ \frac{1}{2}\sum_f F\phi_f^n-\frac{1}{2}\sum_f(\rho\Gamma_\phi)_f S.(\nabla\phi)^n_f+ \frac{1}{2}\sum_f F\phi_f^o-\frac{1}{2}\sum_f(\rho\Gamma_\phi)_f S.(\nabla\phi)^o_f = SuV_p + \frac{1}{2}S_pV_p\phi_p^n+\frac{1}{2}S_pV_p\phi_p^o $$

For every cell we can make equation:

$$ a_p\phi^n_p+\sum_N a_N \phi_N^n = R_p $$

But how we can get elements of matrix $A$ if we don't know $\phi$? $$ [A][\phi]=[R] $$

Here Hrvoje Jasak wrote that every coefficient $ a_p $ includes the contribution from temporal derivative, convection and diffusion terms. But what formula of $a_p = ...$?

http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdf

Answer 1

The second to last equation you have in your question is linear in the $\phi^n_p$, and since you have one equation for each $p$, this is a linear system. It all comes down to just identifying which $\phi^n_p$ appears in which equation.

I don't have this for finite volumes, but if you want to see it happen for finite elements, you may want to watch lecture 4 here: http://www.math.tamu.edu/~bangerth/videos.html

Answer 2

First of all you know $\phi$ from field initialization or from previous iteration(or timestep). When you use field values in such way, considering mass fluxes known at the beggining of each iteration, you create a decoupling in nonlinear term, and Jasak mentioned that in his thesis, also saying that it causes an error, which vanishes when iterations converge. Also he mentioned that this decoupling is less severe than pressure-velocity decoupling.

As for the main diagonal coefficient $a_P$ it is negative sum of neighboring coefficients $a_{nb}$ plus implicit part of non-stationary term. Physically it means that flux going out from one cell, goes into another, neighbor cell (share the same cell face).

Check out Ferziger & Peric book (a recommendation).

Answer 3

The $a_p$ is assembled from all the coefficients that multiply $\phi_p^n$.

Starting with the scalar transport equation in the discretized form

$\dfrac{\frac{3}{2}\rho_P\phi_P^n - 2\rho_P\phi_P^o + \frac{1}{2}\rho_P \phi_P^{oo}}{\Delta t} + \sum_{f}F\phi_f^n - \sum_{f}(\rho\gamma)_f(\nabla{\phi})_f^n = S_P^o $

the key point to understand where the $a_p$ comes from are the face values $\phi_f^n$, that will be interpolated from the values between the neighboring cells. Looking at the equation above, you can already see that you have the $\phi_P^n$ value in the temporal term, but the interpolation of the values from the new time step will lead to the final implicit algebraic equation for cell $P$, and more $\phi_P^n$ terms will appear.

Say you choose the weighted Central Differencing Scheme (CDS) to interpolate $\phi_P^n$ for all $\phi_f^n$ of the above equation, leading to

$\phi_f^n = f_x \phi_P^n + (1 - f_x) \phi_N^n$,

where $f_x$ is the CDS weight $P$ is the owner cell and N is the face neighbor cell. For every occurence of $\phi_f^n$ in the first equation, the interpolation will introduce a $\phi_P^n$ and $\phi_N^n$. If you do this by hand and sort out everything that multiplies $\phi_P^n$, you'll get the "formula" for $a_p$.

Another important thing is also to understand that $P$ and $N$ exchange interpretations based on the orientation of the surface are normal vector of the finite volume face. This orientation is defined by the mesh generation step, and is used to optimize the cell sum calculation - check out Hrv's thesis for more info on that. This has no impact in the symbolic calculation of $a_p$, it's an optimization enhancement.