Matlab pde toolbox, transient thermal model: state.time-values, used as input for a time-dependent internal heat source do not correspond to defined tlist-values

Hello,

I'm using the pde toolbox and the Matlab version R2020b to model the dissipation of heat, induced by a pulsed heat source in the center of a thin sample to see if stationary heating is an issue.

The following problem persists: Even though the heat source function is working as intended and returns the expected values, the state.time (input of heat source function), which is provided by the solver, does not correspond to the defined tlist-values, which the solver should be using as time steps.

When looking at the state.time-values, instead of increasing by a given step size for each iteration, they appear to increase arbitrary, with the distance between later state.time-values being incredibely large (see statetime_saved.mat in the appendix) .

Thus the later intervals in which the pulsed heat source should be "turned on" are not covered, which completely defies the use of this simulation.

Is there a way, to "force" the solver to use the exact tlist-values or is there another way to work around this?

Main code:

%%% Setup: ultra-thin film sample with the following parameters:

k = 400; %thermal conductivity W/(m-K)
rho = 8960; %density, kg/m^3
specificHeat= 386; %J/(kg-K)

thick = 0.1; %film thickness in meters
hCoeff=1; %Convection coefficient, W/(m^2-K)
Tambient= 30; %Ambient Temperature, K
emiss =0.5; %Emissivity of film surface

%%% Create thermal, transient pde model:
model = createpde('thermal','transient');

%%% Define geometry:
load('gd.mat'); load('sf.mat'); load('ns.mat'); %geometry matrices available in the appendix
g = decsg(gd,sf,ns);

g=geometryFromEdges(model,g);
model.Geometry=extrude(g,thick);

%%% Generate Mesh:
generateMesh(model);

%%% Assign thermal properties:
thermalProperties(model,'ThermalConductivity',k, ...
                        'MassDensity',rho, ...
                        'SpecificHeat',specificHeat);


figure; 
pdegplot(model,'CellLabels','on','FaceAlpha',0.5); 
title 'Geometry With Face Labels Displayed';

The heat is applied to the sample (circular geometry) in the form of cell 6 (C6) acting as a time-dependent internal heat source. Excess heat is carried off via faces 8-11 (connections to square geometry).

%%% Specify Boundary Conditions:
    %Heat Flux from faces 8-11:
    thermalBC(model,'Face',8:11, ...
                    'ConvectionCoefficient',1000, ...
                    'AmbientTemperature',30);
    
    %Heat Source:
        internalHeatSource(model,@pulsedHeatSource,'Cell',6);
        
%%% Specify Initial Conditions:
thermalIC(model,30);

%%% Solve PDE:
tlist = [0:10:5000];
R1 = solve(model,tlist);

%%% Plot course of temperature at center of sample:
plot(R1.SolutionTimes,R1.Temperature(5515,:))

Only one temperature peak visible, where there should be peaks periodically every 1000.

Heat source:

function Q = pulsedHeatSource(location,state)
    pulseLength=100;
    repRate=0.001;
    repTime=1/repRate;

    max=200000; %maximal value
    
    Q = zeros(1,numel(location.x));
    if(isnan(state.time))
        % Returning a NaN when time=NaN tells the solver that the heat source is a function of time.
        Q(:,:) = NaN;
    else
        % Periodic square wave, shifted to positive values only
        Q(:,:) = (max/2*square(2*pi*repRate*floor(state.time),100*pulseLength/repTime))+max/2;
    end
    % Show state.time input for each function call:
    display(state.time)
end

Help would be very much appreciated! Thank you in advance,

Marvin

Answer 1

[Jorge Garcia Garcia]

I have a similar problem with a mechanical time dependant problem returning NaN in some state.time. if you find a solution, let me know

Answer 2

[Ben Halliwell]

With a pulsed heat source, there are large periods when the heat source is constant. The PDE toolbox uses a variable step solver (ode15s), so the step size varies based on the error. When the heat source is constant, the ode solver will be very accurate and the timestep will significantly increase. It might increase so much that the time-step is greater than the length of the pulse, effectively skipping over it, therefore the solver does not 'see' the pulse and so the temperature does not increase.

To solve this, you need to reduce the step size by reducing the error, which can be done by reducing the relative tolerance.

model.SolverOptions.RelativeTolerance = 1e-6;

This is the result I got running your code for a relative tolerance of 1e-6 (the default is 1e-3).