Microsoft Word - 48dallavecchia.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 49, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Enrico Bardone, Marco Bravi, Tajalli Keshavarz Copyright © 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-40-2; ISSN 2283-9216 Optimization of an ad hoc Realized Space Frame Structured RVE for FEM Modeling of Nanoporous Biopolymeric Scaffolds Obtained by Supercritical Fluids Assisted Process Francesco Naddeo*a, Lucia Baldinoa, Stefano Cardeaa, Alessandro Naddeoa, Ernesto Reverchona,b a Department of Industrial Engineering, University of Salerno, Via Giovanni Paolo II, 132, 84084, Fisciano (SA), Italy b NANO_MATES, Research Centre for Nanomaterials and Nanotechnology, University of Salerno, Via Giovanni Paolo II, 132, 84084, Fisciano (SA), Italy frnaddeo@unisa.it In this work an hypothesis of modeling nanofibers network of Poly-L-lactide (PLLA) scaffolds loaded with hydroxyapatite (HA) nanoparticles, suitable for tissue engineering applications, is presented to investigate the mechanical properties by FEM analysis. Scaffolds were produced by Supercritical CO2 drying of polymeric gels. FEM modeling of nanoporous biomaterials involves computational problems such as: the reproduction of the nano-morphology by means of different techniques, such as molecular dynamics simulations that hardly lead to results adherent to the experimental evidence; the reverse engineering to extrapolate geometric information; the identification of a periodic representative volume element (RVE) to reach more coherent results and to reduce the simulation computational effort. Basic modeling assumptions are: a) polymer particles are small enough to exhibit Brownian motion; b) scaffold SEM images show that the porous structure consists of curved fibers that depart from punctiform nuclei realizing a space frame; c) scaffold experimental compressive tests show that the porous material behaves as a soft isotropic material. On the basis of these assumptions, a parametric algorithm that creates a cubic RVE, showing a nanofibers network and having the same porosity of the real material, has been written; RVE size has been optimized on the bases of its material isotropic degree measured by an ad hoc created iterative algorithm for generating a rod spaceframe; RVE mechanical behavior has been optimized by curving appropriately each fiber according to the experimental data and on the basis of SEM imaging diagnostic. Linear FEM simulations on mechanical behavior have given qualitatively and quantitative satisfactory results when compared to the experimental ones. 1. Introduction Modeling and simulations offer a way to support the development of scaffolds to be used as human tissue substitutes. Computer-Aided Tissue Engineering (CATE) is based on additive manufacturing techniques for the production of patient-specific scaffolds, starting from geometric data obtained from medical imaging. The definition of the internal morphology of the scaffolds is considered one of the major problems within the application of CATE. In the literature, we can find some studies about non parametric computational models based on the finite element method (FEM) and computational fluid dynamics, developed to analyze composite scaffolds for bone regeneration (Milan et al., 2009). In most cases, a micro-tomography scan of the scaffold has been used (Hu et al., 2009). The results showed that scaffold morphology determines cells migration, when they are treated in a bioreactor (Jungreuthmayer et al., 2009). In particular, fluid perfusion induces a stress distribution that is largely dependent on the distribution of pores within the scaffold (Olivares et al., 2009). Other papers used periodic minimal surfaces and forming processes based on rapid prototyping techniques to optimize the scaffold microstructure (Kapfer et al., 2011). However, modeling hypotheses are rarely supported by experimental micro and nanometric information about the scaffold. DOI: 10.3303/CET1649029 Please cite this article as: Naddeo F., Baldino L., Cardea S., Naddeo A., Reverchon E., 2016, Optimization of an ad hoc realized space frame structured rve for fem modeling of nanoporous biopolymeric scaffolds obtained by supercritical fluids assisted process, Chemical Engineering Transactions, 49, 169-174 DOI: 10.3303/CET1649029 169 Aerogels generally show a great structural variety, from branched to compact clusters, depending on either intrinsic or extrinsic conditions (water content, solvent, pH, reaction time, catalysts, etc.). In the literature, it is possible to find various attempts to model the microstructure of aerogels: Diffusion-Limited Aggregation (DLA) (Ma et al., 2001), Void Expansion Methods (VEM) (Schenker et al., 2009) and Gaussian Random Field methods (GRF) (Quintanilla et al., 2003). In the past, the deviation of aerogels from the predicted values for foams was attributed to the abundance of “dead-ends” (Ma et al., 2002), that are clusters connected to the backbone of the aerogel network at only one point. Other authors used the Diffusion Limited Cluster-cluster Aggregation (DLCA) algorithm (Ma et al., 2001) to generate 3D on-lattice aerogel models. However, DLCA model contains excessive dead ends that lead to an underestimation of the gel Young modulus at a given density. At present, there is no consensus on which simulation strategy is the most satisfactory. In our previous work (Baldino et al., 2015), scaffolds showed a nanometric fibers network that was modeled, in FEM environment, as a space frame composed of basic tetrahedral elements whose beams were arranged randomly enough to make it behave as an isotropic material and appropriately curved to fit the experimental results. The choice of the structure based on tetrahedra, introduced in the model geometric constraints that do not exist in reality and could be an obstacle for the future development of the model. Therefore, in this work is presented a further development of the parametric model to approximate more and more the real mechanical behavior of the scaffold specimen nanostructure. Mechanical characteristics of Poly-L-lactide (PLLA) scaffolds loaded with hydroxyapatite (HA) nanoparticles were modeled with the support of Scanning Electron Microscope (SEM) images. A new interpretation of the nano-porous morphology and of the stress transmission mechanisms was given on the basis of a parametric FEM model managed by an ad hoc created algorithm. 2. Materials and Methods PLLA aerogels at 15 % w/w loaded with HA nanoparticles at 10, 20, 30 and 50 % w/w with respect the polymer, were prepared according to the procedure described in a previous work (Reverchon et al., 2009) and processed by supercritical gel drying at 200 bar, 35 °C for 4 hours (Cardea et al., 2013; Baldino et al., 2015; De Marco et al., 2015). In Figure 1, a SEM image of a PLLA+HA scaffold is reported. Figure 1: PLLA+HA scaffold morphology processed by supercritical gel drying at 200 bar, 35 °C. PLLA+HA scaffold is characterized by a regular nanofibrous morphology, with a mean fiber diameter of about 200 nm. Moreover, HA nanoparticles decoration around PLLA nanofibers is also visible, probably due to electrostatic forces among the polymer and the nanoparticles. This kind of morphology is suitable for tissue engineering applications, since the average fiber dimension is comparable with the extracellular matrix of natural tissues. In Table 1, the Young modulus of each sample is indicated. As expected, the Young modulus of the produced scaffolds increases when the HA amount increases, too. These results are consistent with the ones described in a previous work on the same kind of composite scaffolds (Reverchon et al., 2009). Table 1: Young modulus of produced PLLA and PLLA+HA scaffolds by supercritical gel drying HA, % w/w PLLA 0 10 20 30 50 Young modulus, kPa 79 ± 1.4 108 ± 1.7 118 ± 2.0 121 ± 2.2 124 ± 2.8 In particular, Young modulus changed from 79 kPa for PLLA scaffold to 124 kPa for PLLA+ 50% HA scaffold. 170 2.1 Imaging diagnostic A statistical investigation on the morphological parameters of the aerogel was conducted (diagnostic imaging techniques) on the SEM images of the samples. Using an image analysis software, Adobe Photoshop (Adobe Systems Software Ireland Ltd., UK), SEM micrographs of the scaffolds were converted in grayscale digital images. The investigation was made with a semi-automatic procedure which was based on the assumption that the fibers visible in SEM images, having a homogeneous content of gray, were those lying on the images plane (Cricrì et al., 2012); these fibers were chosen for determining the geometric characteristics with a good approximation. The fibers were mostly arranged in random networks in each sample without apparent dead ends; no clusters or aggregates were found in the fiber network. The statistical investigation led to the determination of the average aspect ratio of the fibers = ≈ 0.05; being the average diameter of the fibers and the average length (i.e., the distance between the two extremes). In the hypothesis that the fibers were sufficiently symmetrical, it was also calculated the average value ≈ 0.305 of the distance between the center point of the generic fiber and the segment connecting the two ends of the same ( ), being the real average length of the fiber. 2.2 Geometric modeling The fibers constituting the nanoscale structure born from the interaction of particles and clusters of colloidal particles that are characterized by a Brownian motion (random walk). This phenomenon leads to the development of randomly oriented fibers which impart an average isotropic mechanical behavior (Baldino et al., 2015). In this work, we imagined an isotropic space frame of cylinder-shaped fibers to simulate the scaffold fibers network. The fibers converge, necessarily, also as a matter of material continuity, in convergence nodes presumably random distributed, because of the isotropic structural nature of the studied porous material (observed experimentally) and since, during the gelling process, there are no directional force fields (except the gravitational, negligible in this case) that may suggest otherwise. Hence, it was possible to assume that in each node converges a given average number of fibers with a standard deviation not very high because of the nature of the random gelling process. Thus, on these assumptions, once randomly generated a given number of nodes in a 3D virtual space, we modeled this fibers network imposing that a given number of fibers converged in each node, defining together a sort of sphere of influence; the radius of this sphere was determined by the maximum length of the fibers converging in the node. A different geometrical configuration would not have been appropriate since the sphere is the one that best lends itself to justify a material isotropic behavior. Obviously, each fiber, which branched off from the convergence node, intercepted a second node contained in the above mentioned sphere at a distance less than or equal to the radius of the sphere. The knowledge of the average aspect ratio , dimensionless, allowed us to calculate the average radius of the sphere of influence of the convergence nodes (assuming that the spherical range of influence is maintained almost constant for each nodes) as a function only of a parameter; i.e., of the average length fiber (assuming that the fibers had almost all the same diameter) being worth the following relationship = that become = . According to what previously assumed, i.e., that the investigated fibers were those lying on the plane of the SEM images, we could assume k = . The nanoporosity experimentally calculated ( = 0.75), gave us information on the number of convergence nodes, thus, on the number of fibers to generate inside the sphere of influence. An ad hoc algorithm was created in the Ansys programming environment that randomly produced nodes in a cubic representative volume element (RVE). Subsequently, a routine cyclically generated the straight beam joining the above mentioned node with each node in its sphere of influence. In this model, the number of fibers introduced in the RVE running in the FEM environment, was calculated according to the following Eq(1): = 4(1 − ) (1) being the length of the cubic RVE ( = 1 for the sake of simplicity); , as mentioned previously, being the real average length of the fiber depending on the curvature that was subsequently imparted to the fibers of the FEM model to converge towards the experimental results; = was connected to the computational burden of the FEM simulations. This parameter was calibrated taking into account the material isotropic degree by means of the isotropic criterion (Naddeo et al., 2014), which allowed us to minimize the computational burden of the simulations (represented by the number of fibers generated per unit volume represented by RVE) while maintaining a high isotropic degree. 171 2.3 Introduction of the curved fibers The fibers observed in SEM images (Figure 1) showed an average curvature that, in our opinion, led to a bending/buckling tendency, determining a higher weakness of the entire structure and the high decrease of the Young modulus (Baldino et al., 2015). This observation was supported by previous studies (Ma et al., 2002, Pirard and Pirard, 1997). For this reason, we introduced a routine in the main algorithm that substituted the beams of the fibers network with parametric curved beams based on spline curve (Baldino et al., 2015) interpolating three points for each fiber, two of which coincident with the endpoints of the fiber of length and a third in a central position at a distance = ℎ from the fiber axis, ℎ 0,1 . 2.4 FEM Element type The algorithm used quadratic three-node beam element in 3-D, suitable for analyzing slender to moderately stubby/thick beam structures. Six degrees of freedom occurred at each node for translations and rotations about the x, y and z directions. 2.5 FEM Mechanical properties The Linear elastic material properties were considered for evaluating Young modulus to compare with experimental Young modulus. The interest was focused on the percentage decrease of the characteristics of compact material compared to the characteristics of nanoporous material. 2.6 Boundary condition In a FEM environment, RVE should be large enough to contain all the intrinsic characteristics of the structure to be simulated, which, in our case, consisted of material properties and geometrical characteristics. Cubic-shaped RVE was chosen to allow easy boundary conditioning; intrinsic statistical characteristics were taken into consideration when choosing its dimensions in order to obtain an acceptable isotropic degree (Cricrì et al., 2012). A FEM routine created identical distribution of nodes between opposite edges to apply periodic boundary conditions (Ching et al., 2009; Cricrì et al., 2012; Naddeo et al., 2014). Assuming small strains and elastic behavior of the material, the algorithm, with a single FEM run characterized by six sequential imposed deformations, was able to calculate all the components of the RVE stiffness matrix. An original isotropic criterion was used by Naddeo et al. (2014) to optimize the size of the RVE. The criterion was based on the minimization of the following function: |∆| = , − , ( , ), (2) in which , is the i, j-th component of the stiffness matrix, extrapolated from FEM calculation, , is the i, j-th component of the unknown isotropic stiffness matrix, and λ and G represent the independent parameters defining the isotropic stiffness matrix (Lamè constants). The minimization of Eq(2) led to the determination of λ and G, that were the parameters defining the behavior of the isotropic material closer to the behavior of the material simulated by the Ansys calculation. Consequently, a routine calculated the parameter δ, that was derived from the ratio between the norm of the “difference tensor” ∆= − and the norm of the tensor . The parameter δ provided information about the material isotropic degree of the model, depending on the size of RVE (in terms of α that was introduced in the Eq(1)). 2.7 Convergence analysis of the FEM results A convergence analysis of the results (in terms of Young modulus) was performed allowing the authors to choose an average element size of the FEM model equal to = (L being the length of a generic fiber) with a corresponding percentage error of 0.075. 3. FEM Modeling results The isotropic criterion showed that the value of the parameter δ tended towards a small constant value for = 0.01; in this way, the RVE size was not too large from the computational point of view, but sufficiently large to make the fibers substantially oriented in a random way imparting to RVE a sufficiently isotropic mechanical behavior. Carrying out a sensitivity analysis by varying the curvature of the beams, an output value in terms of Young modulus equal to that of the experimental result for d ≈ 0.295L was obtained. This value of d was consistent with the imaging diagnostic result equal to d ≈ 0.305L . 4. Comparison of experimental results with FEM results The new algorithm managing the nanostructure FEM model, was implemented in the simulation system reported in our previous work (Baldino et al., 2015) in which we modeled scaffolds produced by supercritical 172 CO2 drying of PLLA gels whose structures showed also a microporous architecture generated by the voids left in the solid material by porogen leaching, modeled in the FEM environment as a close-packing of equal spheres based on the theory of hcp (hexagonal close-packed). These scaffolds were also loaded with HA nanoparticles, from 10 to 50% w/w with respect to the polymer, whose effect on the mechanical properties, was modeled taking into consideration the formation of concentric cylinders of HA nanoparticles around PLLA nanofibers; a second cylinder was modeled as partially covering the inner one, introducing a “loaded surface index” defined as = , being the length of the portion of PLLA cylinder fiber covered by HA and the length of the cylinder fiber, as shown in Figure 2. Figure 2: Explanatory example of a cubic RVE of the nanostructured model: PLLA (fibers); HA (covering elements). The curves in Figure 3 show the comparison of experimental and numerical results in terms of Young modulus. Figure 3: Superimposing of the numerical results on the experimental results. FEM results were obtained by combining both the FEM model of the composite micrometric porous structure and the FEM model of the nanoscale structure; the calibration of parameters of the two FEM models was used only to approximate the extremities of the experimental curve. The other values for the modeling were obtained by varying only the percentage of HA in the FEM algorithm. This comparison confirms the plateau effect for HA loadings larger than 30 % w/w, that was experimentally recorded. In this work the “load surface index” was set at a constant value equal to = 0.225, value comparable to that obtained in previous work = 0.14 confirming the hypothesis that the contact on nanofibers occurs only on small areas between PLLA fiber and quasi spherical HA particles, reducing the amount of reinforcement compared to the case when contact between the two materials is continuous (Baldino et al., 2015). 173 5. Discussion, conclusions and perspectives The experimental results were very well approximated by the FEM results confirming that the curvature of the fibers determines the high softness of such materials. Moreover, the proposed new model of nano-porous structure was released from the use of the tetrahedral based space frame optimized in the previous work (Baldino et al., 2015) by taking into account the phenomena that generally occur during the gelling process from a statistical point of view, the experimental isotropic behavior of the scaffolds and the averages geometric characteristics of the nanofibers network extrapolated by means of imaging diagnostic techniques. Furthermore, the model confirmed the hypothesis that the contact on nanofibers generally occurs only on small areas between PLLA fiber and quasi spherical HA particles, reducing the amount of reinforcement compared to the case when contact between the two materials is continuous (Baldino et al., 2015). Currently, we are working on the introduction of the non-linear behavior of the materials involved and on the possibility of imparting directional characteristics simply by hot-stretching these materials in order to drive the largest amount of fibers along an appropriate direction. Parametric FEM modeling is aimed to find the sensitive parameters to control during the production phase of the scaffolds, in order to optimize scaffold performance. References Baldino L., Naddeo F., Cardea S., Naddeo A., Reverchon E., 2015, FEM modeling of the reinforcement mechanism of Hydroxyapatite in PLLA scaffolds produced by supercritical drying for Tissue Engineering applications, J. Mech. Behav. Biomed. Mater. 51, 225 – 236. Baldino L., Cardea S., Reverchon E., 2015, Natural aerogels production by supercritical gel drying, Chem. Eng. Trans. 43, 739-744. Cardea S., Baldino L., De Marco I., Pisanti P., Reverchon E., 2013, Supercritical gel drying of polymeric hydrogels for tissue engineering applications, Chem. Eng. Trans. 32, 1123-1128. Ching W., Rulis P., Misra A., 2009, Ab initio elastic properties and tensile strength of crystalline hydroxyapatite, Acta. Biomater. 5, 3067–3075. Cricrı` G., Garofalo E., Naddeo F., Incarnato L., 2012, Stiffness constants prediction of nanocomposites using a periodic 3D-FEM model, J. Polym. Sci. Part B Polym. Phys. 50, 207–220. De Marco I., Baldino L., Cardea S., Reverchon E., 2015, Supercritical gel drying for the production of starch aerogels for delivery systems, Chem. Eng. Trans. 43, 307-312. Hu Z., Notarberardino B., Baker M., Tabor G., Hao L., Turner I., Yang L., 2009, On modeling bio-scaffolds: structural and fluid transport characterization based on 3-D imaging data, Tsinghua Sci. Technol. 14, 20– 23. Jungreuthmayer C., Jaasma M.J., Al-Munajjed A.A., Zanghellini J., Kelly D.J., O’Brien F.J., 2009, Deformation simulation of cells seeded on a collagen-GAG scaffold in a flow perfusion bioreactor using a sequential 3D CFD-elastostatics model, Med. Eng. Phys. 31, 420–427. Kapfer S.C., Hyde S.T., Mecke K., Arns C.H., Schroder-Turk G.E., 2011, Minimal surface scaffold designs for tissue engineering, Biomaterials 32, 6875–6882. Ma H.-S., Prévost J.-H., Scherer G.W., 2002, Elasticity of DLCA model gels with loops, Int. J. Solids Struct. 39, 4605–4614. Ma H.-S., Prévost J.-H., Jullien R., Scherer G.W., 2001, Computer simulation of mechanical structure– property relationship of aerogels, J. Non-Cryst. Solids 285, 216–221. Milan J.-L., Planell J.A., Lacroix D., 2009, Computational modelling of the mechanical environment of osteogenesis within a polylactic acid–calcium phosphate glass scaffold, Biomaterials 30, 4219–4226. Naddeo F., Cappetti N., Naddeo A., 2014, Automatic versatile parametric procedure for a complete FEM structural analysis of composites having cylinder-shaped reinforcing fibres, Comput. Mater. Sci. 81, 239– 245. Olivares A.L., Marsal E`., Planell J.A., Lacroix D., 2009, Finite element study of scaffold architecture design and culture conditions for tissue engineering, Biomaterials 30, 6142–6149. Pirard R., Pirard J.-P., 1997, Aerogel compression theoretical analysis, J. Non-Cryst. Solids 212, 262–267. Quintanilla J., Reidy R., Gorman B., Mueller D., 2003, Gaussian random field models of aerogels, J. Appl. Phys. 93, 4584–4589. Reverchon E., Pisanti P., Cardea S., 2009, Nanostructured PLLA_hydroxyapatite scaffolds produced by a supercritical assisted technique, Ind. Eng. Chem. Res. 48, 5310–5316. Schenker I., Filser F.T., Herrmann H.J., Gauckler L.J., 2009, Generation of porous particle structures using the void expansion method, Granul. Matter 11, 201–208. 174