Specify boundary conditions for a thermal model (PDE)

7 views (last 30 days)
Pamela Cambianica on 19 Feb 2021
Commented: Pamela Cambianica on 26 Feb 2021
Hello everyone,
I need to simulate the heat transfer on a surface by using the PDE Toolbox. I have defined the geometry, the mesh, thermal properties of the material, the internal heat source, but I don't know how to define the following boundary conditions:
• Surface boundary condition (z=0): where S is the heat flux proportional to the incidence angle (i), sigma is the Stefan-Boltzmann constant, epsilon is the emissivity, k the thermal conductivity.
• Bottom boundary condition (z=D): where Q is the heat production of the body
Does anyone have any ideas?
Thanks for any help!
Pamela
thermalmodel = createpde('thermal','transient');
%% Geometry
gm = multicuboid(1,1,[0.1 1.9],'ZOffset',[0 0.1]);
thermalmodel.Geometry = gm;
gm2 = rotate(gm,180,[0 0 0],[0 1 0]);
thermalmodel.Geometry = gm2;
%pdegplot(thermalmodel,'CellLabels','on','FaceAlpha',0.5)
%% Generate a mesh for the geometry
msh = generateMesh(thermalmodel,'Hmax', 5);
%% Specify thermal properties of the material
%--------------------------------------------------------------------------
%DENSITY (rho)
rhob = 1000;
rho = @(location,state)rhob*((location.z + 0.122)/(location.z + 0.18));
%SPECIFIC HEAT (cp)
ca=1;
cb=3;
cc=5;
cd=0.4;
ce=9;
cp = @(location,state)ca+cb*state.u+cc*(state.u)^2+cd*(state.u)^3+ce*(state.u)^4;
%THERMAL CONDUCTIVITY (k)
kb = 0.8;
ks = 0.7;
rhos = 2100;
sigma = 8.67e-08;
epsilon = 0.95;
l = 5.5e-05;
beta = 4*sigma*epsilon*l;
kc = @(location,state)kb-(kb-ks)*((rhob-(rho(location,state))))/(rhob-rhos);
k = @(location,state)beta*(state.u)^3+(kc(location,state));
thermalProperties(thermalmodel,'ThermalConductivity',k, 'MassDensity',rho,'SpecificHeat',cp)
%% Internal heat source [W m^-2]
Q = 0.020
internalHeatSource(thermalmodel,Q)

Alan Moses on 23 Feb 2021
You may refer the applyBoundaryCondition function and the solvepde function to see if it fits your requirement. You may also refer this example to use the functions mentioned.
You may refer to a similar solution here.
Pamela Cambianica on 26 Feb 2021
Dear Alan,
thanks a lot for your help. As you can see in the following code, I have included the applyBoundaryCondition function. However, the problem now is that any thermal function, such as
thermalProperties(thermalmodel,'ThermalConductivity',k, 'MassDensity',rho,'SpecificHeat',cp)
or
internalHeatSource(thermalmodel,Q)
gives me the following error: "Check for missing argument or incorrect argument data type in call to function 'internalHeatSource'."
How can I use these functions in a generic model = createpde(N)?
Thanks again for helping me out.
Pamela
thermalmodel = createpde(1);
gm = multicuboid(1,1,[0.1 1.9],'ZOffset',[0 0.1]);
gm = rotate(gm,180,[0 0 0],[0 1 0]);
thermalmodel.Geometry = gm;
pdegplot(thermalmodel,'FaceLabels','on');
axis equal;
%% Specify thermal properties of the material
%--------------------------------------------------------------------------
%DENSITY (rho)
%rhob: density at the bottom boundary [Kg m^-3]
%z is the depth [m]
rhob = 3100;
rho = @(location,state)rhob*((location.z + 0.122)/(location.z + 0.18));
%--------------------------------------------------------------------------
%SPECIFIC HEAT (cp)
ca=-9.05e+05;
cb=5.14e+00;
cc=-6.73e-06;
cd=7.67e-08;
ce=-2.38e-10;
cp = @(location,state)ca+cb*state.u+cc*(state.u.^2.)+cd*(state.u.^3.)+ce*(state.u.^4.);
%--------------------------------------------------------------------------
%THERMAL CONDUCTIVITY (k)
kb = 0.8;
ks = 0.7;
rhos = 2236;
sigma = 5.67e-08;
epsilon = 0.95;
l = 7.5e-05;
beta = 4*sigma*epsilon*l;
kc = @(location,state)kb-((kb-ks)*((rhob-(rho(location,state)))/(rhob-rhos)));
k = @(location,state)beta*(state.u.^3.)+(kc(location,state));
%% Define PDE
d = @(location,state)rho(location,state)*cp(location,state)*(-1);
specifyCoefficients(thermalmodel,'m',0,'d',d,'c',k,'a',0,'f',0);
%% S: absorbed shortwave radiation from the Sun
A = 0.06;
FO = 1367;
C=(1-A)*FO;
%distance Sun-Surface
%Incidence angles
%plot(i,R)
B=zeros(100);
S=zeros(100);
for i = 1:length(in)
for R = 1:length(dist)
B(R,i)=cos(i)/(R.^2);
S(R,i)=B(R,i)*C;
end
end
%% Set boundary conditions
BBC=applyBoundaryCondition(thermalmodel,'neumann','Face',7,'g',0.02);
SurfaceFunction = @(location,state)epsilon*sigma*(state.u.^4)-S(R,i);
SBC=applyBoundaryCondition(thermalmodel,'neumann','Face',1,'g',SurfaceFunction);
%% Generate Mesh
generateMesh(thermalmodel,'Hmax', 0.5);
pdemesh(thermalmodel);
axis equal
%% Set initial conditions
setInitialConditions(thermalmodel,400)
internalHeatSource(thermalmodel,10) %ERROR
tlist = 0:10;
results = solvepde(thermalmodel,tlist);
pdeplot3D(thermalmodel,'ColorMapData',results.NodalSolution(:,2))