4
$\begingroup$

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

enter image description here

I changed delx = 1, the following is observed

enter image description here

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.

enter image description here

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?).

$\endgroup$
  • $\begingroup$ What do the different curves mean on each plot? I couldn't glance it from the code. $\endgroup$ – Javier-Acuna Feb 6 at 18:35
  • 1
    $\begingroup$ The usual way to study the effect of the mesh size is of course to compute the error in the solution. $\endgroup$ – Wolfgang Bangerth Feb 6 at 19:30
  • $\begingroup$ @Javier-Acuna Thanks a lot for the response. The curves denote the transient change in the quantity c observed positions defined in x = 0:delx:xend; % mesh. $\endgroup$ – Natasha Feb 7 at 12:18
  • $\begingroup$ @WolfgangBangerth Thank a lot. I've been looking at absolute and relative errors. Are there instances in which one shouldn't look at absolute/ relative error? For example, if I have 0.5 and 0.2e-25 at point x=1 and time = 0.02, computing(0.5-0.2e-25)/0.2e-25 gives a really high relative error. $\endgroup$ – Natasha Feb 7 at 12:23
  • 1
    $\begingroup$ @Natasha: I would have computed the relative error relative to some "typical size" of the solution -- for example, the mean absolute value of the solution. If you make it relative to the exact value at some point, you're always going to get meaningless numbers anywhere the solution is zero or close to zero. $\endgroup$ – Wolfgang Bangerth Feb 7 at 17:26
2
$\begingroup$

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}})}$$

| cite | improve this answer | |
$\endgroup$
3
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks a lot for the response. From what I understand, let's say we consider the nodal quantities ( at positions x = 1,2,3). At each position, I've two arrays (stores the transient change in c at x positions), one for delx = 0.25 and other for delx = 1. From what has been suggested in the comment above, I can compute absolute and relative errors. And from your suggestions, I understand the time-averaged nodel quantity can also be compared. Could you please explain a bit about the area-averaged quantity and energy norm? $\endgroup$ – Natasha Feb 7 at 12:14
  • 1
    $\begingroup$ I edited my original contribution to answer your questions $\endgroup$ – Ken Grimes Feb 8 at 12:09
  • $\begingroup$ Thank you. Please check my recent edit in the original post $\endgroup$ – Natasha Feb 8 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.