MWR solutions of the steady advection-diffusion equation

Let us explore the steady-state advection diffusion equation and find approximate solutions to it using the MWR. Find \(u(x)\) that satisfies

(67)\[\begin{split}\begin{align} \begin{split} \frac{\partial u}{\partial t} = K \frac{\partial^2 u}{\partial x^2} - c\frac{\partial u}{\partial x}=0\\ c\frac{\partial u}{\partial x} - K \frac{\partial^2 u}{\partial x^2} =0\\ u(0)=0\\ u(1)=1 \end{split} \end{align}\end{split}\]

A scalar quantity \(u\) (e.g. temperature) is advected with velocity \(c\) and undergoes a diffusion dependent on K. The exact solution of this problem is:

(68)\[u^{ex}(x)=C_1 \exp \left ( \frac{c}{K} x \right ) + C_2\]

By using the boundary conditions, we can find the integration constants.

(69)\[\begin{split}\begin{align} \begin{split} C_1 &= -C_2\\ C_2 &= \frac{1}{1-\exp \left (\frac{c}{K} \right )} \end{split} \end{align}\end{split}\]

Tip

Try to derive the constants yourself and check that \(u^{ex}\) really is the exact solution!

Let’s solve by the Method of Weighted Residuals using a polynomial function as a basis. That is, let the approximating function be

(70)\[\tilde{u}(x)=a_0 + a_1 x + a_2 x^2\]

Application of the boundary conditions reveals

(71)\[\begin{split}\begin{align} \begin{split} \tilde{u}(0) &= 0=a_0\\ \tilde{u}(1) &= 1=0+a_1 + a_2\\ a_1 &= 1-a_2 \end{split} \end{align}\end{split}\]

so that the approximating polynomial which also satisfies the boundary conditions is

(72)\[\tilde{u}(x)=(1-a_2)x+a_2x^2\]

To find the residual \(R(x)\), we need the first and second derivative of this function, which are

(73)\[\begin{split}\begin{align} \begin{split} \frac{\partial \tilde{u}}{\partial x} &= (1-a_2)+2a_2 x\\ \frac{\partial^2 \tilde{u}}{\partial x^2} &= 2a_2 \end{split} \end{align}\end{split}\]

So the residual is

(74)\[R(x)=c\left ( (1-a_2)+ 2a_2 x \right ) - K(2a_2)\]

Collocation method

For the collocation method, the residual is forced to zero at a number of discrete points. Since there is only one unknown \(a_2\), only one collocation point is needed. We choose (arbitrarily, but from symmetry considerations) the collocation point \(x = 0.5\). Thus, the equation needed to evaluate the unknown \(a_2\) is

(75)\[R\left(\frac{1}{2}\right)=c\left ( (1-a_2)+ 2a_2 \left(\frac{1}{2}\right) \right ) - K(2a_2)=0\]

and the contant \(a_2\) is

(76)\[a_2=\frac{c}{2K}\]

Least-squares method

The weight function \(W_i\) with respect to the unknown \(a_2\):

(77)\[W_1(x)=\frac{dR}{da_2}=-c+2xc-2K\]

So the weighted residual statement becomes

(78)\[\begin{split}\int_{0}^{1}W_{1}(x)\cdot R(x)dx &= 0 \\ W_1(x)&=\frac{dR}{da_2}=c(2x-1)-2K \\ \int_{0}^{1}\left(c(2x-1)-2K \right )\left[c \left( (1-a_2)+2a_2 x \right ) -K(2a_2)\right ] dx &= 0\end{split}\]

the math is considerably more involved than before, but nothing more than integration of polynomial terms. Direct evaluation leads to the algebraic relation

(79)\[a_2 = \frac{6cK}{\left (c^2 + 12K^2 \right )}\]

Galerkin method

In the Galerkin Method, the weight function \(W_1\) is the derivative of the approximating function \(\tilde{u}\) with respect to the unknown coefficient \(a_2\):

(80)\[W_{1}(x)=\frac{d\widetilde{u}}{da_{2}}=x^{2}-x\]

So the weighted residual statement becomes

(81)\[\begin{split}\int_{0}^{1}W_{1}(x)\cdot R(x)dx &=&0 \\ \int_{0}^{1}\left( x^{2}-x\right) \cdot \left[ c\left( (1-a_2)+2a_2x\right)-K(2a_2) \right] dx &=&0\end{split}\]

Again, the math is straightforward but tedious. Direct evaluation leads to the algebraic equation:

(82)\[a_2=\frac{c}{2K}\]

Note how this solution is, for this special case, exactly the same as the one for the collocation method.

Derivation of the coefficients

We can derive the coefficients for the different methods by using the sympy symbolic math package for python. Install into your course environment:

conda activate "your environment"
conda install sympy

Notebook - derivation of constants

Excercise: Go through the notebook and try to re-derive the functional forms of the \(a_2\) coefficient!

Discussion

Let us evaluate the different approximations by exploring an example. If both the velocity and the diffusion constants are set to 1, the exact solution looks like in Figure 1. We can plot the approximate solutions using the equations above or by quickly re-deriving them using Python/Matlab. For this Matlab’s Symbolic Math Toolbox, sympy for Python are extremely useful.

Before we start let us introduce a useful measure for the quality of the approximation. A reasonable scalar index for the closeness of two functions is the L2 norm, or Euclidian norm. This measure is often called the root-mean-squared (RMS) error in engineering. The RMS error can be defined

(83)\[E_{RMS}=\frac{\sqrt{\int \left( u(x)-\widetilde{u}(x)\right) ^{2}dx}}{\int dx% }\]

which in discrete terms can be evaluated as

(84)\[E_{RMS}=\sqrt{\frac{\sum_{i=1}^{N}\left( u_{i}-\widetilde{u}_{i}\right) ^{2}% }{N}\text{.}}\]

where N is the number of grid points.