PDE Toolbox producing inconsistent solutions

1. The associated geometry and mesh for my problem is given below:

2. The PDE is as follows:

General form accepted by PDE Toolbox (as detailed in documentation):

where such that the PDE reduces to the heat diffusion equation with no generation:

where are the density, specific heat capacity and thermal conductivity respectively (and all are constant). Of course in this scenario, the dependent variable u represents the temperature.

Example main function code for solving this in the PDE toolbox is given below (please excuse its length, and the fact that I have it as a function):

function [model,results,transition,tlist] = lumpedmain()

%% Clear MATLAB
clc;clear;close all;clear global;
format LONGG

%% Declare global variables
global m            %data matrix
global n            %number of orbits
global tlist        %time domain       
global transition   %transition elements
global packX        %pack x dimension
global packY        %pack y dimension
global packZ        %pack z dimension
global qsolar       %solar radiative heat rate
global qalbedo      %albedo radiative heat rate
global qplanetary   %planetary radiative heat rate
global e            %emmitance
global a            %absorptance
global sigma        %stefan-boltzmann constant
global results      %results

%% Create PDE Model
model = createpde;  %N=1 system of scalar equations

%% Create mesh geometry (lumped pack)
% Set global variables
packX = 80;
packY = 40;
packZ = 65;
% Create basic shapes
R1 = [3;4;0;packX;packX;0;0;0;packY;packY];  %pack
%
gd = [R1];  %combine to form gd matrix
% Create names for basic shapes
ns = char('R1');    %create name-space matrix
ns = ns';   %transpose to vector
% Set formula
sf = 'R1';
% Create geometry
g = decsg(gd,sf,ns);
% Append geometry to model
geometryFromEdges(model,g);
% Create mesh
mesh = generateMesh(model,'Hmin',5,'Hmax',20);  %create triangular mesh    
%View geometry w/o mesh
figure
pdegplot(g,'EdgeLabels','on','FaceLabels','on')
xlabel('x (mm)')
ylabel('y (mm)')
grid on
xlim([0 packX])
ylim([0 packY])
% View geometry w/mesh
figure
pdemesh(mesh)
xlabel('x (mm)')
ylabel('y (mm)')
xlim([0 packX])
ylim([0 packY])

%% Set up tlist
% Import Data
c = struct2cell(importdata('200113 - Drive Cycle Test 1_CB1.txt'));
% Convert into Matrix 
c = c(1,1);
fm = cell2mat(c);
% Chop matrix to focus on orbital period
m = fm(11633:24177,:);
% Reset time
x = m(1,2);
for i = 1:length(m)
    m(i,2) = m (i,2) - x;
end
% Chop m to ensure same voltage at start and end of orbit
m = m(1:11577,:);
m(end,3) = m(1,3);

%% Specify time domain by number of orbits
% Set number of orbits
n = 200;
% Define a single orbit
orbit = m(:,2); %time domain is defined by the data points in drive cycle
%
tlistrough = orbit;
% For more than one orbit
if n > 1
    for i = 2:n
        tlistrough = [tlistrough;(orbit+(i-1)*orbit(end))];
    end
end
% Set up vector that monotonically increases across the domain set above
tlist = zeros(length(tlistrough),1);
dt = 0.500000023739631;   %timestep
%
tlist = linspace(tlistrough(1),tlistrough(end),length(tlist));
tlist = tlist';
%
% Combine m array for n orbits
mlong = m;
if n > 1
    for i = 2:n
        mlong = [mlong;m];
    end
end
% Catch element where eclipse transition occurs
j = 1;  %set counter for transition vector
for i = 2:(length(mlong)-1)
    dI = mlong(i,4)-mlong(i-1,4);
    if abs(dI) > 143
        transitionrough(j,1) = i-1;
        j = j + 1;
    end
end
% Add beginning and end elements
transitionrough(1) = 1;
transitionrough(length(transitionrough)+1) = length(tlist);
% Catch any errors
dx = 0;
j = 1;
for i = 1:(length(transitionrough)-1)
    dx = transitionrough(i+1)-transitionrough(i);
    if abs(dx) > 10
        transition(j) = transitionrough(i);
        j = j + 1;
    end
    transition(j) = transitionrough(end);
end
transition = transition';

%% Compute radiation components
% Incident radiation intensities
albedo = 0.35;   %planetary albedo of Earth (Fortescue, P. et al.)
F = 0.7;   %visibility factor
Js = 1371;  %solar radiation intensity (W/m^2)
Ja = Js*albedo*F;   %albedo radiation intensity (W/m^2)
z = 550000; %orbital altitude (m)
Rrad = 6371000; %radius of earth
Rorbit = Rrad+z;
Jp = 237*(Rrad/Rorbit)^2;   %planetary radiation intensity (W/m^2)
% Radiation parameters
sigma = 5.67*(10^(-8)); %Stefan-Boltzman constant
e = 0.63;  %emittance
a = 0.40;  %absorptance
% Radiation components (W/m)
qsolar = a*Js*(packZ/1000);
qalbedo = a*Ja*(packZ/1000);
qplanetary = e*Jp*(packZ/1000);

%% Specify model coefficients
k = 9;
rho = 1792.5;
cp = 1395;
%
specifyCoefficients(model,'m',0,...
                          'd',rho*cp*(packZ/1000),...
                          'c',k*(packZ/1000),...
                          'a',0,...
                          'f',0,...
                          'Face',1); 

%% Specify boundary conditions
% Isolated faces
applyBoundaryCondition(model,'neumann','Edge',1,'q',0,...
                                                'g',0);
%
applyBoundaryCondition(model,'neumann','Edge',3,'q',0,...
                                                'g',0);
% Exposed faces
applyBoundaryCondition(model,'neumann','Edge',2,'q',0,...
                                                'g',@gcoef_sun);    %sun-facing
%
applyBoundaryCondition(model,'neumann','Edge',4,'q',0,...
                                                'g',@gcoef_planet); %planet-facing

%% Specify initial conditions
T0 = 318;   %initial temperature (K)
setInitialConditions(model,T0); %set T0 for whole domain

%% Solve PDE
results = solvepde(model,tlist);

%% Find boundary nodes
nodes_right = findNodes(results.Mesh,'region','Edge',2);
nodes_left = findNodes(results.Mesh,'region','Edge',4);

%% Store nodal solutions according to boundary
% Right boundary (E2)
for i = 1:length(nodes_right)
    T_right(i,:) = results.NodalSolution(nodes_right(i),:);
end
% Left boundary (E4)
for i = 1:length(nodes_left)
    T_left(i,:) = results.NodalSolution(nodes_left(i),:);
end

%% Define new domain
numorb = tlist/5787.47027489007;

%% Plot results
% Average pack temperature
figure
plot(numorb,mean(results.NodalSolution)-273.15,'r')
ylabel('Mean pack temperature (\circC)')
xlabel('Number of Orbits')
xlim([0,numorb(end)]);
grid on
grid minor
% Boundary temperatures
figure
plot(numorb,mean(T_right)-273.15)
ylabel('Mean Temperature (\circC)')
xlabel('Number of Orbits')
xlim([0,numorb(end)]);
grid on
grid minor
hold on
plot(numorb,mean(T_left)-273.15)
lgd2 = legend('Sun-facing','Earth-facing','location','best');

end

NOTE: I have attached the data '200113 - Drive Cycle Test 1_CB1.txt' needed when setting up tlist to this thread.

3. The boundary conditions are as follows:

E1 - Isolated, heat flux = 0

E2 - Neumann condition where

E3 - Isolated, heat flux = 0

E4 - Neumann condition where

where are the emmisivity and Stefan-Boltzmann constant, respectively. is the depth into the page, applied throughout to convert into W/m.

These boundary conditions are supposed to simulate a body in orbit around the Earth, such that the inward radiation flux from the sun on E2, and the inward planetary and albedo radiation flux on E4 vary with position (and thus time). The outward radiation from the body is always present.

Example code for the functions gcoef_sun and gcoef_planet that (attempt) to implement the Neumann conditions are below:

function [g_sun] = gcoef_sun(location,state)

%% Declare global variables
global n            %number of orbits
global tlist        %time domain       
global transition   %transition elements
global qsolar       %solar radiative heat rate
global e            %emmitance
global sigma        %stefan-boltzmann constant
global packZ        %pack z dimension

%% Declare output variable
g_sun = zeros(1,length(location.y)); %output must have the form [NxM]

%% Computation
% The temperature returned depends on the solution time.
if(isnan(state.time))
    % Returning a NaN for g when time=NaN
    % tells the solver that the boundary conditions are functions of time.
    % The PDE Toolbox documentation discusses this requirement in more detail.
    g_sun = NaN;
else
    if state.time == 0
        g_sun(1,:) = (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000);  %eclipse
    else
        for i = 2:2:2*n
            if state.time > tlist(transition(i-1,1),1) && state.time <= tlist(transition(i,1),1)
                g_sun(1,:) = (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000); %eclipse
            elseif state.time > tlist(transition(i,1),1) && state.time <= tlist(transition(i+1,1),1)
                g_sun(1,:) = qsolar + (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000);  %sunlight
            end
        end
    end
end

end

function [g_planet] = gcoef_planet(location,state)

%% Declare global variables
global n            %number of orbits
global tlist        %time domain       
global transition   %transition elements
global qalbedo      %albedo radiative heat rate
global qplanetary   %planetary radiative heat rate
global e            %emmitance
global sigma        %stefan-boltzmann constant
global packZ        %pack z dimension

%% Declare output variable
g_planet = zeros(1,length(location.y)); %output must have the form [NxM]

%% Computation
% The temperature returned depends on the solution time.
if(isnan(state.time))
    % Returning a NaN for g when time=NaN
    % tells the solver that the boundary conditions are functions of time.
    % The PDE Toolbox documentation discusses this requirement in more detail.
    g_planet = NaN;
else
    if state.time == 0
        g_planet(1,:) = qplanetary + (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000);  %eclipse
    else
        for i = 2:2:2*n
            if state.time > tlist(transition(i-1,1),1) && state.time <= tlist(transition(i,1),1)
                g_planet(1,:) = qplanetary + (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000); %eclipse
            elseif state.time > tlist(transition(i,1),1) && state.time <= tlist(transition(i+1,1),1)
                g_planet(1,:) = qplanetary + qalbedo + (-1)*(e)*(sigma)*((state.u)^4)*(packZ/1000);  %sunlight
            end
        end
    end
end
end

4. The results are as follows:

For n=100 and n=200 orbits (Note - for n=100 the simulation takes ~30 minutes, for n=200 it takes in excess of 3 hours. The exact number of orbits is specified within the %% Specify time domain by number of orbits section of the main function, set this accordingly to run the simulation for as long as desired):

I should note that each orbit corresponds to around 5790 s (just under 100 minutes), and here I am solving in steps of ~0.5 s (due to the frequency at which data was sampled - see main function above).

5. Issue

Clearly these two results are not consistent. In the top two plots it seems as if the system is approaching steady-state, but in the bottom two plots there is this strange behaviour that begins around the 70th orbit that is not there for the simulation with n=100. The only thing I have changed between running these two simualtions was the variable n, nothing else. For whatever reason MATLAB is misbehaving.

What could be the cause of this misbehaviour?

I suspect that my boundary conditions may be the issue, but if somebody could please shed some light on this I would greatly appreciate the help. Apologies for the long thread.

Answer 1

[Ravi Kumar]

Hi Atdhe,

Here is a quick fix, tighten the tolerance, this will force ODE solver to take finer steps. Insert the following before solvepde call:

model.SolverOptions.RelativeTolerance = 1E-6;
model.SolverOptions.AbsoluteTolerance =  1E-8;

I got the following results for n = 200.

I think code can be cleaned up quite a bit by:

1. Completely removing globals. You can pass additional variables to g-functions, by wrapping them in an anonymous function. For example, gcoef_sun can be wrapped in gSun = @(location,state) gcoef_sun(location, state, n, transition, ...). This will make the debugging much easier.
2. Use vector operation to speed up and better readability. You can replace most of your for-loop with vector operation, also run 'on' vectorized operation in applyBoundaryCondition, check the doc page for details.

On the modeling aspect, you seem to be combining the continuous ODE system with the discrete transition. A miss in one of the transitions, as ODE takes variable steps, the nonlinear solution can diverge or converge to a different solution which is what you observed with n = 200. I don't know how to fix this modeling aspect, but I think there is a better way to model this instead of two branches in the g-function.

Regards,

Ravi