Analytical solution of 1-D transient heat transfer

Can you please help me code the analytical solution equation for 1D transient heart transfer to get the temperature profile? I need to include multiple terms, I suppose it need "for loop" but I'm not familiar how to code loops or sereies.

[Nour]

@Torsten can you help please

[Walter Roberson]

Any time you see something saying that values that you need are tabulated, then you know that they are not looking for analytic solutions.

[Torsten]

• Plot the function f(x) = x*J1(x)/J0(x) - Bi to get an impression where its zeros are located.
• From the plot, choose good starting guesses to determine the first - say ten - zeros of f using "fzero".
• Compute the array C_n.
• Sum the first - say ten - terms of the series for theta_star for different values of r_star and t.
• Make a surface plot of theta_star over r_star and t (or whatever you want to plot).

[Nour]

I'm asked to do so before we move to numerical method to solve this problem.

Below is my code for the one term method, but it is unvalid for my fourier number so I have to include multiple terms to make it work. I would appreciate if you can help me build the code for multiple terms.

clc 
clear all; 
%-----------------------------Input---------------------------------------- 
% Inputs 
r0=0.05;                         % wall thickness(m)
k=0.1;                           % conductivity(W/m-K)
rho=821.42;                      % density(kg/m^3)                   
cp=1399;                         % specific heat capacity(J/kg-K)
T_ini=400;                       % initial temperature(K) 
T_inf=298;                       % external temperature(K) 
h=5;                             % heat transfer coefficient(w/m^2-K) 
t_f=100                      % Final time (sec)
lampda1=1.694                   % Eigen root
C1=1.3787                       % constant
R=0.05;
T_i=400;
T_inf=298;
r=linspace(0,0.05,20)
t=linspace(0,100,t_f)
%-----------------------------Main parameters -----------------------------
% Parameters
Bi=(h*r0)/k                    % Biot number
alpha=k/(rho*cp)               % Thermal diffusivity
F0=(alpha*t)/(r0*r0)           % Fourier number
r_div=20;
epsilon=1.694;
C_1=1.3787;

T=zeros(1,r_div);
u=zeros(1,r_div);
t=0;
%-----------------------------Temperature calculation ---------------------
for j=1:5
    for i=1:0.05
        u(i)=epsilon*r(i) /R;   
        T(i)=T_inf+(T_i-T_inf)*C_1*exp(-epsilon^2*alpha*t/R^2)*besselj(0,u(i));
    end
end

Answer 1

[Nipun]

Hi Nour,

I understand that you are trying to code the analytical solution equation for 1D transient heart transfer to get the temperature in profile in MATLAB.

Based on the provided system of equations, I propose the following steps to approach the solution:

Start with getting the roots of the transcendental equations. Plot the function and observe the first n zeroes of the equation

syms x
f(x) = x*J1(x)/J0(x) - Bi
fplot(f(x))

Say the interval of the roots is [a,b] then, use "vpasolve" to obtain the roots in that interval

roots = vpasolve(f(x)==0, x, [a,b])

Based on the array of roots, compute Cn using the relation:

C(i) = (2/x)*(J1(x)/(J1(x)^2 + J0(x)^2)) % where x = roots(i)
% You can also represent it in vector format

Compute the temperature profile using sum

temp_profile = sum(C.* exp(-roots.^2 * F0)* J0(roots*r) ) % vector representation

On increasing the size of roots (by appending more zeroes to the array and increasing the interval [a,b]), the accuracy of "temp_profile" will increase.

Link to documentation:

1. symbolic representation : Equation Solving - MATLAB & Simulink - MathWorks India
2. fplot in MATLAB : Plot expression or function - MATLAB fplot - MathWorks India
3. vpasolve : Solve symbolic equations numerically - MATLAB vpasolve - MathWorks India

Hope this helps.

Regards,

Nipun