How to handle the evaluation of functions on staggered ghost nodes?

20 Views Asked by At

I have a convection-diffusion-reaction steady state PDE in the form $$ u(x,z)\frac{\partial C}{\partial x} = \frac{\partial}{\partial z} \left( K(z) \frac{\partial C}{\partial z} \right) - \lambda C - w(z) \frac{\partial C}{\partial z} $$ on $(x,z)\in[0,X]\times[0,H]$ with the following boundary conditions on three of the four boundaries $$ C(0,z) = C_0, \phantom{A}\frac{\partial C}{\partial z}(x,H) = 0, \phantom{A}\textrm{and}\phantom{A} K(0) \frac{\partial C}{\partial z}(x,0) = aC - b(x) $$ I wanted to construct a somewhat conservative form finite difference scheme for the form $$ u(x,z)\frac{\partial C}{\partial z} = \frac{\partial}{\partial z} \left( K(z) \frac{\partial C}{\partial z} - w(z)C \right) + \left(\frac{\partial w}{\partial z} - \lambda \right) C $$ I can apply the Crank-Nicolson scheme for the left-hand side so \begin{multline*} \left[ u(x,z)\frac{\partial C}{\partial x} \right]_{i+1/2,j} \\= \frac{1}{2}\left( \left[ \frac{\partial}{\partial z} \left( K(z) \frac{\partial C}{\partial z} - w(z)C \right) + \left(\frac{\partial w}{\partial z} - \lambda \right) C \right]_{ij} + \left[\frac{\partial}{\partial z} \left( K(z) \frac{\partial C}{\partial z} - w(z)C \right) + \left(\frac{\partial w}{\partial z} - \lambda \right) C \right]_{i+1,j} \right) \end{multline*} if we let $$ q = K(z) \frac{\partial C}{\partial z} - w(z)C \phantom{A}\textrm{and}\phantom{A} R = w_j' -\lambda $$ from this we can write $$ \left[ \frac{\partial q}{\partial z} + R C \right]_{ij} = \left[ \frac{\partial q}{\partial z} \right]_{ij} + R_j C_{ij} $$ next we can write the partial derivative of $C$ in terms of a staggered grid $$ \left[ \frac{\partial q}{\partial z} \right]_{ij} = \frac{q_{i,j+1/2}-q_{i,j-1/2}}{\Delta z} $$ Now the two following terms are \begin{align*} q_{i,j+1/2} &= K_{j+1/2} \left[\frac{\partial C}{\partial z} \right]_{i,j+1/2}-w_{j+1/2} C_{i,j+1/2} = K_{j+1/2} \left( \frac{C_{i,j+1} - C_{ij}}{\Delta z}\right) - w_{j+1/2} \left( \frac{C_{i,j+1}+C_{ij}}{2} \right)\\ q_{i,j-1/2} &= K_{j-1/2} \left[\frac{\partial C}{\partial z} \right]_{i,j-1/2}-w_{j+1/2} C_{i,j-1/2} = K_{j-1/2} \left( \frac{C_{ij} - C_{i,j-1}}{\Delta z}\right) - w_{j-1/2} \left( \frac{C_{ij}+C_{i,j-1/2}}{2} \right) \end{align*} Next we can write \begin{multline*} \left[ \frac{\partial q}{\partial z} \right]_{ij} = \frac{1}{2\Delta z^2} \left( \left(2K_{j+1/2} - \Delta z \cdot w_{j+1/2}\right) C_{i,j+1}\right. \\ -\left.\left( 2\left( K_{j+1/2} + K_{j-1/2}\right) + \Delta z\left( w_{j+1/2} - w_{j-1/2}\right)\right) C_{ij} + \left(2K_{j-1/2} + \Delta z \cdot w_{j-1/2}\right) C_{i,j-1} \right) \end{multline*} Next, we can write the following for the reaction term $$ R_jC_{ij} = \frac{R_j}{2}\left( C_{i,j+1} + C_{i,j-1} \right) $$ Now we can write \begin{multline*} \left[ \frac{\partial q}{\partial z} \right]_{ij} + R_j C_{ij} = \frac{1}{2\Delta z^2} \left( \left(2K_{j+1/2} - \Delta z \cdot w_{j+1/2} + \Delta z^2 R_j\right) C_{i,j+1}\right. - \\ \left.\left( 2\left( K_{j+1/2} + K_{j-1/2}\right) + \Delta z\left( w_{j+1/2} - w_{j-1/2}\right)\right) C_{ij} + \left(2K_{j-1/2} + \Delta z \cdot w_{j-1/2} + \Delta z^2R_j\right) C_{i,j-1} \right)\\ = \mu(\alpha_j C_{i,j+1} + \beta_j C_{ij} + \gamma_j C_{i,j-1}) \end{multline*} Going back to the full scheme, I can now write $$ \left[ u(x,z)\frac{\partial C}{\partial x} \right]_{i+1/2,j} = \frac{\mu}{2}(\alpha_j C_{i,j+1} + \beta_j C_{ij} + \gamma_j C_{i,j-1}) + \frac{\mu}{2}(\alpha_j C_{i+1,j+1} + \beta_j C_{i+1,j} + \gamma_j C_{i+1,j-1}) $$ The left-hand side of this equation can be written as $$ \left[ u(x,z)\frac{\partial C}{\partial x} \right]_{i+1/2,j} = u_{i+1/2,j} \left(\frac{C_{i+1,j} - C_{ij}}{\Delta x}\right) $$

Next I can move everything over to the right-hand side and write $$ 0 = \frac{\mu \Delta x}{2} (\alpha_j C_{i,j+1} + (\beta_j - \frac{2}{\mu \Delta x}u_{i+1/2,j}) C_{ij} + \gamma_j C_{i,j-1}) + \frac{\mu\Delta x}{2}(\alpha_j C_{i+1,j+1} + (\beta_j - \frac{2}{\mu \Delta x}) C_{i+1,j} + \gamma_j C_{i+1,j-1}) $$ I can multiply the factor on the outside in and get the new expression, as well as moving the terms at $i$ to the left to get $$ A_j C_{i,j+1} + B_{ij} C_{ij} + D_j C_{i,j+1} = E_j C_{i+1,j+1} + F_{ij} C_{i+1,j} + G_{j} C_{+i+1,j-1} $$ For convenience, I'll let the \begin{gather*} A_0 = -\frac{ \mu \Delta x}{2}\left(\alpha_0 + T \gamma_0 \right)\\ D_n = -(\alpha_n + \gamma_n)\\ E_0 = \frac{ \mu \Delta x}{2}\left(\alpha_0 + T \gamma_0 \right)\\ G_n = \alpha_n + \gamma_n \end{gather*} Where $$ T = \frac{K_0 - a\Delta z}{K_0 + a\Delta z} $$ which comes from the boundary conditions (I can re-edit to show the derivation if other's find it necessary).\ Lastly, from the boundary conditions I finish with writing $$ A_j C_{i,j+1} + B_{ij} C_{ij} + D_j C_{i,j+1} - S_{ij} = E_j C_{i+1,j+1} + F_{ij} C_{i+1,j} + G_{j} C_{+i+1,j-1} + S_{i+1,j} $$ where $$ S_{ij} = \begin{cases} \gamma_0\frac{2\Delta z}{K_0 + a \Delta z}b_i & \text{if } j = 0\\ 0 & \text{otherwise} \end{cases} $$ This whole writing shows that I can solve the matrix equation $$ U \vec{C}_{i} - \vec{S}_i = V \vec{C}_{i+1} + \vec{S}_{i+1} $$ The issue is that in this problem, when evaluating the $\gamma_0$ or $\beta_0$, I must evaluate $K_{-1/2}$ and $w_{-1/2}$ and my functions are $$ K(z) \propto z^{p_1} \phantom{A}\text{and}\phantom{A} w(z) \propto z^{1-p_1} $$ where $p_1 \in (0,1)$ and it is expected that $z \geq 0$. How would I go about justifying some kind of evaluation of $K_{-1/2}$? I was reading up on Ghost Fluid Methods as they seem to pertain to problems with interfaces, which I believe the boundary of the problem itself is considered to be an interface between the domain and whatever is outside the domain.