A How to validate a code written for solution of 1D diffusion?

[Atr cheema]

Consider the conceptual model presented in the attached image, of heat conduction in a bar.

There is a heat source at left side and heat is observed at point Ho after a distance L from the source. If we consider only heat transfer through conduction then this problem can be modeled by diffusion/heat conduction equation.
$$ \frac{\partial T}{\partial t} =D \frac{\partial T}{\partial x^2} $$
If I have developed a new semi-analytical or numerical solution (and eventually code) of this problem, then how can I validate my code?. For simplicity, suppose the source to be sinusoidal. Is there a standard solution for such an ideal case with sine input source, so that I can compare my code and then apply it to real field conditions?

 

[Chestermiller]

Consider the conceptual model presented in the attached image, of heat conduction in a bar.

There is a heat source at left side and heat is observed at point Ho after a distance L from the source. If we consider only heat transfer through conduction then this problem can be modeled by diffusion/heat conduction equation.
$$ \frac{\partial T}{\partial t} =D \frac{\partial T}{\partial x^2} $$
If I have developed a new semi-analytical or numerical solution (and eventually code) of this problem, then how can I validate my code?. For simplicity, suppose the source to be sinusoidal. Is there a standard solution for such an ideal case with sine input source, so that I can compare my code and then apply it to real field conditions?

Sure. See Conduction of Heat in Solids by Carslaw and Jaeger.

 

[Atr cheema]

Sure. See Conduction of Heat in Solids by Carslaw and Jaeger.

I have gone through some chapters of that book. He provides semi-analytical solutions in the form of equations for which we need to write code for example for example in case of semi infinite rod the solution he provides is
$$ \Delta h_m = \sum\limits_{m=1}^p \Delta H_m \{ erfc \frac{u}{2 \sqrt[]{p-m}} \} $$ p is total number of time steps and u is field variable.
I am unable to write a proper code for this equation. I have written following MATLAB code for this but it is giving me erratic results.

clear all; close all; clc;
deltat = 0.5;
deltaHvecc=sin(linspace(0.3,2.6,50));

deltaHvec=deltaHvecc-deltaHvecc(1);
x=45;
v = 40;
pvec = 1:50;  % total number of discrete values are 50
%suma=linspace(0,0,length(pvec));
u = x/sqrt(v*deltat);  % field variable
for i=1:length(pvec)
    p=pvec(i);
    deltaH=deltaHvec(i);
    suma = 0;
    for m=1:i
        erfcArg =  (u/(2*sqrt(p-m)));
      erfcfunc=erfc(erfcArg);  % calculating error function
      tek = deltaH*erfcfunc;
      suma = suma + tek;          % summation in the equation
    end
summa(i)=suma;
end
plot(summa); hold on; plot(deltaHvecc,'k')

The result of this code is given in the image I have attached with this answer, which is wrong as calculated can not be more than source.
So I can't compare the solution that I have derived with that of Carslaw and Jaeger.
My problem will be solved if either I am able to correct my code or get sample standard solution for the problem posed in original question.

 

[Chestermiller]

This is not the solution I had in mind. The solution I had in mind was for the dynamic steady state at long times, in which the temperature approaches $$T=T_0+A(x)\sin{\omega t}+B(x)\cos{\omega t}$$ This solution is given in Transport Phenomena by Bird, Stewart, and Lightfoot, Chapter 4, Example 4.1-3

 

[hilbert2]

The factor $\frac{\Delta t}{(\Delta x )^2}$, where $\Delta t$ is the timestep and $\Delta x$ the position step in the discretization, needs to be small enough or otherwise the solution will become unstable (small wavelength Fourier components of the solution will "explode" in amplitude). The code can be tested against the self similar solution

$C(x,0) = Ae^{-B(x-x_0 )^2}$
$C(x,t) = A(t)e^{-B(t) (x-x_0 )^2}$,

called "self-similar" because it's Gaussian for all $t \in \mathbb{R}$ (the functions A(t) and B(t) can be found from several sources), or simply the linear constant heat flux solution where the temperature/concentration are forced to stay constant at the endpoints of the domain.

I haven't done this for a sine/cosine time dependent source term (in the wave equation case that kind of term would excite oscillation modes that are in resonance with the frequency of the source term), but I have some related codes here in my blog:

https://physicscomputingblog.wordpress.com/2017/02/11/numerical-solution-of-pdes-part-1-explicit-method/
https://physicscomputingblog.wordpress.com/2017/02/11/numerical-solution-of-pdes-part-2-implicit-method/
https://physicscomputingblog.wordpress.com/2017/02/20/numerical-solution-of-pdes-part-3-2d-diffusion-problem/ .

 

[Chestermiller]

The solution in Bird et al is: $$T=T_0+Ae^{-\sqrt{\omega/2D}x}\cos{(\omega t-\sqrt{\omega/2D}x)}$$where A is the amplitude of the temperature variation at x = 0.