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

16 Ansichten (letzte 30 Tage)
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

Akzeptierte Antwort

Ravi Kumar
Ravi Kumar am 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 Kommentare
Kevin Giles
Kevin Giles am 3 Nov. 2020
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
Ravi Kumar am 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte


Version

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by