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

4 次查看（过去 30 天）
Kevin Giles
Kevin Giles 2020-10-31
评论： Ravi Kumar 2020-11-3
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

请先登录，再回答此问题。

采纳的回答

Ravi Kumar
Ravi Kumar 2020-11-2
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 个评论
Kevin Giles
Kevin Giles 2020-11-3
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 2020-11-3
Yes, bad elements could also be the problem. I would still suspect the STL and now mesh. You can use quality() function on mesh.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by