Transient thermal model (PDE Toolbox) failing to solve due to time dependent heat flux and ambient temperature. "Unable to meet integration tolerance."

Hello,

I am currently trying to develop a thermal model of an object (cryobag) which has a changing ambient temperature (changing freezer temperature). Because of this I have realised I cannot use the convective convective coefficient and emissivity thermalBC because the ambient temperature must be a fixed value. So instead I have created a seperate function for the heat flux and will use in the heat flux thermalBC.

For this, I have the function for changing ambient temperature built in. However the surface of the object is also changing which is required for the convection and radiative heat flux. I am currently using state.u as the temperature of the surface but I dont know if this is correct as I think it may be solving for every tempature point within the object which is not required. The model is solving at extremely small time steps (1e-16) where as I only want it to solve at 5 second intervals for example. Is there any other ways to do this or am I going down the complete wrong route?

Any help would be greatly appreciated.

I am getting these errors when the code is run:

Warning: Failure at t=2.246105e+00. Unable to meet integration tolerances without reducing the step size below the smallest value

allowed (7.105427e-15) at time t.

> In ode15s (line 655)

In pde.EquationModel/solveTimeDependent (line 101)

In pde.ThermalModel/solve (line 91)

In test (line 36)

Error using pde.EquationModel/solveTimeDependent (line 103)

Solution failed to reach the requested end time.

Error in pde.ThermalModel/solve (line 91)

u = self.solveTimeDependent(coefstruct,u0,[],tlist,false);

Here is the heat flux function and script:

function q = heatflux(location, state)
time=state.time;
%disp(time)
%disp(location);
%disp(state.u)

k_air = 24.35e-3 ; % Thermal conductivity, W/(m*C)
rho_air = 1.225; % Density, kg/m^3
cp_air = 1004; % Specific heat, W*s/(kg*C)
mew_air = 1.729e-5; %dynamic viscosity (kg/m*s)
kv_air = 1.338e-5; % kineamtatic viscosity m2/s
g = 9.81; %accelration due to gravity (m/s^2)
emis = 0.1; %emissitivity
L = 1.6; %characteristic length
stefan = 5.670367e-8; % W/(m^2 K^4)
Ta = 193.15; % Freezer temp
Tsur = state.u; % Surface temp
%beta = 1/Tsur; %thermal expansion coefficient (K^-1)

if  (time>=0 && time<20)
    Tamb=277.15;
elseif  (time>=20 && time<40)
    Tamb= 263.15;
elseif  (time>=40 && time<50)
    Tamb = 233.15 - (2/60)*(time-40);
else
    Tamb = 277.15;
end

beta= 2/(Tsur+Tamb);

qr = emis*stefan*(Tsur^4 - Tamb^4);

% qc = convective heat transfer
Pr = cp_air*mew_air/k_air; %Prandtl Number (of fluid should it be air?)
Gr = (g*beta*(Tsur-Tamb)*L^3)/(kv_air^2); % Grashof Number
Ra = Pr*Gr; %Rayleigh Number
Nu = (0.68 + ((0.67*Ra^(1/4))/((1 + (0.492/Pr)^(9/16))^(4/9)))); %Nusselt Number
hc = (Nu*k_air)/L; %convective heat transfer coefficient (W/m2 K)
qc = hc*(Tsur-Tamb); %covective heat transfer
q = qr + qc; %total heat flux/ Boundary condition

end

thermalmodel = createpde('thermal', 'transient');
importGeometry(thermalmodel,'cryobagscaleup2.stl');
pdegplot(thermalmodel,'FaceLabel', 'on', 'CellLabel', 'on', 'FaceAlpha',0.5)
mesh = generateMesh(thermalmodel,'Hmax', 0.08);
pdeplot3D(thermalmodel)

vc = 0.075; % Volume fraction of DMSO
k_bag = 0.22411*vc^2 - 0.64025*vc + 0.61579 ; % Thermal conductivity, W/(m*C)
rho_bag = 1000; % Density , kg/m^3
cp_bag = 4182; % Specific heat, W*s/(kg*C)
g = 9.81; %accelration due to gravity (m/s^2)
emis = 0.1; %emissitivity
L = 1.6; %characteristic length
thermalmodel.StefanBoltzmannConstant = 5.670367e-8; % W/m^2 K^4
stefan = 5.670367e-8; % W/(m^2 K^4)


%location values
%location.x = -0.0855799;
%location.y = -0.258654;
%location.z = 0.0941327;

heatfluxfunc = @(location,state) heatflux(location,state);

thermalProperties(thermalmodel,'ThermalConductivity',k_bag,'MassDensity', rho_bag, 'SpecificHeat',cp_bag);
internalHeatSource(thermalmodel,0);

thermalIC(thermalmodel,293.15);
%thermalBC(thermalmodel,'Face',1:thermalmodel.Geometry.NumFaces,'ConvectionCoefficient',hc ,'AmbientTemperature',Tf);
%thermalBC(thermalmodel,'Face',1:thermalmodel.Geometry.NumFaces,'Emissivity',emis ,'AmbientTemperature',Tf);
thermalBC(thermalmodel,'Face',1:thermalmodel.Geometry.NumFaces,'HeatFlux', heatfluxfunc);

time = [0:5:50];
R = solve(thermalmodel,time);
Tmin = 193.15;
Tmax = max(R.Temperature(:,end));

 
h = figure;
for i = 1:numel(time)
    pdeplot3D(thermalmodel,'ColorMapData',R.Temperature(:,end))
    caxis([Tmin,Tmax])
    view(130,10)
    title(['Temperature at Time (s)' num2str(time(i))]);
    M(i) = getframe;
   
end

Answer 1

[Ravi Kumar]

Hi Kevin,

Make sure your geometry from STL is in SI units too as you have properties specified in SI units. I don't see anything obvious from the code that could trigger instability in transient solution. It might be good idea to solve a steady-state solution to see if everything is put-together correctly.

Regards,

Ravi

[Kevin Giles]

Hi Ravi,

Thanks for the reply. I am also getting a warning about the mesh:

Warning: Found elements with poor shape quality.

> In pde.EquationModel/generateMesh (line 81).

Could this be the problem more so than the code itself?

I will try it in a steady state solution.

Thanks,

Kevin

[Ravi Kumar]

Yes, bad elements could also be the problem. I would still suspect the STL and now mesh. You can use quality() function on mesh.