- 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.
- 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.

# PDE Toolbox producing inconsistent solutions

4 views (last 30 days)

Show older comments

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.

##### 0 Comments

### Accepted Answer

Ravi Kumar
on 18 May 2020

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:

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

##### 0 Comments

### More Answers (0)

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!