Multi-scale Modeling of Polymeric Composites Including Nanoporous Fillers of Milled Anodic Alumina

A polymer composite based on an innovative filler consisting of microscale powder of nanoporous alumina is modeled. The passing-through nanoscale pores in this system—roughly columnar cylindrical, with diameter of the order of 100 nm—are fully penetrated by the resin, which is not bonded to the inner pore walls by any chemical agent. This system, previously assessed by laboratory experiments, is modeled here for the first time, based on a computational multi-scale hierarchical approach. First, microscale representative volume element (RVE) is modeled in two steps using finite element modeling. Then, the macro-scale RVE is characterized, using a combination of micromechanical rules. The elastic response of the composite is simulated to predict its Young’s modulus. This simulation confirms the former experimental results and helps to shed light on the response of the investigated material, which may represent a novel system for use in disparate composite applications. In particular, the nanoporous microfillers composite is compared with a composite material containing the fillers of the same material yet nonporous, bonded to the matrix. It appears that, with respect to this standard concept of three-phase composites, the presence of the nanopores can compensate for the absence of the bonding agent.


Introduction
The field of polymer-based composites is ever growing in engineering applications and-as a consequence-in research striving at improving the respective material performances [1][2][3][4][5]. More and more, in a number of applicative areas such as automotive and aerospace, consumer products and marine industry, the demand grows for replacement of traditional metal-based materials by composites, which should provide the same or even better performance with additional advantages: lightweight (turning into lower energy consumption for transportation and decreased greenhouse gas emission), high recyclability, low cost and-in the biomedical area-increased biocompatibility [6]. One means to obtain performant composites with the above properties has often been devised to be the use of porous reinforcing fillers, which has especially been investigated in the field of dental restorative materials [7][8][9][10][11][12]. In the biomedical area, nanoporous-filler-based composites may also become in per-spective smart, at least for the functionality of drug release, for, e.g., antibacterial tasks [13][14][15].
Actually, restorative dentistry is one of the major areas of research in polymeric composites, where the used polymer matrices are resins. On average, every 5 years all major dental companies introduce to the market one new dental restorative composites per application line that should present higher performances as compared to the previous generation. This has been happening at least during the past three decades, yet it is difficult to detect really innovative concepts and impressive improvements in the performance of any such materials [16]. The most remarkable advancements seem to have occurred in the much claimed aesthetic quality, which has also been the real reason for banning amalgam. However, this aesthetic quality is a questionable property [17], when also considering that its common description as of natural appearance is often erroneously called biomimesis, while the structure and function of artificial restorative materials are quite different from those of the native tooth materials.
Advancements in dental composites may in principle occur in either of two directions, namely the resin or the filler formulation [16,[18][19][20]. However, as a matter of fact the second one has been dominating in the recent history, even giving the name to the different classes of composites. Particularly, for truly potentially innovating composites, more exotic types of fillers have been introduced recently [21], of which class our proposal of nanoporous fillers may also be considered an example.
In the previous work, a dental restorative resin composite based on nanoporous microfillers of aluminum oxide, obtained from anodic porous alumina (APA), was proposed and demonstrated, exhibiting improved response to aging in wet ambient for the elastic modulus as measured by dynamic mechanical analysis (DMA) [22]. The novel experimental composite is based on a standard bisGMA-based matrix used for many years in dental restorative composites [19,23] and an innovative nanoporous microfiller obtained from anodic porous alumina membranes [24,25] separated from the metallic substrate and ball-milled. The composite is only dual component (resin matrix-inorganic filler), without the presence of bonding agent at their interface [22]. It was shown that the mechanical interlocking of the resin through the pores of the microfillers replaces the bonding agent, and the experimental composite exhibits better aging stability in water for the elastic modulus than standard composites, likely due to the lack of chemical bonding agent degradation after hydrolysis occurring in the wet oral environment. In particular, the diminished degradation in elastic modulus on aging was observed at body temperature and at 1 Hz strain frequency. The strength of adhesion of this innovative composite to dental tissue was also characterized in a subsequent experiment [26], resulting in only little lower performance than the average of current ternary restorative composites.
The nanoporous nature of the described composite is also promising for drug loading function, eventually enabling antibacterial or remineralizing capabilities of the material [15,27,28].
In this work, in order to support former work and justify the prosecution of investigations in our concept of a nanoporous microfiller dental resin composite, the mechanical response of the investigated material is simulated through a computational multi-scale modeling procedure. The goal of the work is twofold. On the one hand, to confirm and analyze the experimental results obtained so far with our experimental nanoporous microfiller composites. On the other hand, to provide a simple mean of comparison with traditional composites, based on the same components as to both matrix and filler bulk material properties, yet without porous nanostructuring of the fillers, both with and without the presence of a bonding agent.

Modeling Procedure
Hierarchical multi-scale modeling is developed for the purpose of this study, where the outputs of a scale are fed into the next scale as input data. First, the mechanical properties of isolated APA are obtained using finite element (FE) modeling and the influence of APA size on the results is studied. Then, the FE model of a microscale RVE containing APA and surrounding polymer is constructed. At this scale, the interaction between APA and polymer is taken into account, while the influence of APA size on the interaction is also investigated. Finally, the random orientation of APA fillers inside resin is modeled at the final scale of macro, where laminate analogy technique is combined with generalized methods of cell for obtaining the Young's modulus of the investigated material; the later accounts for the distribution of APA fillers in matrix and the former captures their random orientations.

Properties of Composite Materials and Structure
The parameters of the individual composite components used for the modeling, including both the material properties and the geometrical parameters describing the composite structure, are summarized in Table 1. For the BisGMA-based matrix, we assumed: modulus E m 3 GPa, density ρ m 1.1 g/cm 3 , where "m" stands for matrix. For the bulk amorphous hydrated alumina, as supposed to be resulting from anodization, we assumed: modulus E ba 60 GPa, density ρ ba 2.8 g/cm 3 , where "ba" stands for bulk alumina. For both materials, we assumed a Poisson ratio ν 0.3. Models of the different-sized APA fillers are shown in Fig. 1a.

FE Modeling of the Isolated APA Fillers
The APA filler is constrained at one end and an axial displacement is applied on the other end in order to find the longitudinal and transverse moduli. The reaction forces are extracted from the output of FE modeling, and the moduli are calculated using the following equation: where F R is the summation of reaction forces on the constrained surface and l represents the distance between constrained and deflected surfaces. l denotes the applied displacement and A is the area of the deflected surface. The same boundary conditions are employed on each axis in order to find the APA's different properties. A presentation of the model and applied boundary conditions are presented in Fig. 2 in both directions. FE modeling is conducted using commercial Abaqus CAE package with the standard/explicit model [29].
Decreasing the size of elements, different models with various mesh densities are constructed and analyzed. The size of the elements is chosen small enough to ensure independency of the results from mesh density.

Modeling Microscale RVE
In the second step, microscale RVE consisting of pristine APA and surrounding matrix is modeled in a cubic cell (see Fig. 3). The penetration of matrix into the pores is also taken into account and modeled. It is assumed a shear modulus G 2.50 GPa, and different Poisson's ratios in the different directions (ν x 0.35, ν y ν z 0.26).
When modeling the interaction at the interface between APA fillers and matrix, penalty friction properties using a tangential behavior are applied to the contact surfaces. This option is widely used by researchers in similar contexts, to capture the behavior of this kind of contact surfaces [30][31][32][33]. The friction coefficient is selected between 0.5 and 100 (highest possible amount), as compared to a frictionless interaction, for which it would be set to zero. One end of the model is restricted from any movement and axial displacement is applied to the other end. The effect of this friction coefficient parameter has been evaluated for the case of a 50%wt loading fraction of the fillers (see Sect. 3).
To investigate the effect of different filler sizes and the respective interaction with the matrix on the results, several FE models are constructed. One of these is shown in Fig. 4. The smallest model with an APA filler consisting of seven pores only is chosen among all, for the sake of visibility. Tetrahedral elements were used in a mapped meshing for matrix, and a free mesh was used for the filler. Two types of surface interactions are visible in the representative colored surfaces in the model of Fig. 4: interpore interactions caused by the penetration of the matrix and external interactions at   Fig. 4 is much lower than the implemented model, and Fig. 4 is just a schematic presentation of meshing system and friction surfaces. If the actual model was shown, the meshing system and friction surfaces would not be visible at all because of insignificant ratio of pore size to filler size.
For the case of longitudinal and transverse modulus, the same boundary conditions described in Fig. 2 are applied to the RVE instead of isolated APA and Eq. (1) is used. For the Poisson's ratio, the following equation holds: where ε l and ε t represent longitudinal and transverse strains, respectively.
To obtain the shear modulus, a shear deflection is applied to the free end while the other end is fixed and then following equation is used: where F R represents the total reaction forces on the constrained surface, L is the distance between constrained and deflected surfaces, x is the shear deflection magnitude and A is the area of the deflected surface.

Modeling Composites at Macro-scale
In this step, the macro-scale RVE of the composite is modeled using a combination of micromechanical techniques to estimate the macroscopic Young's modulus of investigated materials. The effects of filler size distribution and different frictional interactions between constituents have already been modeled. Another parameter affecting the macroscopic mechanical properties of the investigated composite is the random orientation of the fillers in the macro-scale RVE. Laminate analogy (LA) is the technique of substituting the randomly oriented composites with a quasi-isotropic laminate, as first developed by Halpin and Pagano [34]. In LA, the randomly oriented short-fiber composite can be replaced with a four-layer quasi-isotropic laminate and the properties are interchangeable. This technique has been proven by several studies on the stiffness and strength properties [35,36].
The generalized method of cells (GMC) has also been developed first by Paley and Aboudi [37]. In GMC, the composite RVE is divided into sub-cells in which the continuity of deflections and tractions are applicable between cells. Among other available averaging methods, the GMC has the advantage of considering mechanical interactions and stress/strain transfer between the sub-cells. The formulation of GMC has been provided briefly in "Appendix A".
LA, GMC and classical lamination theory (CLT) were used in this step to obtain the overall properties of the investigated composite. The combination of these methods has been successfully used previously by Rafiee and Eskandariyun for predicting mechanical properties of graphene/polymer nanocomposites [38]. In this study, the developed strategy was customized for the purpose of the current research. A quasi-isotropic laminate was considered with lay-up configuration of [0/90/ ± 45]. Each of these angle plies included 12 The modeling steps can be briefly outlined as three stages: (i) The modeling starts with extracting the effective modulus of isolate APA filler and analyzing the influence of geometrical parameters; (ii) the interaction between APA and resin is modeled through FE analysis of microscale RVE; (iii) LA is combined with GMC for capturing random orientation and also distribution of APA fillers.

FE Modeling of the Isolated APA filler
As a first step in the model, the pristine APA nanoporous microfiller is designed and analyzed using FE, in order to obtain its equivalent elastic modulus, as compared to the bulk alumina fillers with similar external size. The properties of the starting filler and resin materials are summarized in Table  1 (Experimental section). It was assumed that the size distribution of the APA fillers be Gaussian, with mean diameter of 5 µm and standard deviation of 3 µm. Since such a distribution would show a tail on the short size hand side that would extend down to negative values, which is unphysical, this tail has been cutoff, resulting in the distribution profile shown in Fig. 1b. In Fig. 1a some fillers within this distribution range are shown, which allows one to see the approximate number of pores present inside the typical single filler object.
After single filler modeling, analyses of the same have been done using FE. According to the directions defined in   Table 2 and plot in Fig. 6.
It can be seen that there is no major difference in Young's modulus for the different APA sizes. From Fig. 6, it appears that both moduli undergo the largest change in value for the smallest considered APA filler size of 1 µm diameter. Obviously, on this scale only few pores (just two or three) are span in each transversal direction (meaning that only between four and nine are included in the whole filler), and the boundary effects become quite relevant. Nevertheless, as shown by the average modulus profile, the deviations in the different directions tend to compensate. Anyway, for fillers larger than~5 µm, both the longitudinal and transverse moduli roughly converge, to~43.0 GPa and~35.0 GPa, respectively.

Microscale RVE
The effect of friction coefficient between APA filler walls and fully pore penetrating resin, as variable in the chose range   Fig. 7. Obviously, while for low friction coefficients (approximately below 5), there is a clear deviation toward lower values of the modulus, above 5 the modulus tends asymptotically to a stable value of~9.3 GPa. In any case, even the lowest value (~8.62 GPa) does not present a major deviation from the asymptotic value (~− 7%).
Since no major difference is seen on the final results, thus the results of the model with mild friction coefficient value of 0.5 were considered next and are presented in Table 3. Due to the transversely isotropic behavior of the RVE, five engineering constants are reported, which are defined as follows: E 1 is the longitudinal elastic modulus, whereas E 2 E 3 are the transversal ones; similarly, ν 1 and ν 2 ν 3 are the longitudinal and transversal Poisson ratios, respectively; G is the shear modulus.
It should be noticed that the observed differences for the Young's modulus in the two independent directions are roughly in line with those appearing formerly for the single filler case.
These results were obtained for the case of aligned APA fillers in the composite. Hence, for a realistic model of the composite, where the APA fillers are oriented randomly at all directions in space, the longitudinal modulus mentioned above will be an upper-bound limit.

Macro-scale Composite
Using the LA-GMC approach implemented as described in the experimental section, we obtained the macro-scale modulus of the composite with randomly oriented APA fillers. The results for different weight loading fractions are reported in Table 4 and plotted in Fig. 8. Figure 8 shows that, for a loading fraction of 50%wt, the modulus of the APA composite should fall between~5 and~9 GPa. Actually, previous measurements of the same experimental APA-based composite by means of DMA [22] reported values of storage modulus around 5.2 GPa, so within the range found here, yet closer to the lower limit.
In Fig. 8, the macro-scale results represent the lower bound, and they are compared with the microscale RVE results, which represent the upper bound (case of aligned APA fillers). Along with these two results, the rule of mixture (ROM) estimation is also presented, which has been calculated by using the pristine APA filler modulus. The modulus of the investigated composites with any APA weight fraction is bounded between the lower and upper limits.
For the case of 50%wt loading only, the modeling was repeated for the control material with which the present concept of composite has to be compared with, namely the composite made with same materials but nonporous and yet bonded filler. The resulting modulus is found to be 5.13 GPa, which corresponds to the red circular spot in Fig. 8. While higher than the 50%wt result for the APA-based composite object of the present work, which is found to be 4.89 GPa, this value is not significantly different. Obviously, the APA filler porosity can compensate for most of the action of the bonding agent present in the traditional nonporous and bonded alumina composite.

Discussion
The values of the mechanical properties of the single composite components are probably affected by a rather large uncertainty. This holds especially true for the elastic modulus, and particularly for the filler bulk material. In fact, in the literature, Xia et al. [39] reported an elastic modulus of bulk alumina as high as 140 GPa, whereas Fang et al. [40] reported the quite low value of 5 GPa, when including the porosity as in APA. Ko et al. [41], in turn, reported variable values between 70 and 140 GPa. Here, considering that the alumina resulting from anodization is both mainly amorphous and is additionally hydrated (with water contents around 5%wt [need ref]), we assumed E ba 60 GPa. However, the uncertainty could be quite high, probably up to ± 50%. For the resin modulus, the uncertainty is probably much lower, around ± 15%. In some respect, BisGMA can be thought to be similar to an epoxy, with a typical modulus of~2 GPa. However, the actual formulation of the matrix containing BisGMA usually includes some additional comonomers, to make the paste less thick and viscous, such that it can be actually mixed with the filler powder. In fact, Emami et al. [42] reported, for pure Bis-GMA, a modulus of 4.15 GPa, whereas for pure TEGDMA, a value of 2.37 GPa. We also used TEGDMA as a comonomer, in the wt ratio of 70:30 BisGMA-TEGDMA, which, according to the rules of mixtures, could correspond to an upper limit (weighted mean) of 3.6 GPa. To be conservative, we assumed thus a matrix modulus of 3 GPa. It should also be considered that this additionally depends on the actual degree of conversion reached on photopolymerization [43,44].
It should be reminded here that the APA filler of the experimental composite is not bonded to the matrix, and the only "bonding" present at the interface between the two phases is a friction coefficient assigned to them for the relative displacement of one with respect to another. This has been selected as described in the experimental section. Whereas the upper limit of friction coefficient range considered in the study presented in Fig. 7 may seem unrealistically high, nevertheless we observed that the deviations found in Young's modulus for the highest friction coefficients were only limited (+ 7.5%), and thus we finally assumed a moderate value of 0.5.
In the step from microscale RVE to macro-scale modeling, the random orientation of the APA fillers is introduced. As one possible limitation of our model, it is worth mentioning that the agglomeration of APA fillers is not considered, since APA fillers were very well dispersed in the production process in this study. Therefore, fully dispersed situation is modeled.
It has been stated that the experimental value of 5.2 GPa found previously for the storage modulus of the APA composite for 50%wt loading could probably be increased by further improvement in the mixing process, or by implementing a mixed strategy, also including a bonding agent in support of the resin-filled pores in the filler. This could hopefully result into a modulus of up to 8-10 GPa [22]. The improved filler-matrix mixing, in particular, could reduce the formations of filler agglomerates and thus improve the bonding. This issue affects the modulus similar to the alignment of fillers. Hence, the presented possible area (green area) in Fig. 8 could be reliable even using different experimental parameters. It should be pointed out that currently developed hierarchical modeling very well serves the purpose of this study, since it was intended to predict the modulus. Stochastic or concurrent model might be required to be utilized for more complicated problems like characterizing the fracture behavior [45][46][47].

Conclusion
Previous experiments showed that the modulus of the experimental APA-based composite is not affected significantly by the aging. Therefore, even if still a bit lower than that of traditional (three phases, bonded-fillers) composites, this composite shows the advantage of being more stable over the time. In this work we modeled the material, through the following key stages: (i) modeling isolated APA and studying the influence of geometrical parameters on the results; (ii) modeling the interaction between APA and surrounding polymer through FE analysis as microscale modeling; (iii) combining laminate analogy of generalized method of cell; this is the most important contribution of the current work, which has not been done before. The modeling results presented here supported the statement about an equivalent modulus of the considered composite as that of traditional composites, yet obtained by means of the mechanical interlocking of the resin in the filler nanopores. This is a promising result, in view of the possible advantages of avoiding the bonding agent phase, which may be labile on aging in harsh conditions such as high temperature or penetrating liquids, and of having the inherent capability of the filler pores as carrier for release of active agents into the bulk material, particularly for biomedical applications. In perspective, it will be worth investigating the strain mechanism of the porous filler frame, to check if its elasticity can allow for partial relaxation of stress during deformation, and to model the volumetric shrinkage of the resin, to assess if this could also take advantage of the porous volume reservoir, in case of, apparently defective, incomplete pore filling.
Authors' contributions RR and MS contributed to conceptualization and Supervision; RR and AE contributed to methodology; AE contributed to formal analysis and investigation; MS contributed to writing-original draft preparation; RR, AE and CL contributed to writing-review and editing; RR contributed to funding acquisition and resources.
Funding Open access funding provided by Istituto Italiano di Tecnologia within the CRUI-CARE Agreement. The research was supported by the institutional fund of the COMRESLAB.

Availability of data and material
The data can be made available, on request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: GMC
The GMC is a micromechanical technique developed by Paley and Aboudi to model a multiphase composite consisting of repetitive unit-cells. The GMC presents a generalization of the method of cells, which was previously designed only to model a four-cell network, to an arbitrary number of cells. This micromechanical technique, uses the continuity of displacement and traction rates at the interfaces. The developed homogenization, will results in an effective imitation of the composite's behavior [36].
From the continuity of displacements, we have: Equations A1 to A3 can be replaced in matrix form as: ε (ε 11 · ε 22 · 2ε 12 ) (A5) ε s ε 11 · ε 12 · . . . · ε N β N γ (A6) Similarly, for the continuity of tractions we have: It can be obtained from equation A4 and A7: ε s can be extracted from Eq. A8 as: The average strains to local stresses can be related using stiffness matrices as: The final stiffness matrix will be obtained as [36]: