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

20 views (last 30 days)

Show older comments

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

##### 0 Comments

### Accepted Answer

Ravi Kumar
on 2 Nov 2020

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

##### 2 Comments

Ravi Kumar
on 3 Nov 2020

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!