finite volume coupled equations solving
 

Hello,

How can i solve a coupled two equation system (source term coupling)
the unknow variable are Tf and Ts

k1 (dTf)/dt+k2 (dTf)/dx=k3 (d² Tf/dx²) -k4(Tf-Ts)

k5 (dTs)/dt=k6 (d²Ts/dx²) +k4(Tf-Ts)

After discretising (with fully implicit scheme) :

af,i*Tf,i^(n+1)=bf,i* Tf,i+1^(n+1)+cf,i *Tf,i-1^(n+1)+df,i* Tf,i^n+bf Ts,i^(n+1)

as,i *Ts,i^(n+1)=bs,i *Ts,i+1^(n+1)+cs,i*Ts,i-1^(n+1)+ds,i* Ts,i^n+bs *T(s,i)^(n+1)

i: is for spatial discretisation and n is for the temporal one.

If this is just a problem in one spatial dimension, you can use a block Thomas algorithm to solve the fully coupled linear system very efficiently.

For "1x1 blocks" this reduces to the normal Thomas algorithm or Tridiagonal Matrix Algorithm for 3-point stencils. Just pack the unknowns into the unknown vector cell-by-cell [Tf1,Ts1,Tf2,Ts2...Tfn,Tfs]. And likewise, pack the coefficients into 2x2 blocks. The as and af entries go on the diagonal alternating as the unknowns do in solution vector. But there are a_sf and a_fs coefficients too that complete the block...those arise from the k4(Tf-Ts) term. It is not that complicated, but write out a small...say 4-cell coefficient matrix (yielding and 8x8 matrix) and you will see quickly how things line up.

You could also just use a banded-linear equation solver (a Thomas algorithm with a wider stencil). Data would be packed the same way, but your matrix bandwidth would be 7 (versus the normal TMA's bandwidth of 3), but the solution is done exactly the same. This would be slightly less efficient than the blocked Thomas algorithm because there would be some zero fill-in. You will see it when you if you write out the combined matrix entries by hand.

If the problem is not just one-spatial dimension, go find a pre-conditioned bi-conjuguate gradient stabilized function, learn its interface, generate the coefficient matrix exactly as above, but push the coefficients into the (likely) compressed row storage format. You can likely just use Jacobi preconditioning, ILU would be better if it is available.

Thank you Michael for your answer :)

I set the matrix as you explained and I see what you mean.
I got something like; coef(1,1)= afp,1 ; coef(1,2)= -bf ; coef(1,3)= afE,1 ; Other coefficients in the first line are equal to 0.
For the second line, coef(2,1)= -bs ; coef(2,2)= asp,1 ; coef(2,3)= 0 ; coef(2,4)= asE,1 ; the Other coefficients in the second line are equal to 0.
And so…

So I will have the nulls values in the extremities of the matrix ( as in ordinary tridiagonal matrix) and in addition, the coefficients: (i,i+1) and (i+1,i) are equal to zero (from i=2), is that with no influence to solve the system with a Thomas algorithm?

Best regards!

If you keep going, you will see that away from the boundaries, there will be even wider entries. Let me work with a more simple version of your equations:

Tf[i-1] - 2.0*Tf[i] + Tf[i+1] + C*(Ts[i] - Tf[i]) = 0
Ts[i-1] - 2.0*Ts[i] + Ts[i+1] + D*(Tf[i] - Ts[i]) = 0

So, for rows 2i and 2i+1, the entries will be:
You can see that the source term coupling results in a sparse matrix with bandwidth of 5, so you can use you use a 5-wide bandwidth Thomas algorithm...it is based on Gauss Elimination just like the normal 3-wide version, but with one more band of elimination and an extra term in the back substitution. Note that in the middle, you COULDl have six potentially non-zero entries, but source term coupling doesn't fill out all six. If you had an equation with both laplacian(Ts) and laplacian(Tf), it would fill all six entries.

There are also cyclic redundancy techniques that may be slightly more efficient. But for most 1D transient simulations, the performance of any of these techniques will be more than adequate.

Here is a link to a pentadiagonal solver document:

http://www.hindawi.com/journals/mpe/2015/232456/

Here is link to some sample FORTRAN code:

http://www.math.uakron.edu/~kreider/anpde/penta.f

Thank you for your precious help :)