1D heat equation crank nicholson

Hey, I'm trying to solve a 1d heat equation with the crank nicholson method. A one-step algortihm for the semidiscrete heat equation (generalized trapezoidal method).

Solving this problem:

rho * c * T' = k * T,xx + Q'_v => ( Me * T'e + Ke * Te = Fe )

with:

M * v_n+1 + K * d_n+1 = F_n+1

d_n+1 = d_n + dt * v_n+alpha

v_n+alpha = (1 - alpha) * v_n + alpha * v_n+1

alpha is 0.5 for crank nicholson

Geometrie: it's a rod

Temperature a t=0 0 everywhere

So, i do not really know how to implement a heat flux (W/m^2) on one boundary side.

Please help :)

clear all;
close all;
clc

%%%--- Mesh generation ---%%%
le=0.01;                      % Distance between nodes [m]
units=0.1;                   % Length of rod
p=[0:le:units]';             % Nr of nodes
x = [0:le:units];            % Spatial discretisation
nr_nodes = length(p);       % Nr of nodes

  for   i = 2:nr_nodes        % Creates the element matrix
      e(i-1,:) = [i-1 i];
  end

%%%--- Creating space ---%%%
nr_elements =length(e);     % Nr of elements

K=zeros(nr_nodes,nr_nodes); % Empty stiffness matrix
M=zeros(nr_nodes,nr_nodes); % Empty mass matrix
F=zeros(nr_nodes,1);        % Empty forcing vector
Q=zeros(nr_nodes,1);        % Initial forcing vector

%%%--- Parameters ---%%%
k=50;                      % Heat conductivity
pc=7800*450;               % Heat capacity
A=0.01*0.005;              % Element cross section
q=1e4;                   % Heat flux
Q(1)=q;

%%%--- Assembling matrices ---%%%
for j = 1:nr_elements                        % Loop over all elements
    
        nodes=e(j,:);                        % Identify nodes
        le=p(j+1)-p(j);                       % Calculate distance between nodes
        
        Ke = ((A*k)/le)*[1 -1;-1 1];          % Stiffness entry
        Me = ((pc*A*le)/6)*[2 1;1 2];         % Mass entry       
        Fe = ((A*le)/4)*[2 1;1 2]*Q(nodes,1); % Forcing entry
                            
        K(nodes,nodes)=K(nodes,nodes)+Ke;    % Position stiffness entry at identified nodes
        M(nodes,nodes)=M(nodes,nodes)+Me;    % Position mass entry at identified nodes      
        F(nodes,1)=F(nodes,1)+Fe;            % Position forcing entry at identified nodes
end
%F(1)=q;

%%%--- Initial conditions ---%%%
Told=zeros(nr_nodes,1);     % Initial temperature, 0 everywhere
%F(1)=1;                     % With unit heat production at node one

%%%--- Time parameters ---%%%
dt = 0.01;                  % Size of time step
t=10;                        % Time

alpha   =   0.5;            % Crank Nickolson
%%%--- Preparing matrices ---%%%

Tnew    =   Told;
Fold    =   F;
Fnew    =   Fold;
d0      =   zeros(nr_nodes,1);

    
nu0     =   inv(M)*(F-K*d0);

temp    =   zeros(nr_nodes,t/dt);


for i =1:t/dt
    
    thermometer(i)=Tnew(1);                    % Vector to store temperature data at node 1

    temp(:,i)   =   Tnew;   
    
    if i==1
    dguess  =   Told+(1-alpha)*dt*nu0;
    else
    dguess  =   Told+(1-alpha)*dt*nu;
    end
    nu  =   inv(M+alpha*dt*K)*(F-K*dguess);

    Tnew   =   dguess+alpha*dt*nu; %d

% plot(x,Tnew);
Told   =    Tnew;
% plot(x,Tnew)
% pause(0.1)
end

% figure(1)
plot(x,temp(:,end))
time=[0:0.01:10-0.01];
plot(time,temp(1,:));

[Tatiane Cristine Antonio]

Hi, did you figure how to do this? I'm trying to do the same! Thanks!