Coupled Partial Differential Equations

Question

I'm trying to solve the following system of coupled differential equations, the two-temperature model for $e$ = electrons and $l$ = lattice.

$$ \rho_{e}C_{p,e}\frac{\partial T_{e}}{\partial t} = k_{e}\frac{\partial^{2}T_{e}}{\partial x^{2}} - G(T_{e}-T_{l}) + S(x,t) $$ $$ \rho_{l}C_{p,l}\frac{\partial T_{l}}{\partial t} = k_{l}\frac{\partial^{2}T_{l}}{\partial x^{2}} + G(T_{e}-T_{l}) $$

I was thinking of re-casting them into a form suitable for the Crank-Nicolson method, which would be the following if $T_{e}=u_{i}^{n}$ and $T_{l}=v_{i}^{n}$ where $i$ is the space node and $n$ is the time node.

$$ A_{o}u_{i+1}^{n+1} + B_{o}u_{i}^{n+1} + A_{o}u_{i-1}^{n+1} = A_{1}u_{i+1}^{n} + B_{1}u_{i}^{n} + A_{1}u_{i-1}^{n+1} + C_{1}v_{i}^{n} + D_{1}S_{i}^{n} $$ $$ A_{2}v_{i+1}^{n+1} + B_{2}v_{i}^{n+1} + A_{2}v_{i-1}^{n+1} = A_{3}v_{i+1}^{n} + B_{3}u_{i}^{n} + A_{3}v_{i-1}^{n+1} + C_{3}v_{i}^{n} $$

where $A_{o}=-\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{2dx^{2}}, \ A_{1} = -A_{o}, $ $B_{o}=1+\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{dx^{2}}, \ B_{1} = 1-\frac{dt}{\rho_{e}C_{p,e}}\frac{k_{e}}{dx^{2}} + G, \ C_{1} = \frac{dt}{\rho_{e}C_{p,e}}G, \ D_{1} = \frac{dt}{\rho_{e}C_{p,e}}$

$A_{2}=-\frac{dt}{\rho_{l}C_{p,l}}\frac{k_{l}}{2dx^{2}}, \ A_{3}=-A_{2}$ $B_{2} = 1+\frac{dt}{\rho_{l}C_{p,l}}\frac{k_{l}}{2dx^{2}}, \ B_{3}=1-\frac{dt}{rho_{e}Cp_{e}}\frac{k_{l}}{2dx^{2}}+G, \ C_{3} = -\frac{dt}{\rho_{l}C_{p,l}}G$

Analytically, this works out but I'm drawing a blank on how to implement it (I want to use Python). I can do each of the equation systems by themselves easily enough, but I was trying to solve both in the same matrix system. I've seen certain attempts at this sort of problem by converting the two unknowns at the updated time into a vector $\vec{W}_{i}^{n+1}=(u_{i}^{n+1},v_{i}^{n+1})^{T}$ but due to the coefficients being different for each (apart from the coupling constant G) I wasn't sure how to go about it.

Answer 1

Good afternoon,

If the coefficients are correctly found (I have not checked), you can present the system as @lightxbulb says. Being $N$ the total number of space nodes, let me expand the answer a little bit for completeness, therefore

$u^{n+1} = \left[u_1^{n+1},\,u_2^{n+1},\,...,\,u_N\right]^\top$

$v^{n+1} = \left[v_1^{n+1},\,v_2^{n+1},\,...,\,v_N\right]^\top$

$b(u^n,v^n) = \left[b_1,\,b_2,\,...,b_N\right]^\top$

where $b_i,\forall i=2,...,N-1$ can be computed as Assuming $A_1u_{i-1}^{n+1}$ to be a typo in your question and be $A_1u_{i-1}^{n}$ instead, which makes more sense to me

$b_i = A_{1}u_{i+1}^{n} + B_{1}u_{i}^{n} + A_{1}u_{i-1}^{n} + C_{1}v_{i}^{n} + D_{1}S_{i}^{n}$

and $b_1$ and $b_N$ will depend on the boundary conditions of the problem.

$c(u^n,v^n) = \left[c_1,\,c_2,\,...,c_N\right]^\top$

where $c_i,\forall i=2,...,N-1$ can be computed as

$c_i = A_{3}v_{i+1}^{n} + B_{3}u_{i}^{n} + A_{3}v_{i-1}^{n+1} + C_{3}v_{i}^{n}$

and $c_1$ and $c_N$ will depend on the boundary conditions of the problem.

Finally, we have that

$M = \begin{bmatrix} B_0&A_0&0&\ldots & 0 & 0& 0\\ A_0 &B_0&A_0&\ldots& 0 & 0 & 0\\ \vdots &\vdots&\vdots&\ddots & \vdots & \vdots & \vdots\\ 0 &0&0&\ldots& A_0 &B_0&A_0\\ 0 &0&0&\ldots& 0 &A_0&B_0\\ \end{bmatrix} $

$N = \begin{bmatrix} B_2&A_2&0&\ldots & 0 & 0& 0\\ A_2 &B_2&A_2&\ldots& 0 & 0 & 0\\ \vdots &\vdots&\vdots&\ddots & \vdots & \vdots & \vdots\\ 0 &0&0&\ldots& A_2 &B_2&A_2\\ 0 &0&0&\ldots& 0 &A_2&B_2\\ \end{bmatrix} $

As @lightxbulb says, it becomes a decoupled system maybe you don't want to threat the system matrix as a whole. You can use a huge variety of methods for solving this linear system to obtain the vector of unkowns you mentioned

$\vec{W}^{n+1} = \begin{bmatrix} u^{n+1}\\v^{n+1} \end{bmatrix}$

I hope it helps,