[Alvin1990]

I am trying to solve 1D convection diffusion equation using finite difference (Crank Nicholson Method):









0.0002 to 0.001 incrementing by 0.0001

The question is to find the time it requires and compare that time with alpha values.

For the x-direction, central difference and for time derivative, forward difference was used. For the right boundary, backward difference is used.

I tried my code in python:

import numpy as np
import matplotlib.pyplot as plt

j = 7
u = 0.002*(1+j)
L = 2
ndx = 1001
dx = L/(ndx-1)

alpha = np.arange(0.0002, 0.001+0.0001, 0.0001)
dt = min(dx**2/(2*min(alpha))*.8, 0.8*dx/u)
ndt = int(320/dt)+1

x = np.linspace(0, 2, ndx)
Temperature_of_node_to_check = int(1.5/dx)

t_1_5 = np.zeros(alpha.size)

courant_number = u*dt/dx
if courant_number < 1:
    print("C < 1" + "ok")
peclet_number = u*L/alpha
if any(peclet_number <= 2):
    print("P <= 2" + "ok")
cell_peclet_number = u*dx/(2*alpha)
if any(cell_peclet_number < 2):
    print("C.P < 2" + "ok")
s = alpha*dt/dx**2
if any(s <= 0.5):
    print("s <= 0.5" + "ok")

courant_number = (u*dt/dx)*np.ones(alpha.size)
diffusivity = alpha*dt/dx**2


a = -(courant_number/2.0+diffusivity)
b = 2*(1+diffusivity)
c = (courant_number/2.0-diffusivity)

T = np.zeros(shape=(ndt, ndx))
T[0,:] = 10+j
T[:,0] = 30

D = np.zeros(ndx-2)*1.0

def thomas(a, b, c, D):
    N = D.size
    C = np.ones(N-1)*c
    B = np.ones(N)*b
    A = np.ones(N-1)*a
    
    A[N-2] = a-c/3.0
    B[N-1] = b + (4/3.)*c
    C[0] = C[0]/B[0]
    D[0] = D[0]/B[0]
    for i in range(1, N-1):
        C[i] = C[i]/(B[i-1]-A[i-1]*C[i-1])
        D[i] = (D[i]-A[i-1]*D[i-1])/(B[i]-A[i-1]*C[i-1])
        
    D[N-1] = (D[N-1]-A[N-2]*D[N-2])/(B[N-1]-A[N-2]*C[N-2])
    result = np.zeros(N)
    result[N-1] = D[N-1]
    for i in range(N-2, -1, -1):
        result[i] = D[i] - C[i]*result[i+1]
    return result


k = 0
for ai, bi, ci, alphai, coi, diff in zip(a,b,c,alpha, courant_number, diffusivity):
    for i in range(ndt-1):
        for j in range(1, ndx-1):
            D[j-1] = (coi/2.0+diff)*T[i,j-1] + 2*(1-diff)*T[i,j] + (-coi/2+diff)*T[i,j+1]
        D[0] = D[0] - ai*T[i,0]
        T[i+1,1:ndx-1] = thomas(ai,bi,ci,D)
        T[i+1,ndx-1] = (4/3.)*T[i,ndx-2] - (1/3.)*T[i,ndx-1]
        if T[i+1, Temperature_of_node_to_check] >= 25.:
            t_1_5[k] = (i+1)*dt
            print((i+1)*dt)
            break
    fig, ax = plt.subplots(1, 1, figsize = (24,24))
    label = r'$\alpha$' + '={:.5f} '.format(alphai) + r'$m^{2}/s$' + r' t = {:5f}'.format((i+1)*dt)
    ax.plot(x, T[i+1,:], "--g", label = label)
    ax.set_xlabel("x", fontsize = 30)
    ax.set_ylabel("T", fontsize = 30)
    ax.legend(fontsize = 30)     
    plt.plot(x, T[i+1,:], "--g")
    plt.yticks(np.linspace(16,32, 9))
    plt.plot([0,1.5],[25,25],'k')
    plt.plot([1.5,1.5],[25,17],"k")
    filename = "alpha={:.5f}".format(alphai) + ".jpg"
    plt.savefig(filename)
    k = k + 1
    plt.close()

plt.plot(alpha, t_1_5, "--r")
title = "Time Taken to rech " + r"$25^{o} C$" + " vs " + "diffusivity"
plt.title(title, fontsize = 10)
plt.xlabel(r"$\alpha$", fontsize = 20)
plt.ylabel(r"$T_{1.5}$", fontsize = 20)
plt.yticks(fontsize = 10)
plt.xticks(fontsize = 10)

Here the graph generated from t vs [math]\alpha[math]. As the thermal diffusivity increases from 0.0002 to 0.001, heat transfer decreases as it requires more time.



After 0.002, the heat transfer rate increases which is what I was expecting.



What happens in lower values of alpha that heat transfer decreases? Is it due some numerical error or some other reason?

[Alvin1990]

the images are uploaded in google drive

https://drive.google.com/file/d/1OMU...usp=drive_link

https://drive.google.com/file/d/1KYV...usp=drive_link

[LuckyTran]

Central differencing is not appropriate for the advective term and is terribly inaccurate when the Peclet number is not less than 2.


You can repeat your calculations at different grid sizes and immediately see the impact. And since Pe<2 are diffusion dominated problems, this is why people don't use central differencing for the advective term for convection problems. To be clear, central differencing the diffusion term is 100% okay.

[Alvin1990]

Originally Posted by LuckyTran

Central differencing is not appropriate for the advective term and is terribly inaccurate when the Peclet number is not less than 2.


You can repeat your calculations at different grid sizes and immediately see the impact. And since Pe<2 are diffusion dominated problems, this is why people don't use central differencing for the advective term for convection problems. To be clear, central differencing the diffusion term is 100% okay.

Thanks you for the comment. If I reduce the grid size there is a significant difference. Will upwind scheme circumvent the issue?

[LuckyTran]

The alternative is, you stick with central differencing but using a super dense grid to drive down the Peclet number.

A pure (1st order) upwind scheme has the same issue but the other way. Upwind isn't great for Pe<2. However, since it's very easy to create a coarser grid and jack up the Peclet number, it's favorable and is more optimal. I mention this scheme mostly for educational purposes. A 2nd order upwind scheme or Crank-Nicolson scheme is straightforward to code and works for most problems but the error analysis doesn't give a simple clear rule like Pe>2 versus Pe<2

[FMDenaro]

What you are using is the so called FTCS scheme. It works accurately provided that the cell Peclet number is <2 as previously stated.
Furthermore the stability constraint for such a low cell Peclet value requires a small time step.

I strongly suggest to check first your code using the analytical solution for the steady case with Dirichlet conditions. Then, if it works correctly, check the solution with Neumann BC.

[Alvin1990]

I edited my code to take account of dx and dt, so courant number < 1 and alpha*dt/dx^2 <= 0.5. The cell peclet number stays less than < 2. As for higher values of alpha the results start to diverge. What other methods can be used to get consistent and correct result? I have used implicit method(unconditionally stable) to solve it but difficulty understanding whether this method provides accurate result or not?

[FMDenaro]

Originally Posted by Alvin1990

I edited my code to take account of dx and dt, so courant number < 1 and alpha*dt/dx^2 <= 0.5. The cell peclet number stays less than < 2. As for higher values of alpha the results start to diverge. What other methods can be used to get consistent and correct result? I have used implicit method(unconditionally stable) to solve it but difficulty understanding whether this method provides accurate result or not?

The issue is in the physics. If you prescribe an entering heat flux on the left and adiabatic condition in the right, the temperature can only increase up to infinity

[Alvin1990]

There is no heat flux on the left end. It will have constant 30 deg temperature. Right side is insulated. Maximum temperature possible is 30 deg celcius.

[FMDenaro]

Originally Posted by Alvin1990

There is no heat flux on the left end. It will have constant 30 deg temperature. Right side is insulated. Maximum temperature possible is 30 deg celcius.

You are not award that while prescribing the temperature you will get a heat flux? Compute it during the computation.

[Alvin1990]

There will be heat flux as node 1 temperature is 30 deg and node 2 temperature 17 deg. But this difference/heat flux should go down due to advection and diffusion and when all nodes will reach 30 deg, heat transfer will stop.

[FMDenaro]

Originally Posted by Alvin1990

There will be heat flux as node 1 temperature is 30 deg and node 2 temperature 17 deg. But this difference/heat flux should go down due to advection and diffusion and when all nodes will reach 30 deg, heat transfer will stop.

The steady state is the constant but you need to be sure the the transient has no oscillations. As I wrote, check first the analytical solution.