Micro-Structure Modelling and Electrical Properties Analysis of PZT Matrix Ferroelectric Composites

PZT matrix ferroelectric composite is an important research topic in material science because of its many practical, industrial, and scientific applications. Materials with high dielectric permittivity are used to manufacture electronic devices, particularly capacitors and dynamic random access memory (DRAM). Therefore, the development of reliable and efficient micro models to be utilized in analyzing electrical properties can be of great value in accelerating research in this field. In this paper, a 3D microstructure model for PZT matrix ferroelectric composites has been developed and adopted the finite element method (FEM) to calculate the dielectric constant. The microscopy parameters of developed microstructure model are acquired based on the real composites from X-ray (micro-) diffraction and stereological method. The dielectric constant of different volume ratios of PZT matrix ferroelectric composites can be calculated by accurately controlling the volume of Ferrite particles. At the point of validation, the proposed approach makes visual and numeric comparisons between the morphology of the real microstructure and the model generated by the proposed technique. The simulation results by our method was essentially in agreement with experimental results in other literature. Simulation Experimental results also demonstrate that the dielectric constant of PZT matrix ferroelectric composites is significantly changed while the volume ratio of high dielectric phase particles was below 20%. PZT matrix ferroelectric composites Consequently, this method can be easily extended to composites preparation.


Introduction
A strong dependence between macroscopic properties and micro-structures occurring in polycrystalline materials means that these materials are intensely investigated in various scientific domains [1][2][3]. Piezoelectric ceramic (PZT) has the advantage of promising and robust ferroelectric and piezoelectric properties and it has become extensively applied in electronic industries, such as switching devices, super-capacitors with high energy density characteristics, pyroelectric elements and sensors, micro-electromechanical transducers, piezoelectric sensors and actuators. There are some reports for improving the dielectric properties of the PZT by doping or making composites. An improvement in the properties of PZT has been reported by doping the series of materials such as Nb, Fe, and Mn [4,5]. If the two materials can be combined into dual-phase composites that have the advantages of two materials, the composites have great importance in the potential engineering application point of view where they are used in many electrical devices [6][7][8][9].
However, the main disadvantage of the techniques based on sample preparation and detection results is the complex and time-consuming post-processing required to analyze the amount of data produced. Relying upon the destructive data acquisition methods cannot furbish simulations of the physical response that can be compared against experiments since the specific sample was destroyed during data acquisition on the microstructure itself. In recent years, the introduction of computational simulation techniques has reduced the computational cost and time of the procedure [10][11][12][13][14]. Differing from a physical experiment, the method provides simulation using random numbers to conduct many virtual experiments on a computer and then analyzes the results of these experiments on a statistical basis to draw conclusions. Random close packing of spherical (RCPS), which is confirmed to be an effective technique, has been proven to be widely used in computer simulation [15,16]. Therefore, the model of PZT matrix ferroelectric composites is built by RCPS.
It is well-known that numerical simulations based on the finite element method (FEM) are effective techniques for performing the analysis of properties for composites [17][18][19][20][21], such as Guo, XG, et al., who established a molecular dynamics (MD) simulation model of high chromium alloy based on alloy structure by using atomic random replacement and two bubble algorithms. Thus, the properties of the high chromium alloy model are improved, and it becomes more suitable for precision machining [22]. Sheikholeslami, M uses the finite element method to simulate the unstable process of a complex shape energy storage box. It is concluded that a higher heat transfer rate can be obtained by using nano-reinforced phase change material (NEPCM) instead of pure phase change material (PCM) [23]. The hereby work is focused on polycrystals particles of PZT matrix ferroelectric composites. The numerical method, allowing to calculating dielectric constant on the basis of simulations performed for synthetic microstructures, is proposed for these materials.
In particular, the ability to accurately reproduce microstructural features is fundamental for the valid evaluation of the electrical properties' finite element method of PZT matrix ferroelectric composites. To illustrate the effectiveness of the proposed method, the approach presented in this paper relies solely upon the data collected from the 2D scanning electron microscope (SEM) images of polycrystalline microstructures. The effects of the volume ratio of high dielectric phase particles on the electrical properties of dual-phase composites are investigated. The capability of the simulation results is verified with the results of the experiment available in other literature.

The Overall Approach
All the programs in this paper are implemented through MATLAB and TetGen. To evaluate PZT matrix ferroelectric composites' dielectric constant, the work is divided into four parts, including generating 3D random close packing of spheres, generating grain structure, automatic mesh subdivision and finite unit calculation. The macroscopic dielectric constant of PZT matrix ferroelectric composites is mainly determined by their microstructure and the volume ratios of dual-phase materials. Therefore, the first step is setting the inputs, according to X-ray (micro-) diffraction microscopy provides 3D microstructural phase representations, which contain the number of phases. In addition, the grains of polycrystalline composites have been shown to obey lognormal distribution [24]. In previous work, an effective computer numerical method was established to determine the continuous numerical empirical probability function of the 3D distribution of spherical two-phase particles [25]. Sequentially, the grain model after edge cutting is obtained by counting the high-low dielectric phase grain element distribution of the RCPS model and adjusting boundary conditions.
After that, the grain model is an adaptive subdivision by TetGen software. The data is loaded to MATLAB to process. The characteristics of finite element units are determined by the processed data.
Then, the PZT matrix ferroelectric composites' dielectric constant is carried out to evaluate adopted the finite element method. According to the minimum energy rule, the macroscopic dielectric constant is obtained by reassembling the finite element units and calculating. Finally, the macroscopic dielectric constant under different ratios of dielectric constant can be obtained.
By changing the ratios of the second phase particles in ferrite/PZT dual-phase composite, its electrical properties can be controlled which shows satisfactory behavior on different applications.

Generating 3D Random Close Packing of Spheres
In order to establish a three-dimensional microscopic model and ensure the agreement between the model and the PZT matrix ferroelectric composites, the random close packing of spheres (RCPS) is chosen to use to form the microscopic basal model. The RCPS model, a basic geometric model, has provided great help in many physical and engineering problems and played a guiding role in the field of materials [26][27][28]. The RCPS method is a modified rearrangement algorithm, and there is another method called the sequential generation method. The electrical properties of the simulation model vary with the method of the duty ratio per unit volume, the radius of the generated sphere and distribution pattern [29,30].
Thus, according to X-ray (micro-) diffraction, microscopy provides 3D microstructural phase representations and there is no intermediate or interfacial phase in PZT matrix ferroelectric composites [3,28,31]. Furthermore, by controlling the above-mentioned conditions, the distribution of our simulation model is in agreement with real PZT matrix ferroelectric composites. The particle size distribution generated is shown in Figure 1.

Generating 3D Random Close Packing of Spheres
In order to establish a three-dimensional microscopic model and ensure the agreement between the model and the PZT matrix ferroelectric composites, the random close packing of spheres (RCPS) is chosen to use to form the microscopic basal model. The RCPS model, a basic geometric model, has provided great help in many physical and engineering problems and played a guiding role in the field of materials [26][27][28]. The RCPS method is a modified rearrangement algorithm, and there is another method called the sequential generation method. The electrical properties of the simulation model vary with the method of the duty ratio per unit volume, the radius of the generated sphere and distribution pattern [29,30].
Thus, according to X-ray (micro-) diffraction, microscopy provides 3D microstructural phase representations and there is no intermediate or interfacial phase in PZT matrix ferroelectric composites [3,28,31]. Furthermore, by controlling the above-mentioned conditions, the distribution of our simulation model is in agreement with real PZT matrix ferroelectric composites. The particle size distribution generated is shown in Figure 1. The generated RCPS model is shown in Figure 2.

Generating Grain Structure
After getting the appropriate RCPS model, the Laguerre Voronoi method is used to generate a grain model whose volume fraction could be controlled accurately. To form a grain model without

Generating 3D Random Close Packing of Spheres
In order to establish a three-dimensional microscopic model and ensure the agreement between the model and the PZT matrix ferroelectric composites, the random close packing of spheres (RCPS) is chosen to use to form the microscopic basal model. The RCPS model, a basic geometric model, has provided great help in many physical and engineering problems and played a guiding role in the field of materials [26][27][28]. The RCPS method is a modified rearrangement algorithm, and there is another method called the sequential generation method. The electrical properties of the simulation model vary with the method of the duty ratio per unit volume, the radius of the generated sphere and distribution pattern [29,30].
Thus, according to X-ray (micro-) diffraction, microscopy provides 3D microstructural phase representations and there is no intermediate or interfacial phase in PZT matrix ferroelectric composites [3,28,31]. Furthermore, by controlling the above-mentioned conditions, the distribution of our simulation model is in agreement with real PZT matrix ferroelectric composites. The particle size distribution generated is shown in Figure 1. The generated RCPS model is shown in Figure 2.

Generating Grain Structure
After getting the appropriate RCPS model, the Laguerre Voronoi method is used to generate a grain model whose volume fraction could be controlled accurately. To form a grain model without

Generating Grain Structure
After getting the appropriate RCPS model, the Laguerre Voronoi method is used to generate a grain model whose volume fraction could be controlled accurately. To form a grain model without voids, the Laguerre Voronoi method is adopted to distribute the voids of the RCPS model to each grain. The grain model completely inherits the phase volume ratios of the RCPS model.
The irregular boundary can affect the adaptive mesh subdivision and change the volume fraction of the second-phase particle in this model. In order to eliminate the boundary influence of the 3D microstructure model, it is necessary to cut the edge of the model for forming the model of ferrite/PZT dual-phase composite after getting the appropriate 3D microstructure model.
The count method is as follows: First, the grains are divided into two types: Pa phase grain and Pb phase grain. Next, the total volume of two kinds of grains needs to calculate, called Vt. The total volume of Pa phase, called Va. The total volume of Pb phase, called Vb.
If the ideal volume ratio of high to low dielectric phase grains is Vhigh: Vlow, the total volume of Pa phase grains at this time should be: At this point, the total volume of Pb phase particles should be: By calculating, it can ensure the grains have the same volume ratio as well as RCPS models. The space gap between spheres in RCPS, which is filled when the model of PZT matrix ferroelectric composites formed. In order to make sure of the validity of this model, it's necessary to compare real images of SEM with the model in Figure 3. At this time, the model has begun to take shape, and the grain structure of PZT matrix ferroelectric composites can be seen intuitively at different volume ratios. voids, the Laguerre Voronoi method is adopted to distribute the voids of the RCPS model to each grain. The grain model completely inherits the phase volume ratios of the RCPS model. The irregular boundary can affect the adaptive mesh subdivision and change the volume fraction of the second-phase particle in this model. In order to eliminate the boundary influence of the 3D microstructure model, it is necessary to cut the edge of the model for forming the model of ferrite/PZT dual-phase composite after getting the appropriate 3D microstructure model.
The count method is as follows: First, the grains are divided into two types: Pa phase grain and Pb phase grain. Next, the total volume of two kinds of grains needs to calculate, called Vt. The total volume of Pa phase, called Va. The total volume of Pb phase, called Vb.
If the ideal volume ratio of high to low dielectric phase grains is Vhigh: Vlow, the total volume of Pa phase grains at this time should be: At this point, the total volume of Pb phase particles should be: By calculating, it can ensure the grains have the same volume ratio as well as RCPS models. The space gap between spheres in RCPS, which is filled when the model of PZT matrix ferroelectric composites formed. In order to make sure of the validity of this model, it's necessary to compare real images of SEM with the model in Figure 3. At this time, the model has begun to take shape, and the grain structure of PZT matrix ferroelectric composites can be seen intuitively at different volume ratios.

Automatic Mesh Subdivision
The model established is a three-dimensional model, and the basic unit of three-dimensional division can be divided into three-dimensional tetrahedral units, pentahedral units, hexahedral units, and so on [32][33][34][35]. The grain model established has many irregular polyhedral bodies and many sharp angles before segmentation, so three-dimensional tetrahedral units is chosen to divide them. In order to facilitate subsequent calculation, TetGen software is chosen to achieve the adaptive subdivision of the model. Figure 4 shows that after subdivision of different ferrites volume proportion, in which red particles are for ceramic matrix and green particles represent ferrites. It's clearly found that all basic units are tetrahedron.

Automatic Mesh Subdivision
The model established is a three-dimensional model, and the basic unit of three-dimensional division can be divided into three-dimensional tetrahedral units, pentahedral units, hexahedral units, and so on [32][33][34][35]. The grain model established has many irregular polyhedral bodies and many sharp angles before segmentation, so three-dimensional tetrahedral units is chosen to divide them. In order to facilitate subsequent calculation, TetGen software is chosen to achieve the adaptive subdivision of the model. Figure 4 shows that after subdivision of different ferrites volume proportion, in which red particles are for ceramic matrix and green particles represent ferrites. It's clearly found that all basic units are tetrahedron.

Finite Unit Calculation
The using basic unit is a three-dimensional tetrahedral unit, as Figure 5 shows, so the finite element adopted is also a four-nodes tetrahedral unit. For the convenience of description, it's better to choose the center of the tetrahedral shape to locate the origin of the local coordinate. The nodes and coordinates of four corners are recorded as Angle1 (x1, y1, z1), Angle2 (x2, y2, z2), Angle3 (x3, y3, z3), Angle4 (x4, y4, z4).

Finite Unit Calculation
The using basic unit is a three-dimensional tetrahedral unit, as Figure 5 shows, so the finite element adopted is also a four-nodes tetrahedral unit. For the convenience of description, it's better to choose the center of the tetrahedral shape to locate the origin of the local coordinate. The nodes and coordinates of four corners are recorded as Angle1 (x1, y1, z1), Angle2 (x2, y2, z2), Angle3 (x3, y3, z3), Angle4 (x4, y4, z4).

Automatic Mesh Subdivision
The model established is a three-dimensional model, and the basic unit of three-dimensional division can be divided into three-dimensional tetrahedral units, pentahedral units, hexahedral units, and so on [32][33][34][35]. The grain model established has many irregular polyhedral bodies and many sharp angles before segmentation, so three-dimensional tetrahedral units is chosen to divide them. In order to facilitate subsequent calculation, TetGen software is chosen to achieve the adaptive subdivision of the model. Figure 4 shows that after subdivision of different ferrites volume proportion, in which red particles are for ceramic matrix and green particles represent ferrites. It's clearly found that all basic units are tetrahedron.

Finite Unit Calculation
The using basic unit is a three-dimensional tetrahedral unit, as Figure 5 shows, so the finite element adopted is also a four-nodes tetrahedral unit. For the convenience of description, it's better to choose the center of the tetrahedral shape to locate the origin of the local coordinate. The nodes and coordinates of four corners are recorded as Angle1 (x1, y1, z1), Angle2 (x2, y2, z2), Angle3 (x3, y3, z3), Angle4 (x4, y4, z4).  There are two angular nodes on each side. The potential function is linearly distributed on each side, and the tetrahedron of the function also must be linear on each face. In summary, here is the following potential function which is approximated by interpolation functions: By transforming it into matrix form and using Cramer's Rule, A 1 , A 2 , A 3 , A 4 in the potential function can be obtained.
At this point, the coordinates of each node on the tetrahedral element could be calculated and A 1 , A 2 , A 3 , A 4 can be obtained with the potential function. In other words, all the potentials in a tetrahedral element can be obtained by the potentials of four nodes, so the next step is to obtain the potential values of all unknown nodes.
We need to define the characteristics of the units, this paper uses the energy minimum method, so the first calculation is the energy in each unit.
According to the energy of Poynting theorem, in a three-dimensional static electric field, the energy of the electric field is as follows: where ε is the dielectric constant of the electric field acting area, E is the electric field strength of the area, and the relation with the electric potential U is as follows: where S is a differential matrix for the vertices of a tetrahedron. After obtaining the energy of each tetrahedral unit and assembling all the units together, the total energy of the established model can be obtained, denoted as W. In this paper, the description of the total energy W of the model is the total energy of solving the regional electric field, which can be regarded as the sum matrix of each unit matrix-that is to say, the U matrix and S matrix described in Equation (8). U matrix is the potential matrix, which is composed of the node potential of all tetrahedral units, and the S matrix is composed of the S e matrix of each unit in the model.
As shown in Figure 6, the small model consists of three tetrahedral units. For the tetrahedral Unit A, it consists of four nodes: 1, 2, 4, 5; for the tetrahedral Unit B, it includes four nodes: 1, 2, 3, 4; for the tetrahedral Units C, it includes four nodes: 1, 2, 3, 6. Then, based on the fixed-solution condition, the linear equations of the potential of the model nodes can be built. The solution conditions are set herein: Materials 2020, 13, x FOR PEER REVIEW 6 of 11 There are two angular nodes on each side. The potential function is linearly distributed on each side, and the tetrahedron of the function also must be linear on each face. In summary, here is the following potential function which is approximated by interpolation functions: U = 1 + A 2 x + A 3 y + A 4 z. (3) By transforming it into matrix form and using Cramer's Rule, A 1 , A 2 , A 3 , A 4 in the potential function can be obtained.
At this point, the coordinates of each node on the tetrahedral element could be calculated and A 1 , A 2 , A 3 , A 4 can be obtained with the potential function. In other words, all the potentials in a tetrahedral element can be obtained by the potentials of four nodes, so the next step is to obtain the potential values of all unknown nodes.
We need to define the characteristics of the units, this paper uses the energy minimum method, so the first calculation is the energy in each unit.
According to the energy of Poynting theorem, in a three-dimensional static electric field, the energy of the electric field is as follows: where ε is the dielectric constant of the electric field acting area, E is the electric field strength of the area, and the relation with the electric potential U is as follows: where S is a differential matrix for the vertices of a tetrahedron. After obtaining the energy of each tetrahedral unit and assembling all the units together, the total energy of the established model can be obtained, denoted as W. In this paper, the description of the total energy W of the model is the total energy of solving the regional electric field, which can be regarded as the sum matrix of each unit matrix-that is to say, the U matrix and S matrix described in Equation (8). U matrix is the potential matrix, which is composed of the node potential of all tetrahedral units, and the S matrix is composed of the matrix of each unit in the model. As shown in Figure 6, the small model consists of three tetrahedral units. For the tetrahedral Unit A, it consists of four nodes: 1, 2, 4, 5; for the tetrahedral Unit B, it includes four nodes: 1, 2, 3, 4; for the tetrahedral Units C, it includes four nodes: 1, 2, 3, 6. Then, based on the fixed-solution condition, the linear equations of the potential of the model nodes can be built. The solution conditions are set herein: The overall energy is minimal. According to the Thomson principle in the electromagnetic field, under boundary and initial conditions, the distribution of the electromagnetic field must meet the minimum energy of the electromagnetic field. Therefore, after obtaining the total energy of the electromagnetic field of the solution region, according to the principle of finding the extreme value, the potential of each node in the solution region can be offset and set to zero. The overall energy is minimal. According to the Thomson principle in the electromagnetic field, under boundary and initial conditions, the distribution of the electromagnetic field must meet the minimum energy of the electromagnetic field. Therefore, after obtaining the total energy of the electromagnetic field of the solution region, according to the principle of finding the extreme value, the potential of each node in the solution region can be offset and set to zero.
When the model is set up, the boundary needs to be processed. Therefore, a certain number of node potential is needed to obtain. The total number of nodes called sum_p, the number of known potential nodes known_p, call it kp, Unknown node number unknown_p, called unkp.
U p denotes the known node potential matrix U f denotes the unknown node potential matrix. Therefore, the total energy is calculated in such a way as to: The U kp T are known, and S c and S b can be obtained by calculation. Only the U kp T are unknown.
According to the S matrix property can be obtained: For easy to calculate, the potential value of the upper interface is set to 10, and the potential value of the lower interface is set to 0. The potential values of all unknown nodes can be obtained by solving the above equations.
The general finite element equation constructed in this paper is a large sparse matrix symmetric algebraic equation group. For solving its solution, a lot of software packages have been used as a solution tool. In this paper, the potential values of all unknown nodes by using the 'taucs' software package can be calculated. According to the model established by the proposed, any potential values in the model can be calculated by obtaining the parameters of each node and the interpolation function of each unit. Several cross-sections are selected, and their potential distribution is shown in Figure 7. When the model is set up, the boundary needs to be processed. Therefore, a certain number of node potential is needed to obtain. The total number of nodes called sum_p, the number of known potential nodes known_p, call it kp, Unknown node number unknown_p, called unkp.
denotes the known node potential matrix denotes the unknown node potential matrix. Therefore, the total energy is calculated in such a way as to: The are known, and and can be obtained by calculation. Only the are unknown. According to the S matrix property can be obtained: For easy to calculate, the potential value of the upper interface is set to 10, and the potential value of the lower interface is set to 0. The potential values of all unknown nodes can be obtained by solving the above equations.
The general finite element equation constructed in this paper is a large sparse matrix symmetric algebraic equation group. For solving its solution, a lot of software packages have been used as a solution tool. In this paper, the potential values of all unknown nodes by using the 'taucs' software package can be calculated. According to the model established by the proposed, any potential values in the model can be calculated by obtaining the parameters of each node and the interpolation function of each unit. Several cross-sections are selected, and their potential distribution is shown in Figure 7. The relation between the macroscopic dielectric constant M and W of a cube is: The U in the equation represents the voltage at both ends of the plate, and L represents the side length of the cube. Therefore, for our three-dimensional model, the expression of macroscopic dielectric constant is as follows: The relation between the macroscopic dielectric constant M and W of a cube is: The U in the equation represents the voltage at both ends of the plate, and L represents the side length of the cube. Therefore, for our three-dimensional model, the expression of macroscopic dielectric constant ε M is as follows: The expected macroscopic dielectric constant epsilon ε M can be obtained by solving the above equations.

Results
According to relevant information, it is known that the permittivity of nickel-zinc ferrite is generally 10-1000, and that of PZT is generally 460-3400 [36][37][38][39][40]. Therefore, the dielectric constant ratios of dual-phase particles are set as 1:500. The results are shown in Figure 8: The expected macroscopic dielectric constant epsilon can be obtained by solving the above equations.

Results
According to relevant information, it is known that the permittivity of nickel-zinc ferrite is generally 10-1000, and that of PZT is generally 460-3400 [36][37][38][39][40]. Therefore, the dielectric constant ratios of dual-phase particles are set as 1:500. The results are shown in Figure 8: According to the above results, it can intuitively see that the larger the volume ratio of Ferrite particle is, the smaller the macroscopic dielectric constant is. Where the black line is Shut, V. N's experimental result [31].
The simulation results reveal that mathematical method can be a new way to solve the preparation dual-phase composites: The volume ratio of Ferrite particles less than 20%, dielectric constant of PZT matrix ferroelectric composites decreases significantly. This method can be played a guide to make composites preparation when it is necessary to customize the PZT matrix ferroelectric composites with dielectric constant requirements.

Conclusions
In this article, a 3D microstructure model of PZT matrix ferroelectric composites is established. And a numerical method is developed to analyze the electrical properties for composites by a threedimensional microscopic model. The dielectric constant and cross-section potential can be calculated by controlling the volume ratio of the two-phase particles in the model. By building a model of PZT matrix ferroelectric composites, might lead experimenters to have a target for preparation without the unnecessary experimental process.
We come to the following conclusions: 1. According to X-ray (micro-) diffraction microscopy provides 3D microstructural phase representations, and the inputs (such as the number of spheres in unit volume, the duty ratio per unit volume, the radius of the generated sphere and distribution pattern) were obtained by stereological techniques. The inputs can be adjusted to simulate the PZT matrix ferroelectric According to the above results, it can intuitively see that the larger the volume ratio of Ferrite particle is, the smaller the macroscopic dielectric constant is. Where the black line is Shut, V. N's experimental result [31].
The simulation results reveal that mathematical method can be a new way to solve the preparation dual-phase composites: The volume ratio of Ferrite particles less than 20%, dielectric constant of PZT matrix ferroelectric composites decreases significantly. This method can be played a guide to make composites preparation when it is necessary to customize the PZT matrix ferroelectric composites with dielectric constant requirements.

Conclusions
In this article, a 3D microstructure model of PZT matrix ferroelectric composites is established. And a numerical method is developed to analyze the electrical properties for composites by a three-dimensional microscopic model. The dielectric constant and cross-section potential can be calculated by controlling the volume ratio of the two-phase particles in the model. By building a model of PZT matrix ferroelectric composites, might lead experimenters to have a target for preparation without the unnecessary experimental process.
We come to the following conclusions: 1.
According to X-ray (micro-) diffraction microscopy provides 3D microstructural phase representations, and the inputs (such as the number of spheres in unit volume, the duty ratio per unit volume, the radius of the generated sphere and distribution pattern) were obtained by stereological techniques. The inputs can be adjusted to simulate the PZT matrix ferroelectric composites model's characteristics in different conditions. 2.
The ferrite/PZT dual-phase composite model generated by the proposed method is similar to the SEM observation results. This article makes a comparison between simulated particle size and actual particle size (i.e., Figures 1 and 2) to increase persuasion. 3. The 3D model generated by the proposed technique, and the comparison of the effective dielectric constant with experimental and numerical results from the literature (i.e., Figure 8) can be accounted as promising indicators for the resemblance of the model with the real materials.

Conflicts of Interest:
The authors declare no conflict of interest.