Modeling analytical solution of 1D heat equation

[Atr cheema]

I am trying to write code for analytical solution of 1D heat conduction equation in semi-infinite rod. The analytical solution is given by Carslaw and Jaeger 1959 (p305) as
$$
h(x,t) = \Delta H .erfc( \frac{x}{2 \sqrt[]{vt} } )
$$
where x is distance, v is diffusivity (material property) and t is time. I have written following code in MATLAB for to find heat distribution at 40 m from the input with time. The input is a sine wave at x=0.

clear; close all; clc;
x = 40;       % distance from input to observation point
v = 100;      % diffusivity
p = 300;        % number of steps
Ti = 0;       % Initial head
t =  1:p;    % time steps at which calculation is done, lenth(t) should be
equal to length(To)
delH = sin(linspace(0.1,60,p));

h = zeros(1,p);

for m = 1:p
    func = x/(2*sqrt(v*t(m)));
    h(m) = delH(m) * erfc(func);

end
plot(delH)
hold on;
plot(h)
legend('input','output')
xlabel('time(sec)')
ylabel('temperature')

I am not understanding why the output sinewave keeps its amplitude increasing with time? Should it just mimic the input sine wave with reduced amplitude. The yellow line in

this image shows result from finite difference and I expect result from Carslaw to mimic finite difference result.

 

[mjc123]

I am not understanding why the output sinewave keeps its amplitude increasing with time?

Look at your equation. Is there a clue there? ΔH is your input sine wave; you are multiplying it by a factor that increases with time.
I am wondering whether the equation is the right one for the scenario - is it perhaps for constant ΔH? Someone correct me if I'm wrong. I remember doing something like this some years ago, and I think the equation was of the form
ΔT(x,t) = ΔT₀e^i(ωt-βx)e^-αx
This (or the real part of it) gives at x a sine wave of reduced but constant amplitude, and phase shifted relative to x=0 (which your solution isn't).

 

[Atr cheema]

Look at your equation. Is there a clue there? ΔH is your input sine wave; you are multiplying it by a factor that increases with time.
I am wondering whether the equation is the right one for the scenario - is it perhaps for constant ΔH? Someone correct me if I'm wrong. I remember doing something like this some years ago, and I think the equation was of the form
ΔT(x,t) = ΔT₀e^i(ωt-βx)e^-αx
This (or the real part of it) gives at x a sine wave of reduced but constant amplitude, and phase shifted relative to x=0 (which your solution isn't).

Thank you for your response. You are exactly right that I should have a sine wave at x with constant but reduced amplitude and with a relative phas shift. I am using equation 15 on page 357 of Jan Taler's book. Page 15 is here. Can you please tell me if equation 15 is for this scenario or not?

 

[mjc123]

It's a bit small, and I can't scroll up and down, but it looks to me as if ΔT(x=0) is supposed to be constant. From your table of contents, this section is dealing with "Formula Derivation for Temperature Distribution in a Half-Space with a Step Increase in Surface Temperature", i.e. T(x=0) increases at time t=0 from T₀ to T₀ + ΔT₀ (constant). It looks as if you want p.366 "Formula Derivation for a Quasi-Steady State Temperature Field in a Half-Space when Surface Temperature Changes Periodically".

 

[Chestermiller]

You're using the equation incorrectly. In your system, $\Delta H$ the temperature change at the boundary is not constant with time. You need to express the solution to your problem as a convolution integral with appropriate time delays in the erfc term. What would the solution look like if expressed as an integral in this way?

 

[Atr cheema]

Thank you again for your insightful comments.

It's a bit small, and I can't scroll up and down, but it looks to me as if ΔT(x=0) is supposed to be constant. From your table of contents, this section is dealing with "Formula Derivation for Temperature Distribution in a Half-Space with a Step Increase in Surface Temperature", i.e. T(x=0) increases at time t=0 from T₀ to T₀ + ΔT₀ (constant). It looks as if you want p.366 "Formula Derivation for a Quasi-Steady State Temperature Field in a Half-Space when Surface Temperature Changes Periodically".

The problem is that the input function in my real case is not completely sinosoidal but p.366 refers to cases where input is completely periodic. Here is complete Jan Taler's book. Actually in this paper, the author has shown how he used the solution of Carslaw and Jaeger to derive to the final equation. He uses equation equation of Carsla (2a in paper) to come to the conclusion (3a in paper). Equation 2a is same as I wrote above and as given by Carslaw on page 305 (equation 5).

 

[Atr cheema]

You're using the equation incorrectly. In your system, $\Delta H$ the temperature change at the boundary is not constant with time. You need to express the solution to your problem as a convolution integral with appropriate time delays in the erfc term. What would the solution look like if expressed as an integral in this way?

Isn't $\Delta H$ changing with every step in my code? I am sorry I am not very good at math and its been quite a time I am looking at this problem. Can you tell what changes I have to make in my code?

 

[Chestermiller]

The equation should be $$h(x,t)=\int_0^t{\frac{d\Delta H(t-\tau)}{d\tau}\left[erfc{\frac{x}{2\sqrt{\nu \tau}}}\right]d\tau}$$