Singular Jacobian matrix in steady state problem - PDE toolbox - Error using pde.intern​al.FESolve​rUtilities​/solveStat​ionaryNonl​inear Nonlinear solution failed due to singular

I am trying to solve a heat conduction based steady state problem, however every time I run it, I get the following error, no matter what uniform initial condition I give:

Error using pde.internal.FESolverUtilities/solveStationaryNonlinear

Nonlinear solution failed due to singular Jacobian matrix.

To be honest, my code is more complex than that. I have to spatial domains coupled together. One is solved using a transient solver and the other is solve using a steady state solver. They are coupled by continuity of the heat flux via a function that I created for each time step. The idea is to make calculations faster considering part of the geometry at steady state for each time step of the other part of the geometry.

At a ceratin time, the residual in the steady state region becomes higher than the residual tolerance set by PDE toolbox (1e-4), therefore the steady state solver should start iterating trying to bring it down and converge, but it does not happen, since it crashes.

Here's my code. Geometry and material properties are contained in specific functions:

clear all

clc

close all

% Useful: https://it.mathworks.com/help/pde/ug/pdemodeler-app.html

% ---- Load universal physical constants ---- %

BoltzmannConstant = 1.380649*1e-23; % [J/K]

R = 8.314; % [J/mol*K]

Avog = 6.02214076*1e23; % [atoms/mol]

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

% -------- Start of User input -------- %

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

% ---- Simulation time setting ---- %

tlist = 0:0.2:60*17.3; %0:.5:60*17.3;

% tlist = 0:0.02:60*17.3;

%--- Choose geometry ----%

% Possibilities: 'custom', 'Ponnappan', 'Faghri'

geom = 'Faghri';

[rOut_Wall, eps, rOut_wick, por, rOut_Na, l_first, l_evap, l_adiab,...

l_cond] = chooseGeometry(geom);

% ---- Mesh refinement ----%

% NB: the mesh is triangular over the 2D geometry.

% In this model, the edge of each triangle is fixed by the value of this

% variable.

meshAccuracy = 5;

x = 0:meshAccuracy:l_first+l_evap+l_adiab+l_cond;

% ---- Thermodynamic conditions ----%

% Possibilities: 'custom', 'Ponnappan', 'Faghri'

[T0, Q_heat] = TD_IC(geom);

% ---- Display geometry plots ----%

% When this variable is true, the geometry and mesh plots are displaied for

% the wick-wall region and for the sodium vapor region

display_geometry = 'True';

% ---- Stainless Steel SS304 properties (Wall and solid structure of the wick) ----%

[rho_Wall, k_Wall, Cp_Wall] = wall_properties ();

% ---- Sodium properties (Solid-Liquid-Vapour) ---- % Reference: https://inis.iaea.org/collection/NCLCollectionStore/_Public/43/095/43095088.pdf

[MmNa, P, k_NaS, k_NaL, rho_NaS, rho_NaL, rho_NaG, rho_Namelt,...

Cp_NaS, Cp_NaL, Cp_NaG, Cp_Namelt, visc_Na, Tmelt_Na, Tevap_Na,...

deltaT, H, Hvap, hs, hl, hv] = sodium_properties();

% ---- Surrogate model for the sodium thermal conductivity representing the phase change ----%

% Diffusion coefficients

vBar = @(~,state) sqrt(8*BoltzmannConstant.*state.u./(pi*MmNa))*1e3; % [mm/s] CHECKED

Dk = @(~,state) 2/3*rOut_Na.*vBar(0,state); % [mm^2/s] CHECKED

Dvisc_adiab = @(~,state) rOut_Na^2./(8.*visc_Na(0,state)).*P(0,state); % [mm^2/s] CHECKED

Dvisc_evap = @(~,state) rOut_Na^2./(4.*visc_Na(0,state)).*P(0,state); % [mm^2/s] CHECKED

D_Na_adiab = @(~,state) Dk(0,state) + Dvisc_adiab(0,state); % [mm^2/s] CHECKED

D_Na_evap = @(~,state) Dk(0,state) + Dvisc_evap(0,state); % [mm^2/s] CHECKED

% Effective thermal conductivities

k_evap = @(~,state) D_Na_evap(0,state).*(hl(0,state)*Hvap*MmNa.*P(0,state)*1e-6)./(R/MmNa*Avog*BoltzmannConstant*state.u.^3); % [W/mm*K]

k_adiab = @(~,state) D_Na_adiab(0,state).*(hv(0,state)*Hvap*MmNa.*P(0,state)*1e-6)./(R/MmNa*Avog*BoltzmannConstant*state.u.^3); % [W/mm*K]

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

% -------- End of User input -------- %

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

%% Transient and steadystate models creation

% PDE model creation - Transient region

thermalModelHP = createpde('thermal','transient-axisymmetric'); % page 5-243 of the guide

thermalModelHP.StefanBoltzmannConstant = 5.670373e-8*1e-6; % Stefan-Boltzmann constant, W/(mm^2-K^4)

[faceDomainsTrans, nodes_wallWick] = build_wickWall(thermalModelHP, meshAccuracy, rOut_Na, rOut_wick,...

rOut_Wall, l_first, l_evap, l_adiab, l_cond, display_geometry);

% PDE model creation - Steadystate region

thermalModelHP_steady = createpde('thermal','steadystate-axisymmetric'); % page 5-243 of the guide

% Create customized mesh

[nodes_Na, triangles_Na, faceDomainsSteady] = customMesh(thermalModelHP_steady, meshAccuracy, l_first+l_evap+l_adiab+l_cond, rOut_Na, display_geometry);

%% Define thermal properties %%

% Wall

thermalConductivity.Wall = k_Wall;

spHeat.Wall = Cp_Wall;

density.Wall = rho_Wall;

% Wick

kNa_melt = @(~,state) k_NaS + ...

(k_NaL(0,state) - k_NaS).*(state.u - Tmelt_Na - deltaT)./(2*deltaT);

CpMelt = @(~,state) (Cp_NaS + Cp_NaL(0,state))/2 + H/(2*deltaT);

% CpEvap = @(~,state) (Cp_NaL(0,state) + Cp_NaG(0,state))/2 + Hvap/(2*deltaT); In theory, we assume that no evaporation happens in the wick

thermalConductivity.Wick = @(~,state) (k_NaS*(k_NaS + k_Wall(0,state) - (1-por)*(k_NaS - k_Wall(0,state)))./...

(k_NaS + k_Wall(0,state) + (1-por)*(k_NaS - k_Wall(0,state)))).*(state.u < (Tmelt_Na - deltaT))...

+ (kNa_melt(0,state).*(kNa_melt(0,state) + k_Wall(0,state) -...

(1-por)*(kNa_melt(0,state) - k_Wall(0,state)))./...

(kNa_melt(0,state) + k_Wall(0,state) + (1-por)*(kNa_melt(0,state) - k_Wall(0,state))))...

.*(state.u >= (Tmelt_Na - deltaT) & state.u <= (Tmelt_Na + deltaT))...

+ (k_NaL(0,state).*(k_NaL(0,state) + k_Wall(0,state)...

- (1-por)*(k_NaL(0,state) - k_Wall(0,state)))./...

(k_NaL(0,state) + k_Wall(0,state) + (1-por)*(k_NaL(0,state) - k_Wall(0,state))))...

.*(state.u > (Tmelt_Na + deltaT));

spHeat.Wick = @(~,state) (por*Cp_NaS + (1-por)*Cp_Wall(0,state)).*(state.u < (Tmelt_Na - deltaT)) + ...

(por*CpMelt(0,state) + (1-por)*Cp_Wall(0,state)).*(state.u >= (Tmelt_Na - deltaT) & state.u <= (Tmelt_Na + deltaT)) + ...

(por*Cp_NaL(0,state) + (1-por)*Cp_Wall(0,state)).*(state.u > Tmelt_Na + deltaT);

% (por*Cp_NaL(0,state) + (1-por)*Cp_Wall(0,state)).*(state.u > Tmelt_Na + deltaT & state.u < Tevap_Na - deltaT) + ...

% (por*CpEvap(0,state) + (1-por)*Cp_Wall(0,state)).*(state.u >= (Tevap_Na - deltaT) & state.u <= (Tevap_Na + deltaT)) + ...

% (por*Cp_NaG(0,state) + (1-por)*Cp_Wall(0,state)).*(state.u > (Tevap_Na + deltaT));

density.Wick = @(~,state) (por*rho_NaS + (1-por)*rho_Wall).*(state.u < Tmelt_Na) + ...

(por*rho_Namelt + (1-por)*rho_Wall).*(state.u >= (Tmelt_Na - deltaT) & state.u <= (Tmelt_Na + deltaT)) + ...

(por*rho_NaL(0,state) + (1-por)*rho_Wall).*(state.u > Tmelt_Na + deltaT);

% Sodium

thermalConductivity.Na = @(location,state) ...

k_adiab(0,state).*(location.y >= 0 & location.y < l_first) + ...

k_evap(0,state).*(location.y >= l_first & location.y <= (l_first+l_evap) & location.y >= (l_first+l_evap+l_adiab) & location.y <= (l_first+l_evap+l_adiab+l_cond)) + ...

k_adiab(0,state).*(location.y > l_first+l_evap & location.y < (l_first+l_evap+l_adiab));

spHeat.Na = @(~,state) Cp_NaG(0,state);

density.Na = @(~,state) rho_NaG(0,state);

%% Define constant boundary conditions %%

% Heat flux on evaporator wall

fluxSS = @(~,state) Q_heat/(2*pi*rOut_Wall*l_evap); % [W/mm^2]

% Boundary condition to be applied for the first iteration ONLY on the left

% boundary of the wick-wall structure: it assumes steady state heat

% conduction (same heat flux on the right and left surfaces). It can be

% improved most probably, but I have run out of ideas.

flux = @(location,state) -Q_heat/(2*pi*rOut_Wall*l_evap).*(location.y >= l_first & location.y <= (l_first+l_evap)) + ...

0.*(location.y > (l_first+l_evap));

% BC definition for the PDE toolbox - TRANSIENT REGION

thermalBC(thermalModelHP,'Edge',7,'AmbientTemperature',293,...

'Emissivity', eps); % Condenser's edge

thermalBC(thermalModelHP,'Edge',4, 'HeatFlux', 0); % First adiabatic's edge

thermalBC(thermalModelHP,'Edge',6, 'HeatFlux', 0); % Adiabatic's edge

thermalBC(thermalModelHP,'Edge',5, 'HeatFlux', fluxSS); % Evaporator's edge

thermalBC(thermalModelHP,'Edge',9, 'HeatFlux', 0); % Condenser end plate - Wall

thermalBC(thermalModelHP,'Edge',8, 'HeatFlux', 0); % Condenser end plate - Wick

thermalBC(thermalModelHP,'Edge',3, 'HeatFlux', 0); % Evaporator end plate - Wall

thermalBC(thermalModelHP,'Edge',2, 'HeatFlux', 0); % Evaporator end plate - Wick

% BC definition for the PDE toolbox - STEADY STATE REGION

thermalBC(thermalModelHP_steady,'Edge',2, 'HeatFlux', 0); % Evaporator end plate - Na

thermalBC(thermalModelHP_steady,'Edge',3, 'HeatFlux', 0); % Condenser end plate - Na

%% Definition of the weakly-coupled solver %%

% Initialize solution as a struct

tempDistr = struct();

tempDistr.T_wickWall = T0*ones(size(nodes_wallWick,2),length(tlist)); %zeros(nodes_trans,length(tlist))% CHECK, MAYBE NOT OK TO CALL THIS AS AN ARRAY, MAYBE IT IS A STRUCT

tempDistr.T_Na = T0*ones(size(nodes_Na,2),length(tlist));

% It takes long and it is prone to instabilities...

for ii = 1:length(tlist)-1

% Set the non-constant BC & IC definition - Transient region

if (ii==1)

thermalIC(thermalModelHP, unique(tempDistr.T_wickWall(:,1)));

else

thermalIC(thermalModelHP, result_wickWall);

end

% WHAT CONDITION TO GIVE HERE? MAYBE ASSUMING CONTINUOUS FLUX AS FIRST

% INPUT? WATCH OUT ALSO FOR AXISYMMETRIC BC OF THE MODEL: IS IT RIGHT

% FOR AN HOLLOW CYCLINDER?

thermalBC(thermalModelHP,'Edge',1, 'HeatFlux', flux); % Symmetry (dudx = 0): not needed because the problem is specified as axisymmetric

IDX = 1;

for domain = faceDomainsTrans

thermalProperties(thermalModelHP, ...

ThermalConductivity=thermalConductivity.(domain), ...

SpecificHeat=spHeat.(domain), ...

MassDensity=density.(domain), ...

Face=IDX);

IDX = IDX+1;

end

% Solve the transient case

thermalModelHP.SolverOptions.RelativeTolerance = 1e-2;

thermalModelHP.SolverOptions.AbsoluteTolerance = 1e-4;

thermalModelHP.SolverOptions.ResidualTolerance = 1e-3;

thermalModelHP.SolverOptions.MaxIterations = 50;

thermalModelHP.SolverOptions.ReportStatistics = 'on';

fprintf('Statistics for wick-wall transient calculations at %f\n', tlist(ii))

result_wickWall = solve(thermalModelHP,[tlist(ii) tlist(ii+1)]);

fprintf('\n')

% Update the struct containing the temperature distribution for each

% time step - transient region

tempDistr.T_wickWall(:,ii+1) = result_wickWall.Temperature(:,2); % MAYBE HERE A STRUCT IS NEEDED, NOT AN ARRAY

% -------------------------------------------------------------------------------------------------------------------------%

% Set the non-constant BC & IC definition - Steady state region

if (ii==1)

thermalIC(thermalModelHP_steady, unique(tempDistr.T_Na(:,1)));

else

thermalIC(thermalModelHP_steady, result_Na);

end

xPoints = thermalModelHP.Mesh.Nodes(2,find(thermalModelHP.Mesh.Nodes(1,:)==rOut_Na));

status.u = result_wickWall.Temperature(find(thermalModelHP.Mesh.Nodes(1,:)==rOut_Na),2);

KAPPA = k_evap(0,status);

gradT_Na = result_wickWall.RGradients(find(thermalModelHP.Mesh.Nodes(1,:)==rOut_Na),2);

iter_flux = -KAPPA.*gradT_Na;

flux = @(location,~) interp1(xPoints, iter_flux, location.y, 'linear', 'extrap');

% Spatial non-constant BC for steady state region

% AM I UPDATING RIGHT THE FLUX?

thermalBC(thermalModelHP_steady,'Edge',1, 'HeatFlux', flux,...

'Vectorized', 'on'); % Na-wick interface

% thermalBC(thermalModelHP_steady,'Edge',4, 'HeatFlux', 0); % Symmetry (dudx = 0): not needed because the problem is specified as axisymmetric

IDX = 1;

for domain = faceDomainsSteady

thermalProperties(thermalModelHP_steady, ...

ThermalConductivity=thermalConductivity.(domain), ...

Face=IDX);

IDX = IDX+1;

end

thermalModelHP_steady.SolverOptions.ReportStatistics = 'on';

fprintf('Statistics for Na vapor core steady state calculations at %f\n', tlist(ii));

result_Na = solve(thermalModelHP_steady);

fprintf('\n')

% Update the struct containing the temperature distribution for each

% time step - steady state region

tempDistr.T_Na(:,ii+1) = result_Na.Temperature(:,1);

xPointsSS = thermalModelHP_steady.Mesh.Nodes(2,find(thermalModelHP_steady.Mesh.Nodes(1,:)==rOut_Na));

statusSS.u = result_Na.Temperature(find(thermalModelHP_steady.Mesh.Nodes(1,:)==rOut_Na),1);

KAPPASS = k_evap(0,statusSS);

gradT_NaSS = result_Na.RGradients(find(thermalModelHP_steady.Mesh.Nodes(1,:)==rOut_Na),1);

iter_fluxSS = +KAPPASS.*gradT_NaSS;

flux = @(location,~) interp1(xPointsSS, iter_fluxSS, location.y, 'linear', 'extrap');

end

%%

T_wickWall = result_wickWall.Temperature;

figure

pdeplot(thermalModelHP,'XYData',T_wickWall(:,2),'ZData',T_wickWall(:,2))

T_Na = result_Na.Temperature;

figure

pdeplot(thermalModelHP_steady,'XYData',T_Na(:,1),'ZData',T_Na(:,1))

[Torsten]

And the code works if you model all domains as unsteady ?