1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Jul 8, 2017 #1
    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?

    Attached Files:

  2. jcsd
  3. Jul 8, 2017 #2
    Sure. See Conduction of Heat in Solids by Carslaw and Jaeger.
  4. Jul 8, 2017 #3
    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.
    Code (Text):
    clear all; close all; clc;
    deltat = 0.5;

    v = 40;
    pvec = 1:50;  % total number of discrete values are 50
    u = x/sqrt(v*deltat);  % field variable
    for i=1:length(pvec)
        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
    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.

    Attached Files:

  5. Jul 9, 2017 #4
    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
  6. Jul 9, 2017 #5


    User Avatar
    Science Advisor
    Gold Member

    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.wordpr...solution-of-pdes-part-3-2d-diffusion-problem/ .
  7. Jul 9, 2017 #6
    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted