Effect of mesh size on solution curves for a 1D problem

Question

I'm interested in studying the effect of mesh size on the behavior of the solution curves of 1D convection-diffusion problem.

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2} - v\frac{\partial C}{\partial x}$$

I could use MATLAB's pdepe solver to solve the above system.

function sol=so()
format short
global D nnode init_co find_index v
m = 0;
delx = 0.25;
xend = 5; 
D = 500;
v = 200;
x = 0:delx:xend; % mesh
find_index  = 0:delx:xend;
t = 0:0.00001:0.5;
init_co = [1 ; zeros(length(x)-1,1)];
nnode = length(x);
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
plot(t,sol)
xlabel('time')
ylabel('c')
grid on
function [g,f,s] = pdefun(x,t,c,DcDx)
g = 1;
f = D*DcDx;
s = -v*DcDx;
end

function c0 = icfun(x)
c0 = init_co(find(find_index==x));
end

function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = cl - 3;
    ql = 0;
    pr = 0;
    qr = 1;
end
end

I changed delx = 1, the following is observed

By just observing the above plots, I find the behavior (i.e the steepness of the curves) is the same and doesn't vary much with mesh size.

I think it's a crude way to observe just the steepness of the curves. I'd like to ask for suggestions on other quantitative measures that can be computed from solution curves to study the effect of the mesh size.

EDIT: I've compared one quantity(absolute error) from the suggestions received below.

I understand the error increases with an increase in mesh size. I'd like to know how to strike a balance -- fine mesh will increase computational cost and coarse mesh will produce errors and even lead to oscillations in solution curves for high delx. For instance, in the literature, I came across the CFL number, which should be less than 1 to obtain stable solutions. However, I'm not sure how to compute the CFL number for my system since I don't know the step size that the solver(ode15s) takes. Also, I'd like to know what is the acceptable range of error (what percent of deviation from true solution is acceptable? is there a thumb rule?).

Answer 1

I'm a bit surprised that nobody gave you an analytical solution for comparing your simulation results with it and measure your mesh convergence. It's always good to specify the boundary conditions explicitly, because it's really difficult to understand it from your code unless you dig deep into the MATLAB documentation of pdepe. Your boundary conditions are:

$$C(0,t) = C_{b}$$

and

$$-D \frac{\partial C}{\partial x}|_{x=L} = 0$$

and you have this advection-diffusion equation:

$$\frac{\partial C}{\partial t} = D \frac{\partial^{2} C}{\partial x^{2}} - v \frac{\partial C}{\partial x}$$

Let's define these dimensionless parameters:

$$\theta = \frac{C - C_{b}}{C_{b}}$$

$$\hat{x} = \frac{x}{L}$$

$$\hat{t} = \frac{Dt}{L^{2}}$$

$$Pe = \frac{vL}{D}$$

So, your advection-diffusion equation becomes:

$$\frac{\partial \theta}{\partial \hat{t}} = \frac{\partial^{2} \theta}{\partial \hat{x}^{2}} - Pe \frac{\partial \theta}{\partial \hat{x}}$$

And your boundary conditions:

$$\theta(0,\hat{t}) = 0$$

$$-\frac{\partial \theta}{\partial \hat{x}}|_{\hat{x} = 1} = 0$$

Your initial condition is:

$$\theta(\hat{x},0) = -1$$

By using separation of variable technique you have:

$$\theta(\hat{x},\hat{t}) = X(\hat{x})T(\hat{t})$$

So:

$$XT^{'} = X^{''}T - PeX^{'}T$$

or:

$$\frac{T^{'}}{T} = \frac{X^{''}}{X} - Pe \frac{X^{'}}{X} = -\lambda^{2}$$

For spatial part:

$$X^{''} - Pe X^{'} + \lambda^{2} X = 0$$

The general solution is:

$$X(\hat{x}) = \exp{(\frac{Pe \cdot \hat{x}}{2})} (A\cos{(\eta \hat{x})} + B \sin{(\eta \hat{x})})$$

If put $\hat{x} = 0$:

$$A = 0$$

So:

$$X(\hat{x}) = \exp{(\frac{Pe \cdot \hat{x}}{2})} (B \sin{(\eta \hat{x})})$$

if we put $\hat{x} = 1$:

$$X^{'}(\hat{x}) = \frac{Pe}{2} \exp{(\frac{Pe \cdot \hat{x}}{2})} (B \sin{(\eta \hat{x})}) + \exp{(\frac{Pe \cdot \hat{x}}{2})} B \eta \cos{(\eta \hat{x})}$$

$$X^{'}(\hat{x}) = B \exp{(\frac{Pe \cdot \hat{x}}{2})} (\frac{Pe}{2} \sin{(\eta \hat{x})} + \eta \cos{(\eta \hat{x})})$$

So:

$$X^{'}(1) = B \exp{(\frac{Pe}{2})} (\frac{Pe}{2} \sin{(\eta)} + \eta \cos{(\eta)}) = 0$$

So the eigenvalue equation is:

$$\frac{Pe}{2} \sin{(\eta)} + \eta \cos{(\eta)} = 0$$

or:

$$\tan{(\eta_{n})} = -\frac{2 \eta_{n}}{Pe}$$

We know that $\eta$ and $\lambda$ should be related as:

$$\lambda^{2}_{n} = \frac{Pe^{2}}{4} +\eta^{2}_{n}$$

The temporal solution is:

$$T(\hat{t}) = \exp{(-\lambda^{2}\hat{t})}$$

So:

$$\theta_{n}(\hat{x},\hat{t}) = B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} \exp{\Bigg(-\Big(\frac{Pe^{2}}{4}+\eta_{n}^{2}\Big) \hat{t}\Bigg)}$$

Final solution is:

$$\theta(\hat{x},\hat{t}) = \sum_{n=0}^{\infty} B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} \exp{\Bigg(-\Big(\frac{Pe^{2}}{4}+\eta_{n}^{2}\Big) \hat{t}\Bigg)}$$

$B_{n}$s are extracted by using initial condition:

$$ -1 = \sum_{n=0}^{\infty} B_{n} \exp{(\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})}$$

Or:

$$-\exp{(-\frac{Pe \cdot \hat{x}}{2})} = \sum_{n=0}^{\infty} B_{n} \sin{(\eta_{n}\hat{x})}$$

So:

$$-\exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} = \sum_{n=0}^{\infty} B_{n} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})}$$

Taking integration from $\hat{x} = 0$ to $\hat{x} = 1$:

$$-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \sum_{n=0}^{\infty} B_{n} \int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x}$$

But:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

Where:

$$\gamma_{n} = \int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d \hat{x}$$

So:

$$-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{m}\hat{x})} d\hat{x} = B_{m}\int_{0}^{1} \sin^{2}{(\eta_{m}\hat{x})} d \hat{x}$$

So:

$$B_{n} = \frac{-\int_{0}^{1} \exp{(-\frac{Pe \cdot \hat{x}}{2})} \sin{(\eta_{n}\hat{x})} d\hat{x}}{\int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d \hat{x}}$$

Finally, you can compare your results with different $\Delta x$ or mesh size with this analytical solution.

Update:

Just a note here. I was thinking why we should have:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

We know that:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \frac{\cos{(\eta_{n})} \cos{(\eta_{m})} (\eta_{m} \tan{(\eta_{n})} - \eta_{n} \tan{(\eta_{m})})}{(\eta_{n}^{2}- \eta_{m}^{2})}$$

When $n \neq m$. But we know that: $\tan{(\eta_{n})} = -\frac{2\eta_{n}}{Pe}$ and $\tan{(\eta_{m})} = -\frac{2\eta_{m}}{Pe}$:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = -\frac{2}{Pe}\frac{\cos{(\eta_{n})} \cos{(\eta_{m})} (\eta_{m} \eta_{n} - \eta_{n} \eta_{m})}{(\eta_{n}^{2}- \eta_{m}^{2})} = 0$$

So, we must have:

$$\int_{0}^{1} \sin{(\eta_{n}\hat{x})} \sin{(\eta_{m}\hat{x})} d\hat{x} = \gamma_{n} \delta_{mn}$$

Where:

$$\gamma_{n} = \int_{0}^{1} \sin^{2}{(\eta_{n}\hat{x})} d\hat{x} = \frac{1}{2} - \frac{\sin{(2\eta_{n})}}{4 \eta_{n}}$$

Finally:

$$B_{n} = -\frac{\eta_{n}}{\lambda_{n}^{2}(\frac{1}{2} - \frac{\sin{(2\eta_{n})}}{4 \eta_{n}})}$$

Answer 2

What you're looking for is an a-posteriori error estimate for your mesh study. Normally, these quantitative measures are line/area/volume- and/or time- averaged nodal or element quantities (e.g. avg. temperature of some body integrated over time). In structural analysis it is also customary to use the energy norm as a quantitative measure.

A line/area or volume-averaged quantity is simply to take the average value across a line, area or a volume. The idea is to avoid that your analysis depends on a single node or element, which might be not representative. In your case, since it's a 1D problem, you may want to integrate $c(x)$ from $x_0$ to $x_L$.

Energy norm is a concept used in structural analysis and thus might not be of applicability to you. You can read more about the energy norm here: https://blogs.solidworks.com/tech/2018/05/manage-solidworks-simulation-mesh-refinement-using-an-energy-norm-error-plot.html