Specify boundary conditions for a thermal model (PDE)

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)

Answer 1

[Alan Moses]

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]

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
[dist]=xlsread('distance24NovAU.xls'); 
%Incidence angles
[in]=xlsread('incidence');
%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))