Non Linear Heat Conduction - Crank Nicolson

I am trying to model the temperature development of a metal during friction welding (initially just in 1D). I would really appreciate it if someone could give me an idea of where I may be going wrong.

So, I have written/modified a Crank-Nicolson finite difference scheme to tackle the problem and have been comparing my results with a model that I have produced using COMSOL . When I use constant material properties the two models produce results that are within 0.1% of each other. When I use variable material properties the MATLAB model produces results that are ~10% lower than the COMSOL model (100C lower over a 1000C temperature range). I'm inclined to trust COMSOL and am concerned about my implementation of the variable properties and the validity of the iterative method I use.

How I worked through the problem (the code is at the bottom):

-------------------------

1. Set up my mesh spacing, initial temperature profile and my heat flux

2. Produce my sparse tridiagonal matrix and the coefficients on the RHS of the equation (i.e. the A and d in Ax=d)

3. Solve that problem using the backslash operator

I now have a set of temperatures at a time of t+1 but due to the non linearity of the problem I believe that these are incorrect as the coefficients in the tridiagonal matrix should be based upon these temperatures(and they're not!). What I therefore do is:

4. Use this new set of temperatures to create a new tridiagonal matrix and then solve the problem once again

5. I re-input these new values and repeat the process until convergence.

So, any thoughts?

Kind Regards

Chris

%Crank Nicolson scheme for 1D friction welding

clear

%------Mesh Description and Initial Conditions------%

nt=10; tmax=10; Dt=tmax/(nt); tt=0:Dt:tmax;     % Time spacing (s)
Dz=0.02;  L=20; z=(0:Dz:L); nz=length(z);       % Grid spacing (mm)
T=(20+(z)*40).';                                % Initial temperature profile
q=2*ones(length(tt),1);                         % Heat flux

Tstore=zeros(length(z),length(tt)-1);           %Set up matrix to store all temperature values
Tstore(:,1)=T;

% --- Loop over time steps

for k=1:nt-1

% --- Temperature varying material properties

kk =(T/120 + 9.8)/1000;
rho= 8.2;
Cp= (T/3000 +0.4)/1000;
alfa=kk./(rho.*Cp);

% --- Tridiagonal matrix

b = [(-alfa(2:end)/(2*Dz^2));0];          % Super diagonal coefficients
c = [0;(-alfa(1:end-1)/(2*Dz^2))];        % Subdiagonal coefficients
a = ((1/Dt) +(alfa/Dz^2));                % Main Diagonal
at = (1/Dt) -(alfa/Dz^2);                 % Coefficient of u_i^k on RHS

b(1) = -(alfa(2)/Dz^2);                   % Coefficients of boundary nodes
c(1)=0;
b(end)=0;
c(end) = (-alfa(end-1)/(Dz^2));

A=(spdiags(a,0,length(a),length(a)));
A=(spdiags([0;b],1,A));                   % Super diagonal with spacer added to correctly position vector
A=(spdiags(c(2:end),-1,A));               % Sub-diagonal with values taken from c2 onwards (as c1 is outside matrix)

% --- Update RHS for all equations, including those on boundary

d = - [0;c(2:end).*T(1:end-1)] ...
+ at(1:end).*T(1:end) ...
- [b(1:end-1).*T(2:end);0];

d(1)=d(1)+ (alfa(1)./(kk(1).*Dz)).*(q(k+1)+q(k)); %Include heat flux

Titt=(A\d);
Tit=ones(length(T),10);
Tit(:,1)=Titt;

% The tridiagonal matrix was formed using the material properties based
% upon the initial temperature - this is incorrect as the tridiagonal matrix 
% is a function of the final (or timstep t+1) temperature. 

% The (incorrect) temperature profile calculated above can be used to update 
% the material properties and the same calculations can be performed to 
% iterate to the 'correct' new (t+1) temperature.

for it=2:10 %In practice only ~5 iterations are required for convergence

kkit =(Titt/120 + 9.8)/1000;
rhoit= 8.2;
Cpit= (Titt/3000 +0.4)/1000;
alfait=kkit./(rhoit.*Cpit);

bit = [(-alfait(2:end)/(2*Dz^2));0];
cit = [0;(-alfait(1:end-1)/(2*Dz^2))];
ait = ((1/Dt) +(alfait/Dz^2)); 
atit = (1/Dt) -(alfait/Dz^2); 

bit(1) = -(alfait(2)/Dz^2); 
cit(1)=0;
bit(end)=0;
cit(end) = (-alfait(end-1)/(Dz^2));

Ait=(spdiags(ait,0,length(ait),length(ait)));
Ait=(spdiags([0;bit],1,Ait));                             
Ait=(spdiags(cit(2:end),-1,Ait));                         

d(1)=at(1)*T(1)-b(1)*T(2)+ (alfa(1)./(kk(1).*Dz)).*q(k)...
   +(alfait(1)./(kkit(1).*Dz)).*(q(k+1));

Titt=(Ait\d);
Tit(:,it)=Titt;

end

T=Titt;
Tstore(:,k+1)=T;

end

[Chris Simpson]

I thought I'd update this with an answer to my own question!

The problem was the finite difference discretisations that I was using. They needed to be more complex to accommodate the variable properties -using an average heat diffusivity is just incorrect! The below image should hopefully show you how I dealt with the variable heat capacity and thermal conductivity. I used this as the basis of a FTCS explicit scheme, but I'm sure the principles hold for implicit schemes:

https://lh4.googleusercontent.com/-UyhK-80w5Es/UJ1jhBqmCNI/AAAAAAAAZt8/VnmZO_TQ0Ik/s691/Nonlinear%2520Heat%2520Equation.gif

If that doesn't work try this web page (not mine!):

http://jupiter.ethz.ch/~kausb/Teaching/CompGeoFinDiff/Finite_Differerence_NonLinear.pdf

This only looks at a variable thermal conductivity type problem. It should be noted that problems with variable heat capacity or density would require an iterative process (linearising the non linear problem) as they rely on the temperature at the time-step t+1 as well as time-step t:

Hope this helps someone!

Chris

[mtl]

Hi Guys,

do you have any idea about how to build the tridiagonal matrix for the head conduction problem with two different materials?

how to get the k(t, m+1/2), C(t+1/2, m), density (t+1/2, m)

Many thanks.