The Finite Element Method of Flow and Heat Transfer in
Heterogeneous Materials
ROSLIANA ESO1, ARMAN2
1Department of Physics Education, Universitas Halu Oleo, INDONESIA
2Department of Mathematics, Universitas Halu Oleo, INDONESIA
Abstract: - This study aims to explore the heat flow transfer on materials (i.e., homogenous material, particle
material, and sandwich material) by using an open-source simulation. The heat flow occurs due to the
conduction process equation with the 2T model of the source. We use the Finite Element Method (FEM) to
obtain the global heat transfer solution without heat interaction between the walls or layers. The results showed
that each domain has a different temperature value according to the point and time used. So further research is
expected to research other types of heterogeneous materials.
Key-Words: - Finite Element Method, Heat Flow, particle Material, and sandwich material.
Received: May 19, 2022. Revised: November 16, 2022. Accepted: December 27, 2022. Published: February 7, 2023.
1 Introduction
Recently, there has been a tendency to design the
material structure such as composite for a particular
application, where various heterogeneities formed
from particles or sandwich material are also present
due to the manufacturing procedure. Thermal
application in heterogeneous materials, [1] range for
insulation, heat exchangers, and heat sinks had
characterized by the thermal behavior and heat flow
processes [2]. Thermal behavior is essential in
numerous applications, affecting any device's
lifetime and safe operations [3]. The recent
experimental findings [4], [5] suggest that
heterogeneous materials can show thermal behavior
different from Fourier's law at room temperature in
a macroscale object due to parallel heat conduction
channels. The use of approaches two-temperature
model, compared with the experimental data under
certain conditions for physical background and
characteristic behavior, found that the 2T model
offers different insights about the observed heat
conduction phenomenon [6].
The difference in the heat conduction channels
originates in the heterogeneity. For instance, the
homogenous material in a metal sample has good
conductivity properties, while the particle material
in the inclusions is more similar to an insulator.
Assume the modeling that thermal radiation and
heat convection could not contribute to the thermal
behavior [7].
The two-temperature (2T) approach, [6], [7], [8]
however, it has not been tested on heat pulse
experiments so far but has a good background and
could influence, e.g., the thermoelectric conversion
processes [9], [10]. The 2T model has been
successfully used for heat transport in metals under
ultrashort heat pulses, [11] but they are not tested or
not applicable for macroscale solids at room
temperature. Furthermore, many other heat
equations can be found in the literature [12], [13],
[14], but type equations are derived as exceptional
cases.
In the present paper, we focus on heat transfer
and the thermal characterization of heterogeneous
materials, considering only isotropic, heat
conduction channels, and constant material
parameters using two thermal approaches developed
for the materials. We needed the assumption that
the initial and boundary conditions can differ from
the subsystem and that both subsystems receive
different amounts of heat. We also assumed that the
heat capacities are equal, and the model needs
information about the material structure, especially
the constituents. Illustrations or splitting the
computational domain into individual small patches
and finding local solutions continued on these
patches to obtain global heat transfer solutions using
the finite element methods (FEM), [15], [16], [17].
However, we need a good command of matrix
algebra and computer programming to apply a
powerful tool to engineering problems and obtain
valuable solutions.
2 Numerical Formulation
In this study, we suppose that the heterogeneous
material in heat conduction channels (i.e., particle
material and sandwich material) divide into two
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
1
Volume 18, 2023
subsystems that obey the Fourier law. The system is
characterized by two diffusivities, presenting two
distinct characteristic time scales. The heat transfer
model shows three spatial dimensions as
x ; (1)
(2)
where ρ, c, , T, , f (x) are the mass density
(kg/m3), specific heat (J/kg K), thermal conductivity
(W/m K), temperature (K) and heat flux (W/m) and
source term ((W/m2), for the corresponding
subsystem, respectively. Moreover, ∂ it denotes the
partial derivative concerning time (t) or space (x). If
there are no further the source terms f (x)= 0, the
system will be in the balance of internal energy.
The temperature on the boundary of the body is
prescribed by
where T0 is a known function, simply a constant.
Equation (2) is a constitutive, called Fourier's law.
If the heat flux out of the body (perpendicular
to the surface) is determined, then Fourier's law
helps us in deciding the partial derivative of the
temperature concerning the outward normal vector
n: ;
For the perfect insulation in
ГN, the equation becomes a homogeneous Neumann
boundary condition
. Transformation of the
model to the non-dimensional system was using the
finite element methods (FEM) with assumptions
that
,
then equation (1) became
;
N ;
in D (3)
The entire system is characterized by two
diffusivities, presenting two distinct characteristic
time scales, where a situation can occur in two
specific components, such as sandwich material and
particle material. The internal heat generation may
emerge due to chemical reactions or a heat current
flowing through the rod; the left end is kept at a
constant temperature, while at the right boundary,
no heat is allowed, the source term
, the
heat exchange between the subsystems and the
interface area is uniform in the domain of interest.
Assuming that this functional system is also
appropriate for the temperature function T and can
express the existence of parallel heat conduction
channels, the problem may now be written in the
weak formulation and T a solution in the weak
sense. For this purpose, we multiply the differential
in equation (3) with the so-called test function v,
which
and and integrate both sides
over the whole domain Ω = (0, L). Therefore,
.
(4)
If we require the test function v to satisfy the same
Dirichlet boundary condition, as the
temperature T, and remember the homogeneous
Neumann boundary condition
, by using the
Divergence Theorem, equation (4) becomes
(5)
In conduction channels
, so equation (5) becomes
(6)
The first derivative of v appears inside an integral.
Since the integration provides a smoothing effect,
Need to be continuous, so the calculus comes into
play for a differentiable function of the weak
derivative that is equal to the usual result.
If time discretization was using implicit schema and
assuming that
, qk(x)
and
which is
time interval ( , the heat
transfer model show as
(7)
where
is a given function and
An unknown
function.
The numerical solution simulation uses
FreeFem++ software, an open-source software
combined with Gnuplot, to perform a chart. One
must define initial and boundary conditions to solve
these heat conduction models. In this study, we
will restrict the material to three domains (i.e.,
a homogeneous domain, a particle domain, and a
sandwich domain, with a specific boundary surface
to make it easier to simulate. Heterogeneous
materials define as materials with dramatic
heterogeneity in composition from one domain area
to another. The domain sizes could range from
micrometers to millimeters, and the domain
geometry can vary to form very diverse material
systems.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
2
Volume 18, 2023
One must define the initial temperature
distribution and boundary surface for a homogenous
domain material model. The boundary conditions
are defined through the constitutive equations to
preserve physical consistency. In particle material,
the subsystems are composed of large-size particles
scattered in the material. The boundary conditions
do not require special situations, whereas, in
sandwich material with a multi-layer domain, each
layer has different characteristics, which could be a
restrictive property. The field is rectangular on the
x-axis and y-axis, as shown in figure1.
First, we need to perform boundary coordinates in
the x-y plane, which surface one is on
= t
, surface II on
, surface III on
,
and surface IV on
. In
FreeFem++ programming, each coordinate point
will be generating mesh and connecting each line in
a triangle form, for example
int n =5.;
border C1(t=0,1){x = 1+t;y=1;label=1;}; /Surface I
border C2(t=0,1){x = 2; y=1+t;label=2;};/Surface II
border C3(t=0,1){x = 2-t;y=2;label=3;}; /Surface III
border C4(t=0,1){x = 1; y=2-t;label=4;}; /Surface IV
mesh Th = build mesh(C1(n)+C2(n)+C3(n)+C4(n));
plot(Th,ps="mesh.eps");
On the particle domain displayed in figure 1, we
need to perform 14 circles with a diameter of 0.15;
for example, the center point of circle 1 is on x
=
1.1
+
0.075 cos
θ
and y
=
1.2
+
0.075 sin
θ
, circle
two on x
=
1.3
+
0.075 cos
θ
and y
=
1.2
+
0.075sin
θ
, and so forth.
On the other side, for sandwich-structured
composite, the domain comprises four layers of
material consisting of two types of materials with
13 boundary surfaces. On the x-y plane, surface
one is lying on
= t
, or x
=
1
+
t and y
=
1, surface 2 on
= t
, or x
= 2
and y
=
1
+
0.25t, surface
three on
= t
, or x
= 2
and y
=
1.25
+
0.25 t, and so on.
3 Results and Discussion
3.1 Numerical Simulation on Homogeneous
Materials
The system has an initial temperature of 0; on the
top layer at y = 2, the temperature is set at 300.
Otherwise, at the bottom layer at y = 1, the
temperature is set at 20. The interval used for
iteration time is 0.01. The iteration process starts
from t = 0, 0.01, 0.02, up to t = 5. Regarding the
boundary conditions, the heat will not come out
from the left and the suitable material because that
part is the insulin surface which will not penetrate
heat and does not transfer heat, so the temperature
inside changes faster. The illustration of heat flow
in the homogenous material is present in figure 2.
We will consider on three-point where the first
point (1.5,1.8), the second point (1.5,1.5), and the
third point (1.5,1.2) have the same temperature
initial 00C at the t = 0. Furthermore, the
temperature at the first point, the nearest
temperature of 3000C, has changed at all times,
where for 100 times iteration at t = 1, it is 243° C;
for 200 times iteration at t = 2, it is 244° C, and at t
= 3, t = 4, and t = 5 the temperature value is is the
same as when t = 2 which is 244° C. Furthermore,
at the second point, the temperature value has
changed at each time, t = 1 of 159° C, at time of t =
a. Homogen domain
b. Particle Domain
c. Sandwich Domain
Fig. 1: Illustrating the materials: a. Homogen domain, b. Particle domain, and c. Sandwich domain
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
3
Volume 18, 2023
2 of 160°, and at time of t = 3, t = 4, and t = 5 the
value the temperature is the same as when t = 2
which is 160° C. For the third point near the bottom
surface with a temperature of 200 C, the temperature
has changed, where for 100 times iteration at t = 1,
the value was 75° C. At 200 times iteration, t=2, its
temperature was 76° C and stable not change for t =
3 (300 iterations), t = 4 (400 iterations), and t = 5
(500 iterations).
The dramatic temperature change appears for 0
to 200 times iterations and then tends to be stable. It
is clear that along the y-axis, the temperature
changes vertically where at the first point, nearest
the source at 300 0C, the temperature was stable at
244 0C; at the second point, in the middle, the
temperature was stable at 160 0C and the third point,
in the bottom nearest temperature 20 0C, the
temperature was steady at 760C. On the other side,
from figure 2, the shape of the heat flow is flat,
which shows that the x-axis will have the same
temperature at each of the same y-coordinate.
Fig. 2: Temperature distribution on homogeneous
materials at different iteration times ; (a) when t = 0,
(b) when t = 1, (c) when t = 2, (d) when t = 3, (e)
when t = 4, (f) when t = 5
So, we can see for all points, such as points
"1.5,1.8, 1.5,1.5, and 1.5,1.2", the temperature only
changes for maximal 200 times iteration from t=
0.01 to t = 2, then at t = 3, t = 4, and t = 5 tend to be
stable and unchanging as shown by figure 3.
Fig. 3: Temperature graph of the observable
point in the homogeneous material simulation
3.2 Numerical Simulation of Particle
Material
For particle material, we consider on three-point, the
first point (1.4,1.8), the second point (1.4, 1.5), and
the third point (1.4, 1.2), which have the same
temperature initial 00C at the t = 0. After 100
times iteration at t = 1, the temperature at the first
point near the heat source 3000C is 247,010C, in the
middle for the second point is 159,751° C, and for
the third near the heat source (20 0C) is 72.728 °C.
Furthermore, at 200 times iteration t = 2, the first
point temperature is 247.150C, the second point is
159.9970C, and the third point temperature is
72.8570C. Furthermore, in the next iteration, with
t = 3, t = 4, and t = 5, the temperature of the three-
point tends to be stable and does not far differ from
t = 2, as shown in figure 4.
Fig. 4: Temperature distribution on material
particles at different times iteration; (a) when t = 0,
(b) t = 1, (c) t = 2, (d) t = 3, (e) t = 4, and (f) t = 5
There are differences in the distribution of
temperature values at each point reviewed of the
particle material compared to the homogeneous
material. For the y-axis, the temperature changes
a
b
c
d
e
f
a
b
c
d
e
f
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
4
Volume 18, 2023
vertically where at the first point, nearest the source
at 300 0C, the temperature was stable at 247.15 0C;
at the second point, in the middle, the temperature
was stable at 159.997 0C, and the third point, in the
bottom nearest temperature 20 0C, the temperature
was steady at 72.8570C.
For homogeneous material, the shape of the heat
flow is flat. In contrast, for particle material, the
heat flow makes a ripple which shows that the x-
axis will have different temperatures at each of the
same y-coordinate, so they have a different pattern.
The temperature of particle material is at a parallel
point; however, the difference in temperature values
is small enough between points along the x-axis to
cause the appearance of small ripples, such as
shown in figure 4.
The graphic simulation in figure 3 and figure 5
shows the case that the exact coordinates point has
little difference in temperature value for the two
sample materials.
Fig. 5: Temperature graph of the observable
point in particle material simulation
Considering the two coordinate points (1.4,1.8)
and (1.5, 1.8) near the heat source 3000C, only one
line appears because overlain both mean almost the
same temperature values. Likewise, in coordinates
(1.4,1.2) and (1.5,1.2), point near the heat source
200C, only one line appears because of the overline
of the same temperature value. Minor changes to
the temperature values along the x-axis are only
observed when plotting the data for all points and
will appear like ripples, as shown in Figure 4.
3.3 Numerical Simulation on Sandwich
Materials
The results of numerical simulations on sandwich
materials obtained from the FreeeFEM++ are shown
in Figures 6 and 7. For sandwich material, we
consider four points where the first point (1.5,1.85),
the second point (1.5, 1.65), the third point (1.5,
1.4), and the fourth point (1.5, 1.15) have the same
with initial temperature 00C at the t = 0. 100 times
iteration at t = 1, the temperature at the review point
becomes, the first point at 243.818 0C, the second
at 187.692°C, the third at 122,393°C, and the
fourth at 47,906°C. Furthermore, at 200 times
iteration, t=2, the temperature increase slightly from
t =2, the first point temperature at 2440C, the
second point at 187,90C, the third point at
122,660C, and the fourth point at 47,99°C. Then, in
the 300 iterations, t = 3, the first point temperature is
2440C, tends to be the same with t =2, whereas the
second point, the third point, and the fourth point,
the average temperature rises a little, respectively at
1880C, 122,670C, and 48°C. For 400 iteration t = 4
and 500 iteration t = 5, the average temperature of
all points considered tends to be stable. It does not
differ from t=3, as shown in figure 6 and precisely
the d section of the illustration.
Fig. 6: Temperature distribution on sandwich
material at different times; (a) when t = 0, (b) when t
= 1, (c) when t = 2, (d) when t = 3, (e) when t = 4,
(f) when t = 5.
Based on the selected at the three coordinate
points (1.5,1.8), (1.5,1.5), and (1.5,1.2) for the 2T
model, which T1at the top of the material is 300°C,
and T2 on the bottom is 20°C, both homogenous
material and heterogenous material approaches can
model the same phenomenon. However, there are
a
b
a
b
c
d
e
f
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
5
Volume 18, 2023
several differences, such as the shape of heat flux
and the amount of heat received by each point at the
same iteration time. First, with the same initial and
boundary conditions, we needed the assumption
that all material domains receive the same amount
of heat. It appears that from the three points, the
homogeneous domain temperature average of point
1 at t = 1s is 247,0060C, t = 2, t = 3, t = 4, and t = 5
at 247.1480C and the material particle domain at t
= 1 is 247.010C, t = 2 at 247.150C, t = 3s t = 4s and
t = 5 at 247.150C. Successively, the temperature
average of point 2 for homogenous material and the
material particle domain at t = 1 is 159°C and
159.751°C, t = 2s is 160 °C and 159.9972°C, t = 3s
is 160 °C and 159.9973°C when t = 4s is 160 and
159.9974°C, and t = 5s is 160 °C and 159.9975°C.
Point 3, the temperature average for homogenous
material at time t = 1 is 75°C, at time t = 2, t = 3, t =
4, and t = 5 is 76°C. In contrast, for the material
particle domain at time t = 1 is 72.7283°C, at time t
= 2 is 72.8566°C, at time t = 3 is 72.85683°C, at
time t = 4 is 72.85684°C, and at when t = 5 is
72.85685°C.
With the assumption that the heat capacities are
equal, information about the material structure,
especially about the constituents, is necessary for
heat flux to determine the difference in the heat
received due to differences in the composition of the
material. Figure 7a, figure 7b, and 7c show where
the grey line represents the homogenous material,
and the purple line represents the heterogeneous
material. It shows that a model should be predictive
for practical applications. On designs of composite
structures, it is easier to estimate the characteristic
time scales using the average temperature model
that accommodates the Fourier heat equations.
However, creating heterogeneous material would be
challenging as it depends on several factors, such as
the material structure and a specific heat transfer
coecient for a constituent in the future
investigation concerning the presented physical and
mathematical aspects.
The graph in figure 7 shows that each point has a
different temperature level. It is also clear that the
increase in temperature occurs from time t = 0 to t =
1 (iteration 100 times) and towards t = 2 and tends
to stable in more than 300 times iterations (t = 3, 4,
and 5). From the point observable shown in the
graph in fig. 7, the homogenous materials have a
uniform temperature along the x-axis. They are also
higher than the particle material and the sandwich
material.
Fig. 7: Temperature distribution of homogenous,
particle, and sandwich materials at different times
iteration and point coordinates
4 Conclusion
Thermal conductivity of heterogenous and
homogenous materials have been successfully
simulated using open-source software. The heat
flow has several differences, such as the shape of
heat flux and the amount of heat each point receives
at the exact iteration times. For homogeneous
material, the form of the heat flow is flat, while for
particle material, the heat flow makes a ripple, so
they have a different pattern. The heat flow of
sandwich material has distributed horizontally
where the temperature change rate maximum only
occurs from 0 to 200 times iteration and tends to be
the fastest instability. It would be interesting to
extend the analysis of temperature eciency and
distribution in both homogenous and heterogenous
domain material and the usefulness of a practical
system of heat transfer processing
References:
[1] X. Wu and Y. Zhu, "Heterogeneous
materials: a new class of materials with
unprecedented mechanical properties",
Mater. Res. Lett., vol. 5, no. 8, pp. 527–532,
2017.
[2] J. Blackwell, "Microstructure
Characterization," Struct. Form. Polym.
Fibers, pp. 457–520, 2001.
[3] Y. T. Zhu, T. C. Lowe, and T. G. Langdon,
"Performance and applications of
nanostructured materials produced by severe
plastic deformation," vol. 51, pp. 825–830,
2004.
[4] M. Wang, N. Yang, and Z. Y. Guo, "Non-
Fourier heat conductions in nanomaterials,"
J. Appl. Phys., vol. 110, no. 6, 2011.
a
b
c
d
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
6
Volume 18, 2023
[5] R. E. Gonzalez-Narvaez, M. López De Haro,
and F. Vázquez, "Internal Structure and Heat
Conduction in Rigid Solids: A Two-
Temperature Approach," J. Non-Equilibrium
Thermodyn., vol. 47, no. 1, pp. 13–30, 2022.
[6] R. Kovács, A. Fehér, and S. Sobolev, "On
the two-temperature description of
heterogeneous materials," Int. J. Heat Mass
Transf., vol. 194, 2022.
[7] Y. Hua, and B. Cao, "International Journal of
Thermal Sciences Ballistic-diffusive heat
conduction in multiply-constrained
nanostructures," Int. J. Therm. Sci., vol. 101,
pp. XXV–XXV, 2016.
[8] M. Di Domenico, "Nonlocal and nonlinear
effects in hyperbolic heat transfer in a two-
temperature model," Zeitschrift für Angew.
Math. Und Phys., vol. 72, no. 1, pp. 1–15,
2021.
[9] A. A. Abood, K. Amarray, M. Garoum, F.
Pirman, and A. Awang, "An extended hot
plate method for measurement of thermal
conductivity variation with temperature of
building materials An extended hot plate
method for measurement of thermal
conductivity variation with temperature of
building materials," 2018.
[10] H. Shibata, A. Suzuki, and H. Ohta,
"Measurement of Thermal Transport
Properties for Molten Silicate Glasses at
High Temperatures by Means of a Novel
Laser Flash Technique," vol. 46, no. 8, pp.
1877–1881, 2005.
[11] S. L. Sobolev, and W. Dai, "Heat Transport
on Ultrashort Time and Space Scales in
Nanosized Systems: Diffusive or Wave-
like?," no.15(12), pp. 4287, 2022.
[12] A. Lovas, "Guyer-Krumhansl-type heat
conduction at room temperature' 1," pp. 4–7,
2017.
[13] E. P. Scott, M. Avenue, and R. Hall, "The
Question of Thermal Waves in Materials,"
vol. 131, pp. 1–6, 2016.
[14] K. O. V Acs, "Generalized Heat Conduction
In Heat Pulse," pp. 1–15, 2014.
[15] I. Erhunmwun, and U. Ikponmwosa,
“Review on finite element method,” 2018.
[16] Z. Qing, Z. Jiashou, and X. Xiaozhou, "The
Partitioned Mixed Model of Finite Element
Method and Interface Stress Element Method
with Arbitrary Shape of Discrete Block
Element," vol. 2013, 2013.
[17] L. Wang, S. Li, G. Zhang, Z. Ma, and L.
Zhang, "A GPU-Based Parallel Procedure
for Nonlinear Analysis of Complex
Structures Using a Coupled FEM / DEM
Approach," vol. 2013, no. Iii, 2013.
WSEAS TRANSACTIONS on HEAT and MASS TRANSFER
DOI: 10.37394/232012.2023.18.1
Rosliana Eso, Arman
E-ISSN: 2224-3461
7
Volume 18, 2023
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The authors equally contributed in the present
research, at all stages from the formulation of the
problem to the final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflict of Interest
The authors have no conflicts of interest to declare
that are relevant to the content of this article.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
