How to solve the ODE based 1D heat equation with arbitrary optical heat source term?

hello, I want to solve the ODE-based 1D heat equation when the arbitrary light source is injected into material system.

previously, I have done this calculation using Euler method, however, there is no information whether this kind of calculation method is valid or not.

I attach the my matlab code in here.

%%%%%%%%%%%%%%% %%%%%%%%%%%%%%

% Heat Equation (Differential Equation --> Using Recurrence Relation)

clear,clc;

% Fixed value

q = 1.602e-19; % Electric charge, coulomb

h = 6.62607004081e-34; % Plank constant, [J*sec]

hbar = h/(2*pi); % Plank constant, [J*sec]

Vf = 1e6; % Fermi velocity, [m/s]

Lcool = 1e-6 ; % Electron cooling length, 1 um

Kb = 1.3806488e-23; % Boltzman constant, Fixed value, [ev/K] % Kb = 1.3806488e-23 [J/K] (8.6173e-5) [eV/K]

Tl = 300 ; % Lattice temperature, 300 K (Room temperature, Fixed value)

L0 = (pi^2)*(Kb^2)/(3*q^2); % 2.4435e-8 [V/K]^2 (Fixed Value)

sig0 = 0.23e-3; % Minimum conductivity, 0.23 [mS] (Fixed value) // sig0 = 5*(q^2/h) = 0.19366 [mS];

mob = 0.2; % [m^2/(V.sec)]

Ef = 0.3; % Intrinsic Fermi Level of graphene on SiO2 at no bias [eV] --> Variable

Vg = 0.5; % Vg = Vg1 = Vg2, Gate voltage

Vdirac = 0.7; % Dirac voltage : 0.7 V

dg = 10e-9; % h-BN thickness : 13 nm

e0 = 8.854e-12; % Vacuum permittivity [F/m]

eg = 6.9; % h-BN's relative permittivity

W = 80e-9; % Slot width

L = 30e-6; % Length of graphene

Pfiber = 0.8e-3; % Incident optical power in the fiber, 0.8 [mW]

CE = 0.4 ; % Coupling efficiency of grating coupler

alpha = 0.28; % Absorption at 30 um propagation (l = 30 um) alpha = 0.28;

% Variables

Cg = eg*e0/dg ; % Gate capacitance

sig_DC = sqrt((sig0)^2 + (mob*Cg*(Vg-Vdirac))^2) % DC conductivity of graphene

Kel = L0*Tl*sig_DC ; % Electronic thermal conductivity [W/K]

Pin = CE*Pfiber; % 0.32 mW

Pload = load('PowerDistributionSlot_max0.02_Uniform.txt');

xpos = Pload(:,1);

% Determination of X step size - Version 1

xpos = -2e-6 : 1e-9 : 2e-6;

dx = 1e-9; % [m]

Nx = length(xpos); % Number of X grid

x = (-(Nx-1)/2 : (Nx-1)/2)*dx; % Define x from -2 um to 2 um with 1 nm mesh

Px(1:Nx) = 0;

Px(1960:2040)=1;

P = alpha*Pin*(Px/sum(Px))/(W*L); % T(2) = 300.00202; Vg = 0.5V

% Initialization

T(1:Nx) = 0 ; % Function T

T(1) = 300; % Initial Conditions for T(0)

T(2) = 300.00202; % Initial Conditions for T'(0) ~ (T(2)-T(1))/dx

% Recurrence process to solve heat equation - Version 1

for n = 2:Nx-1;

T(n+1) = (2+((dx^2)/(Lcool^2)))*T(n) - T(n-1) - ((Tl/(Lcool^2))+(P(n)/(Kel)))*(dx^2);

T(Nx) = 300;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In upper method, Only initial value at the most left position was set to 300 K. We have to set the 300 K at the most right position samely, however, I cannot set to 300K using this previous Euler method.

So the new calculation method should be used, considering Boundary Value Problem at the most left and right position set to 300 K.

I'm considering the solvers such as ode45, bvp4c.

I don't know how to include optical heat source term which I uploaded into these solvers.

please help me.

Answer 1

[Alan Stevens]

Here's a possible solution. I've assumed the temperatures are symmetric about the midpoint. I've extracted data from your file so that there is a unique P value for each position (which isn't the case for your data file) - see attached xPdata file. I've assumed the units are all consistent.

% Heat Equation (Differential Equation --> Using Recurrence Relation)
clear,clc;
% Fixed value  
q = 1.602e-19;                                   % Electric charge, coulomb
h = 6.62607004081e-34;                           % Plank constant, [J*sec]
hbar = h/(2*pi);                                 % Plank constant, [J*sec]
Vf = 1e6;                                        % Fermi velocity, [m/s]
Lcool = 1e-6 ;                                   % Electron cooling length, 1 um
Kb = 1.3806488e-23;                              % Boltzman constant, Fixed value, [ev/K] % Kb = 1.3806488e-23 [J/K] (8.6173e-5) [eV/K]
Tl = 300 ;                                       % Lattice temperature, 300 K (Room temperature, Fixed value)
L0 = (pi^2)*(Kb^2)/(3*q^2);                      % 2.4435e-8 [V/K]^2 (Fixed Value)
sig0 = 0.23e-3;                                  % Minimum conductivity, 0.23 [mS] (Fixed value) // sig0 = 5*(q^2/h) = 0.19366 [mS];
mob = 0.2;                                       % [m^2/(V.sec)]
Ef = 0.3;                                        % Intrinsic Fermi Level of graphene on SiO2 at no bias [eV] --> Variable  
Vg = 0.5;                                        % Vg = Vg1 = Vg2, Gate voltage
Vdirac = 0.7;                                    % Dirac voltage : 0.7 V
dg = 10e-9;                                      % h-BN thickness : 13 nm
e0 = 8.854e-12;                                  % Vacuum permittivity [F/m]
eg = 6.9;                                        % h-BN's relative permittivity
W = 80e-9;                                       % Slot width
L = 30e-6;                                       % Length of graphene
Pfiber = 0.8e-3;                                 % Incident optical power in the fiber, 0.8 [mW]
CE = 0.4 ;                                       % Coupling efficiency of grating coupler
alpha = 0.28;                                    % Absorption at 30 um propagation (l = 30 um) alpha = 0.28;
% Variables
Cg = eg*e0/dg ;                                  % Gate capacitance
sig_DC = sqrt((sig0)^2 + (mob*Cg*(Vg-Vdirac))^2); % DC conductivity of graphene
Kel = L0*Tl*sig_DC ;                             % Electronic thermal conductivity [W/K]
Pin = CE*Pfiber;                                 % 0.32 mW
% Data above this point unchanged from the original
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% xPdata contains data from PowerDistribution_max0.02_Uniform.txt but
% has averaged the P values where they were repeated for a given x value.
load('xPdata');                            % xPdata has no repetition
x = xP(:,1); P = xP(:,2);                  % Extract x and P values 
Nx = numel(x);                             % Number of X grid
nmid = (Nx+1)/2;                           % midpoint number

% Initialization
T(1:Nx) = 0 ;                              % Function T
T(1) = 300;                                % Initial Conditions for T(0)
T(2) = 300.00202;                          % Initial Conditions for T'(0) ~ (T(2)-T(1))/dx

% Recurrence process to solve heat equation 
% P is symmetric about midpoint, boundary conditions are symmetric,
% so assume temperatures are symmetric about midpoint. 
for n = 2:nmid-1
    % x-spacing is not constant, so:
    dx1 = x(n)-x(n-1);  
    dx2 = x(n+1)-x(n);
    dxav = (dx1 + dx2)/2;
    %dTdx2 = (T(n+1) - T(n))/dx2;
    %d2Tdx2 = (dTdx2-dTdx1)/dxav;
    %d2Tdx2 = (T(n) - Tl)/Lcool^2 + P(n)/Kel;
    dTdx1 = (T(n )- T(n-1))/dx1;
    T(n+1) = T(n) + (dTdx1 + ( (T(n) - Tl)/Lcool^2 + P(n)/Kel )*dxav)*dx2;
end

T(nmid+1:end)=T(nmid-1:-1:1);  % Because of symmetry

subplot(2,1,1)
plot(x,P),grid
xlabel('x'),ylabel('P')
subplot(2,1,2)
plot(x,T), grid
xlabel('x'),ylabel('temperature')

This results in the following:

[Jungwoo Lee]

Thank you for your answer.