The Finite Element Method for the Elasticity of Micro-pores Material
ARMAN ARMAN1, ROSLIANA ESO2,*
1Department of Mathematics,
Universitas Halu Oleo,
INDONESIA
2Department of Physics Education,
Universitas Halu Oleo,
INDONESIA
*Corresponding Author
Abstract: -This study aimed to investigate micro-pores' impact on the elastic energy of materials using open-
source simulation. The finite element method analyzed linearly elastic, homogeneous, and isotropic materials.
The domain is modeled as an isotropic elastic material with a constant thickness and is treated as a two-
dimensional rectangular plate. The positions and diameters of the micro-pores were varied. The results indicate
that the position and diameter of the micro-pores within the domain significantly influence the magnitude of
elastic energy. Larger pore diameters correspond to higher values of elastic energy. Additionally, the elastic
energy value of the micro-pores remains relatively constant, though it is smaller when the pores are located at
the ends of the material. This suggests that the elasticity at the end of the material is reduced due to pressure-
induced movement. Furthermore, varying the diameter of the micro-pores leads to more noticeable changes in
material shape, with larger diameters resulting in more pronounced alterations.
Key-Words: - Finite element method, Elastic energy, Isotropis, Rectangular plate, Open-source simulation, and
Micro-porous material.
Received: September 9, 2023. Revised: May 15, 2024. Accepted: June 11, 2024. Published: July 26, 2024.
1 Introduction
Mechanical properties such as elastic modulus of
materials can be affected by variations in pore shape
and distribution, [1]. Micromechanical modeling
mostly focuses on practical changes in the shape of
a material, which is usually determined by various
factors such as ultimate strength, Young's modulus,
and Poisson's ratio, [2]. Micromechanical modeling
for predicting the linearly elastic limit properties
also extends to effective stress-strain behavior
beyond the tensile flow stress in anisotropic sheet
samples, [3]. Apart from the material's length,
width, and thickness, the micro-pore effect is also
part of the geometric effect, which contributes
significantly to the change in the shape of the
material. It should be noted that the microstructure
of porous materials in practical existence consists of
pores distributed irregularly in various sizes and
shapes. In particular, the distribution mode of these
pores also clearly affects the overall elasticity of the
material, [4]. The difference in the arrangement and
distribution of the pore structure plays a vital role in
optimizing the pore topology and achieving the
desired effective material properties or tailoring
specific textures, [5].
Recent experimental findings show that stress
transfer increases with increasing porosity [6],
suggesting that numerical modeling only considers
idealized pore arrangements. Hence, the stress
induced by the external excitations can be
transferred more smoothly in this microstructure but
provides physical insight into how microscopic
porous shapes affect macroscopic mechanical
responses through systematic analysis. This feature
may give novel ideas for designing porous
structures or bio-implants with specific
requirements, [7]. The use of approaches on a
uniaxial compression of micro-porous rubber [8],
compared with the experimental data for the stress–
stretch and volume change–stretch data of
elastomers in uniaxial tension [9] under certain
conditions for physical background and
characteristic behavior found that the constitutive
response of an elastic material offers different
insights about the observed elasticity micro-pores
phenomenon.
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
72
Volume 21, 2024
In this study, we focus on the shape of the
material and the elastic properties of the micro-
porous material, which are only isotropic, with
uniaxial stress tensors considered, [10]. The initial
and boundary conditions could be the same for
subsystems, which requires the assumption that both
subsystems receive the exact amounts of pressure.
Numerical techniques can be conveniently utilized
to analyze porous structural materials with
vacancies, [11]. We also assume that the stress
tensors are equal and that the model requires
information about the material structure, especially
its components. We can find local solutions
following these patches to divide the computational
region into illustrations and individual small
patches, obtaining a deformation solution using the
Finite Element Method (FEM), [12], [13], [14],
[15]. However, one needs proficiency in matrix
algebra and good computer commands to apply
powerful tools to engineering problems and obtain
valuable solutions.
2 Numerical Formulation
In this paper, we denote the gradient and divergence
concerning x = (x1,x2)T 2 by div,
respectively. We denote the Lebesgue space on Ω
by L2(Ω) and the 2-valued Lebesgue space by L2
(Ω; 2). We also often use the Sobolev space
defined by H1(Ω) = {u L2(Ω); u L2(Ω; 2)},
and its trace space on the boundary ΓD, u denotes the
unknown displacement field u(x)=((u1(x),u2(x))T
2. Then, the strain tensor e[u](x) is defined by
and the stress
tensor is denoted by σ[u](x) = (σij[u](x))
and
defined as σ[u] = Ce[u], the elasticity tensor
C(x)=(cijkl(x)) asumed in symmetry condition, so the
overall constitutive equation follows Hooke's law
(1)
If the representative volume element is assumed
homogeneous, isotropic and linearly elastic, then
or anisotropic elastic tensor expressed as
Here, λ
and μ are Lame constant and δij is a delta cronecker,
so eq (1) become:
(2)
For each elastic fracture mechanics, a small
displacement represents deformation[16] where the
superscript T denotes the transposition, equation 2
be written as:
(3)
Therefore, the strain component can be
expressed as a function of the stress component in
eq.3 by combining:
If the stress tensor is uniaxial,
σ11 = σ, and σ22 = σ33 = σ12 = σ23 = σ31 = 0, the
corresponding strain comp e12 = e23 = e31 = 0, then
Young's modulus and Poisson's ratio defined as:
the Lame constant form:
(4)
Assume that an elastic body occupies the
domain volume Ω and let ω be a subdomain of Ω
with boundary and normal vector n, concerning
the elastic equilibrium equation defined as [17],
=0 (5)
Then, in purely elastic materials assuming, by
using the divergence theorem the linear elasticity
equation in eq (5) can generally be expressed in
terms of partial differential equations governing
deformation and stress a slight displacement
gradient (6)
u = g (x), , ,
,
Then, ;
,
where f(x) is a given external load perpendicular to
the plate on Ω, g(x) is a given anti-plane
displacement on ΓD, and q(x) is a given boundary
load in the x3-direction on ΓN . The outward normal
derivative on the boundary of Ω is denoted by ∂/∂n .
Transformation of the model to the non-
dimensional system was using the Finite Element
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
73
Volume 21, 2024
Methods (FEM) with assumptions,
,
then,
Assuming that this functional system is also
appropriate for the pore diameter and can express
the existence of parallel channels with variations in
the location/position of the material pores, the
problem may now be written in the weak
formulation, [18]. For this purpose, we multiply the
differential in equation (6) with the so-called test
function v, which,
and if g , and integrate both sides
over the whole domain Ω , therefore,
. (7)
If we require the test function v to satisfy the
same boundary condition, , by using the
Divergence Theorem, equation (6) becomes
(8)
The displacement 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.
(9)
where
is an unknown
function, and
is a
given function.
At time t [0, T], we denote the weak solution
to (8) by u(t). Then, the weak solution u(t) is
characterized by the following variational principle,
i.e., it is obtained as a unique minimizer of the
following elastic energy including the body and
surface forces under a suitable boundary condition
(10)
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 with the
plate is supposed to be an isotropic elastic material
with a constant thickness and is treated as a two
dimensional rectangular on the x-axis and y-axis in
Ω field with the border in cartesian domain
border A(t=a1,a2) {x = t; y = b1;
border B(t=b1,b2) {x = a2; y = t;
border C(t=a2,a1) {x = t; y = b2;
border D(t=b2,b1) {x = a1; y = t;
The micropores domain restrict the border L
(t=0, 2*pi) and the domain position can vary to
very diverse diameter systems showed in Figure 1.
Fig. 1: Numerical domain on the cartesian plane
The first design is to determine different pore
positions (h), namely h1= (0.2) μm h2= (0.4,) μm,
h3= (0.6) μm and h4= (0.8) μm. After that, the pore
diameter (d) was varied at each different pore
position from the value range 0.01 μm ≤d≤ 0.09 μm.
The second design used four materials with different
pore diameters (d): 0.02 μm, 0.04 μm, 0.06 μm, and
d4 = 0.08 μm. The pore position (h) was varied for
each pore diameter. The parameters Young's
modulus (λ) and Poisson's ratio (μ) [19], are written
with the following code: real mu = 7.617e10; real
lambda = 9.69e10; real rho = 7700; real E = 1.84e6;
and real nu = 0. . Whereas the elastic linear
equations are written as solvelame([u1,u2],[v1,v2])
=
int2d(Th)(lambda*div(u1,u2)*div(v1,v2)
+ 2.*mu*(eps(u1,u2)'*eps(v1,v2)))
- int2d(Th)(gravity*v2)
- int1d(Th,3)(23.75*alpha*
-sin(0.3*pi*t)*v2)+on(4, u1=0, u2=0)+on(2, u2=0).
3 Results and Discussion
3.1 Numerical Simulation of Displacement
on a Variety of Micropore Positions
Because the pores in porous materials are frequently
irregularly distributed and formed in the solid phase,
it is challenging to ascertain the total elastic
characteristics using analytical or experimental
means. For this reason, the numerical analysis
becomes essential. A representative volume element
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
74
Volume 21, 2024
(R) is established if the porous material is
considered homogenous at the macroscopic scale.
The local response of the absorbent material, as
determined by micromechanics, may be used in
conjunction with the particular boundary condition
to derive the overall moduli of porous materials.
Below are the results of a displacement simulation
in location (H = 0,1m and 0,2m) where the micro-
pore diameter is varied. The results are shown in
Figure 2, Figure 3, Figure 4 and Figure 5.
d = 0,02m
H=0,1m
H=0,2m
Fig. 2: Displacement simulation at micro-pore
diameter of 0.2 μm
Fig. 3: Displacement simulation at micro-pore
diameter of 0.4 μm
In the mixed position configuration, small
circular pores oriented along the load direction play
an important role compared with large ones.
Discrete plastic deformation bands are sharper and
essentially bypass these pores. The earlier the
damage starts and, therefore, the lower elasticity can
be expected.
d = 0,06m
H= 0,1m
H=0,2m
Fig. 4: Displacement simulation at micro-pore
diameter of 0.6 μm
d = 0,08m
H= 0,1m
H=0,2m
Fig. 5: Displacement simulation at micro-pore
diameter of 0.8 μm
3.2 Numerical Simulation of Displacement
due to Change in Micropore Diameter
Displacement simulations with various micropore
diameters determined by the design location
micropores of H1=(0.2)µm, H2=(0.4) µm,
H3=(0.6)µm, and H4=(0.8) µm are provided. The
diameter (d) of the micropores was then adjusted for
each separate pore size in the range of 0.01 µm to
0.09 µm. The results are shown in Figure 6, Figure
7, Figure 8 and Figure 9.
d= 0,04m
H= 0,1m
H=0,2m
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
75
Volume 21, 2024
Fig. 6: Simulated displacement of micropores at
H1=0.1 µm, diameter varied from 0.01 to 0.09 µm
H2 =0,4 μm
d=0,02 m
d=0,05 m
d=0,09m
Fig. 7: Simulated displacement of micropores at
H2=0.4 µm, diameter varied from 0.01 to 0.09 µm
A rectangular material with micropores is seen
in the above photograph. The top and lower ends of
the micropores in the material broaden as the
micropores' diameter varies; the higher the
micropores' diameter, the larger the shape of the
material changes. Based on the strain energy
function of microporous materials under uniaxial
compression, a volumetric strain energy function is
produced, taking into account the mathematical link
between strain energy and the impact of porosity on
the relationship between stress and strain. The
simulation results for this area may be used to
calculate the energy level sufficient to achieve the
elastic value, and a graphical representation of the
results is required. The graph in Figure 10 shows the
elastic energy level as a consequence of the
microporous region's pore location and pore size.
H3 = 0,6 μm
d= 0,02m
d=0,05m
d=0,09m
Fig. 8: Simulated displacement of micropores at
H3=0.6 µm, diameter varied from 0.01 to 0.09 µm
H4 = 0,8m
d= 0,02m
d=0,05m
d=0,09m
Fig. 9: Simulated displacement of micropores at
H4=0.8 µm, diameter varied from 0.01 to 0.09 µm
A rectangular material with micropores is seen
in the above photograph. The top and lower ends of
the micropores in the material broaden as the
micropores' diameter varies; the higher the
micropores' diameter, the larger the shape of the
H1 = 0,2 μm
d=0,02 μm
d=0,05m
d=0,09m
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
76
Volume 21, 2024
material changes. Based on the strain energy
function of microporous materials under uniaxial
compression, a volumetric strain energy function is
produced, taking into account the mathematical link
between strain energy and the impact of porosity on
the relationship between stress and strain. The
simulation results for this area may be used to
calculate the energy level sufficient to achieve the
elastic value, and a graphical representation of the
results is required. The graph in Figure 10 shows the
elastic energy level as a consequence of the
microporous region's pore location and pore size.
(a)
(b)
Fig. 10: Elastic energy of micropores: (a) varied in
pore displacement and (b) varied in pore diameter
Elasticity has been shown to increase with the
diameter of the micropores, indicating that
differences in the size of the micropores affect the
elastic value of the material. The elastic value is
almost constant when the micropores shift to
different positions, but the value decreases when the
micropores are at both ends of the material,
indicating that the elasticity at the ends of the
rectangular material is smaller than when the
micropores are located in the middle of the content.
The greater the pore diameter, the greater the
elasticity energy value. However, with variations in
position, the elasticity value of the micropores is
nearly constant, even though the value is small
when the micropores are at the ends of the material,
indicating that the elasticity at the ends of the
rectangular material is smaller than when the
micropores are in the middle. However, creating
elasticity material would be challenging as it
depends on several factors, such as the microporous
position and the microporous diameter for a
constituent in future investigations concerning the
presented physical and mathematical aspects.
4 Conclusion
Rectangular domains with micropores have been
successfully simulated using open-source software.
The presence of micropores affects the ability of the
material to change its shape. The greater the pore
diameter, the greater the elasticity energy value.
Moreover, the elasticity value of the micropores is
nearly constant, even though the value is small
when the micropores are at the ends of the material,
indicating that the elasticity at the ends of the
rectangular material is smaller than when the
micropores are in the middle. It will close due to the
movement of the material induced by the pressure.
The larger the diameter of the microporous, the
more the shape of the material changes, that is, the
wider the top and bottom surfaces of the material. It
would be interesting to extend the analysis of the
diameter of the microporous and its distribution in
homogenous domain material and the usefulness of
a practical system of the material’s shape change
processing.
References:
[1] Y. L. Shen and N. Fathi, “Numerical study of
elastic-plastic behaviour of pore-containing
materials: effects of pore arrangement,” Int. J.
Theor. Appl. Multiscale Mech., vol. 3, no. 4,
p. 262, 2021. DOI:
10.1504/ijtamm.2021.120795.
[2] Ochiai, Shojiro, Satoshi Nakano, Yuya
Fukazawa, Mohamed Shehata Aly, Hiroshi
Okuda, Komei Kato, Takeshi Isobe, Koichi
Kita, and Keiichi Honma. 2010. “Change of
Young’s Modulus with Increasing Applied
Tensile Strain in Open Cell Nickel and
Copper Foams.” Materials Transactions, vol.
51 (5), pp 925–932.
https://doi.org/10.2320/matertrans.M2009384.
[3] F. Barlat, H. Aretz, J. W. Yoon, M. E.
Karabin, J. C. Brem, and R. E. Dick, “Linear
transfomation-based anisotropic yield
functions,” Int. J. Plast., vol. 21, no. 5, pp.
1009–1039, 2005. DOI:
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
77
Volume 21, 2024
10.1016/j.ijplas.2004.06.004.
[4] A. Srivastava, S. Osovski, and A. Needleman,
“Engineering the crack path by controlling the
microstructure,” J. Mech. Phys. Solids, vol.
100, August, pp. 1-20, 2017. DOI:
10.1016/j.jmps.2016.12.006.
[5] O. Cazacu and J. A. Rodríguez-martínez,
“Effects of plastic anisotropy on localization
in orthotropic materials: New explicit
expressions for the orientation of localization
bands in flat specimens subjected to uniaxial
tension,” J. Mech. Phys. Solids, vol. 126, pp.
272–284, 2019.
https://doi.org/10.1016/j.jmps.2019.03.002.
[6] T. Hung, D. Duc, T. Tu, L. Dung, and C. Bao,
“Static and dynamic behaviour of sandwich
beams with porous core: Experiment and
moving least squares mesh-free analysis,”
Journal of Science and Technology in Civil
Engineering (JSTCE) - HUCE., vol. 18(1), pp.
39–54 2024. http://dx.doi.org/
10.31814/stce.huce2024-18(1)-04.
[7] H. Li, S. Dong, J. Liu, Y. Yu, and M. Wu,
“Finite Element Modeling of Porous
Microstructures With Random Holes of
Different-Shapes and -Sizes to Predict Their
Effective Elastic Behavior,” appl. sci., vol. 9,
pp. 1–20, 2019. Doi: 10.3390/app9214536.
[8] N. Li, H. Xu, L. Hu, X. Wang, and L. Lu,
“Elastic constitutive model and compression
experiment based on microporous silica
composites,” Mech. Adv. Mater. Struct., vol.
July, pp. 1–5, 2020.
https://doi.org/10.1080/15376494.2020.17803
53.
[9] Moerman, Kevin M. Fereidoonnezhad,
Behrooz McGarry, J. Patrick, “Novel
hyperelastic models for large volumetric
deformations,” International Journal of Solids
and Structures, vol. 193-194, pp. 474–
491,2020.
https://doi.org/10.1016/j.ijsolstr.2020.01.019.
[10] Y. Yu and W. Huang, “Selection of
regularization parameter in the Ambrosio-
Tortorelli approximation of the Mumford-
Shah functional for image segmentation,”
Numer. Math., vol. 11, no. 2, pp. 211–234,
2018. DOI: 10.4208/nmtma.OA-2017-0074.
[11] M. Repka, V. Sládek, and J. Sládek,
“Numerical Analysis of Poro-elastic Materials
Described by the Micro-dilatation Theory,”
Procedia Eng., vol. 190, pp. 248–254, 2017.
DOI: 10.1016/j.proeng.2017.05.334.
[12] M. S. H. Mojumder, M. N. Haque, and M. J.
Alam, “Efficient Finite Difference Methods
for the Numerical Analysis of One-
Dimensional Heat Equation,” J. Appl. Math.
Phys., vol. 11, no. 10, pp. 3099–3123, 2023.
DOI : 10.4236/jamp.2023.1110204.
[13] Erhunmwun, Iredia, and Usiosefe
Ikponmwosa,“Review on Finite Element
Method.” J. Appl.Sci.Environ.Manage , vol.
21, no. 5, pp. 999-1002, 2018.
https://doi.org/10.4314/jasem.v21i5.30.
[14] 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,” Mathematical Problems in
Engineering, vol. September, pp. 1-6, 2013.
http://dx.doi.org/10.1155/2013/950696.
[15] 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,”
Mathematical Problems in Engineering, vol.
13, pp. 1-15, 2013.
http://dx.doi.org/10.1155/2013/618980.
[16] M. Kimura, H. Notsu, Y. Tanaka, and H.
Yamamoto, “The Gradient Flow Structure of
an Extended Maxwell Viscoelastic Model and
a Structure-Preserving Finite Element
Scheme,” J. Sci. Comput., vol. 78, no. 2, pp.
1111–1131, 2019.
https://doi.org/10.1007/s10915-018-0799-2.
[17] M. Kimura, T. Takaishi, S. Alfat, T. Nakano,
and Y. Tanaka, “Irreversible phase field
models for crack growth in industrial
applications : thermal stress , viscoelasticity ,
hydrogen embrittlement,” SN Appl. Sci., vol.
June, 2021. DOI: 10.1007/s42452-021-04593-
6.
[18] J.Sladek, V. Sladek, M.Reka, P.L. Bishay,
“Static and dynamic behavior of porous
elastic materials based on micro-dilatation
theory: A numerical study using the MLPG
method,” International Journal of Solids and
Structures, vol. 96, pp. 126–135, 2016.
https://doi.org/10.1016/j.ijsolstr.2016.06.016.
[19] Bischoff, J.E., Arruda, E.M., Grosh, K., “A
new constitutive model for the compressibility
of elastomers at finite deformations,” Rubber
Chem. Technol, vol. 74 (4), pp.541–559.
2001. DOI: 10.5254/1.3544956.
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
78
Volume 21, 2024
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The authors equally contributed in the 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.
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
WSEAS TRANSACTIONS on ADVANCES in ENGINEERING EDUCATION
DOI: 10.37394/232010.2024.21.10
Arman Arman, Rosliana Eso
E-ISSN: 2224-3410
79
Volume 21, 2024
