Transient thermal model (PDE Toolbox) failing to solve due to time dependent heat flux and ambient temperature. "Unable to meet integration tolerance."
12 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
Yes, bad elements could also be the problem. I would still suspect the STL and now mesh. You can use quality() function on mesh.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!