Specify boundary conditions for a thermal model (PDE)

8 次查看(过去 30 天)
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
Alan Moses 2021-2-23
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.
  1 个评论
Pamela Cambianica
Pamela Cambianica 2021-2-26
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))

请先登录,再进行评论。

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by