The Elastic Behaviour of Sintered Metallic Fibre Networks: A Finite Element Study by Beam Theory

Background The finite element method has complimented research in the field of network mechanics in the past years in numerous studies about various materials. Numerical predictions and the planning efficiency of experimental procedures are two of the motivational aspects for these numerical studies. The widespread availability of high performance computing facilities has been the enabler for the simulation of sufficiently large systems. Objectives and Motivation In the present study, finite element models were built for sintered, metallic fibre networks and validated by previously published experimental stiffness measurements. The validated models were the basis for predictions about so far unknown properties. Materials and Methods The finite element models were built by transferring previously published skeletons of fibre networks into finite element models. Beam theory was applied as simplification method. Results and Conclusions The obtained material stiffness isn’t a constant but rather a function of variables such as sample size and boundary conditions. Beam theory offers an efficient finite element method for the simulated fibre networks. The experimental results can be approximated by the simulated systems. Two worthwhile aspects for future work will be the influence of size and shape and the mechanical interaction with matrix materials.

The mechanics of sintered, metallic fibre networks as used for the present study have been investigated by [1,[30][31][32][33][34][35]. In [36], an architectural characterization of the six network samples of the present study was published, together with experimentally obtained mechanical properties.
For finite element (FE) analyses of structures whose dimensions are dominated by their extension along only one axis, such as screws or fibres, beam theory can offer an efficient simulation method [37]. The foundations of beam theory were laid in [38][39][40] and are available in today's textbooks [41]. For the application of boundary conditions (BC) to random fibre networks and for the determination of the representative volume element (RVE), [29,42] have provided a much referred to concept [43][44][45][46]. By this concept, the RVE is obtained for a defined relative error for each physical property. Image acquisition by computed tomography (CT) scanning which uses the principle of Röntgen radiation [47] is a commonly used method for measuring the dimensions of metallic three-dimensional (3D) structures [48].

Motivation and scope of the present study
The motivation of the present study was to predict previously unknown mechanical properties for metallic, sintered fibre networks. For that purpose, FE models were developed. The input geometries for the FE models were based on CT scans, acquired from real network samples in [36]. Experimental values for the network Young's modulus from that same study were used for a validation of the FE models. The properties of these six network samples define the scope and at the same time also the limitations of the present study. Other metallic, sintered networks with similar but different properties will require a re-evaluation of the proposed FE models.

Mathematical notation
Throughout the present study, the following notation is used: x for scalars, x for vectors, and x for 2 nd rank tensors. Cube faces are written X. The vector product is given as "×" and dot product as "Á". Relations which are greater-than and approximately equal are written "≳". If the relation is greater-than or equal "!" is used.

Network samples and meshing step
The present study uses as geometry input for the FE models the dimensions of six AISI (American Iron and Steel Institute) 316L network samples. The architectural network values have been published in [36] (see Table 1). Each scanned sample section is a cube of volume V = 4 3 mm 3 (see Fig 1a). Following the two-phase model in [29], V is split into fibre volume and void volume: Equally, the volume boundary @V is defined to consist of a section of fibre boundary and a section of void boundary on the cube surface S: Six quadratic cube faces form S (see Fig 1b), two of them perpendicular each to one of the three axes: Two samples each are manufactured with a fibre volume fraction f of 10, 15, and 20%. In [36], the fibre production process by bundle drawing [49], the network manufacturing at N.V. Bekaert S.A. (Belgium), the CT scan acquisition at General Electric (resolution 7.75 μm), and the applied skeletonisation algorithm [50][51][52] are discussed in detail. The skeletonisation algorithm reduces the 3D fibre bodies to their medial axes which are strings of voxels (see Fig 2). The medial axis models obtained for [36] are transferred in the present study into beam assemblies and run as FE models. The approximately hexagonal fibre cross-section is simulated as a round cross-section of radius R = 20 μm.

FE models
For the present study, two FE models are implemented (see Table 2). Model A and Model B are based on five specifications: 1. Linear elasticity and static equilibrium, 2. Euler-Bernoulli or Timoshenko beam elements as AISI 316L fibres, 3. Rigid joints or torsional springs as inter-fibre joint models at contact points between fibre paths (see Fig 2), 4. BC as defined in Eqs (8) and (11), and 5. the new introduced model parameter of BC-depth h BC . Linear elasticity and static equilibrium conditions. For mechanical FE analyses, the simulated system is transferred into a representation by the global stiffness matrix K; the most common assembly method being the direct stiffness method [53]. K links the global load vector F and the global displacement vector u: The meshed geometries and material specifications define the variables of K. In the present study, the material is assumed to behave linearly elastic (i.e. fulfil Hook's law that stress and strain are linked by the Young's modulus: σ = E [54]). In this case, the variables of K become constants. The AISI 316L material stiffness is simulated as E fibre = 200GPa with a Poisson's ratio ν = 0.3 which expresses the material's lateral contraction as fraction of the axial extension. For plastic or non-elastic simulations, the variables of K are functions of e.g. force or displacement.
BC are imposed in FE by prescribed values for entries of F and u. The non-prescribed entries of F and u are determined by the FE solver (present study: Abaqus 6.13 [55]) with solutions respecting the equilibrium conditions of forces in Eq (5) and of moments in Eq (6) [56]. Both equations adopt the Lagrangian reference frame [57]. (A comparison to the Eulerian reference frame and its advantages for the modelling of fluids is available in [58] The surface traction vector t is obtained as the product of the Cauchy stress tensor s and the unit outward normal n out . (A detailed discussion of s is available in today's textbooks [59].) In the mechanical simulations of the present study, body force per volume f is neglected; gravitational forces or magnetic forces being typical examples for f . The point vector x specifies the location of a point relative to the origin; i.e. also the equilibrium of moments in Eq (6) being taken about the origin.
Beam elements. One Euler-Bernoulli beam element (B33) and two Timoshenko beam elements, one with linear interpolation (B31) and one with quadratic interpolation (B32), are implemented for the present study. Table 3 contains the complete list of implemented Abaqus elements. The neglected shear strain of the Euler-Bernoulli beam and the advanced Timoshenko beam are further discussed in [41] for the simplified 2D case. Literature recommends in general for beam assemblies which include short beam elements the Timoshenko beam [28,[60][61][62]. For this case, [62] has documented an overestimation of the structural stiffness by Euler-Bernoulli beam elements. Inter-fibre joint models. Model A is simulated for rigid inter-fibre joints. Torsional spring elements (CONN3D2) are inserted into the medial axis models in Model B between intersecting fibre paths (see Fig 2). This concept for adaptive joint strength of Model B is proposed for comparable cases in publications such as [63]. In the present study, the inter-fibre joint stiffness K joint is defined: The scaling factor s [m] in Eq (7) allows the directed variation of K joint for an approximation of the experimental values of [36]. The cross-sectional area A fibre and E fibre are included as they relate the modelled value of K joint to the fibre dimensions and to the fibre stiffness. Alternative model variables could be chosen too.
Kinematic uniform BC and load cases. In [42], a simulation set of kinematic uniform BC (KUBC), static uniform BC (SUBC), and periodic BC (PBC) was proposed for a random heterogeneous composite. KUBC are applied in the present study as defined in Eq (8). In the case of KUBC, a macroscopic strain tensor E imposes the displacement vector u on all x located on @V: Six independent load cases are implemented in the present study by KUBC for the network samples: i = 1 to 3 for the simulation of tensile tests along the axes x, y, z and i = 4 to 6 for the corresponding shear tests. The combination of all six load cases leads to the symmetrical stiffness matrix C, where the six entries C ii on the main diagonal stand for the Young's moduli E and Shear moduli G [64,65]: with : The general 4 th rank stiffness tensor contains in total 81 components. The simplification to only 21 independent variables in C is achieved by energy considerations and symmetry [66].
Mixed BC. The application of further BC to random fibre networks requires additional considerations [29]. In the case of SUBC, a macroscopic traction vector G is imposed on @V. PBC add a periodic fluctuation term v. In [29], it is discussed for random fibre networks that considering: G can't be prescribed on @V void , only on @V fibre . Due to the randomness of the fibre network geometry, the equilibrium conditions in Eqs (5) and (6) are not fulfilled when G is imposed. For overcoming this, mixed BC (MBC) are proposed by [29] and adopted in the present study for the simulation of tensile testing along the axes x and y by one independent load case each (load case i = 1 and 2, definition see above). A macroscopic strain tensor E prim imposes the primary displacement u prim along the axis of tension only: Along the two non-prescribed axes, a secondary displacement u sec is caused as a consequence of the imposed u prim . The total deformation of the sample is obtained under MBC as: The results obtained in [29] for the elastic modulus of computer generated random fibre networks without preferred orientation direction predict greater stiffness for KUBC: BC-depth. The present study extends the BC model of Eqs (8) and (11) for the investigated material by the variable of h BC . The value of h BC prescribes to which depth into the material along the inward normal n in BC are imposed on @V. The reason for the model modification is the obtained structural response which will be discussed further in the following section by Model A in Figs 3 and 4.

Representative volume element determination
For the determination of the RVE size V RVE , the present study relied on an algorithm proposed by [42]. This algorithm has been applied since its publication to random fibre networks in [29,32] and other materials [43][44][45][46].
In [42], V RVE is obtained as a probability estimate. The RVE is understood as the minimum cube volume V RVE for which a particular physical property Z after n realisations can be determined within the margin of the relative error rel . Due to this definition, V RVE is in each case a function of these three variables (see Eq (20)).
In order to acquire the input for the RVE algorithm, the available material sample is split into sets of sub-samples of V k (with L k ¼ V 1=3 k ) of n k realisations. The size function of the RVE is then determined by a regression analysis for the standard deviation D Z(V k ) of Z over V k . In [45], the mean for V k of a physical property Z (V k ) and the respective variation D 2 ZðV k Þ are calculated as follows: Using the expressions gained in Eqs (12) and (13), the absolute sampling error abs and rel are defined for this algorithm, following [29,42,67]: V is linked to the variation D 2 ZðVÞ and the point variance D 2 Z [29, 42, 68]: Z fibre takes in the present study the value of the respective physical property in V fibre , Z void of the property in V void (with E void = 0). Eq (16) requires from the sample volume V [29]: The term integral range A 3 [69][70][71][72] provides a mathematical measure for the characterisation of random structures.

FE modelling of experimental in-plane Young's modulus
The transverse Young's modulus E T (¼ 1 2 ðE x þ E y Þ) is obtained through the in-plane load cases i = 1 and 2 by Model A and B. The experimental results of [36] are approximated by Model B.
Model A-Rigid joint model. The value of h BC has a non-negligible influence on the obtained E T in particular for h BC < 77.50 μm (see Fig 3 for the exemplary values of sample 316L-10%-No.1, and 316L-20%-No.1). Whether h BC still changes E T for h BC ! 77.50 μm depends on the imposed BC type. The deformation plots demonstrate for MBC that for low h BC the applied BC cause a mechanical response almost exclusively in the close proximity to the two constrained cube faces F x0 and F x1 (see Fig 4). For greater h BC , a nearly linear displacement increase through the sample along the axis of tension can be observed in analogy to the increased E T in Fig 3. The change from linear to quadratic interpolation decreases the E T obtained by Timoshenko beams only marginally (see Tables 4 and 5). E T increases considerably when Euler-Bernoulli beams are simulated: The reduced number of prescribed degrees of freedom (DOF) in the case of MBC reduces also the value of E T when compared to KUBC: The findings in Eqs (22) and (23) hold true irrespectively of the h BC value and confirm the prediction of [62] and [29] respectively.
Model B-Spring joint model. E T is obtained by Model B for values of s defined in Table 2 with h BC = 77.50 μm. This h BC value is chosen because of the obtained E T magnitude under MBC (similar to the experimental values of [36]) and because of the considerably reduced dE T /dh BC at this value. As required, results for Model B tend towards Model A for greater s independently of the applied BC (see Fig 5): Under MBC with only two constrained cube faces, E T results match better the experimentally obtained stiffness values of in-plane tensile testing of [36]. The best match of Model B under MBC is obtained for s = 5 μm (see Fig 6) which is used for the following analyses.

Stiffness matrix and transverse isotropy
The dominantly in-plane fibre orientation of the network samples (see [36,73] for detailed analyses) influences the mechanical behaviour to a great extent when the remaining load cases i = 3 to 6 (see Fig 7 for sample 316L-20%-No.2) are obtained. It can be described as transversely isotropic.
Eqs (25)-(27) present the elastic moduli together with the confidence intervals [C ij ± 2D C ij ] under KUBC for B32 of each f. All non-bold entries would be expected to equal zero. However, they only nearly vanish which is caused by the numerical artefacts. C 11 deviates from C 22 by a maximum error of 6.13% for f = 10%. The conditions that C 44 % C 55 and C 13 % C 23 show the biggest errors for f = 15% (5.67% and 4.80% respectively). The obtained error for C 66 % 1 2 ðC 11 À C 12 Þ varies between 3.20% for f = 10% and 11.27% for f = 20%. GPa ð25Þ  GPa ð27Þ Size effect and prediction of RVE size The sample set of V k and respective n k in Table 6 with k between 1 and 6 was extracted randomly for each f from the network samples. n k replicates the scheme found in [45]. The original publication [29] relates to the measurement error which would be rather difficult to quantify for the presented meshed fibres network systems. E T shows for the two BC types two different size effects, decreasing for KUBC and increasing for MBC (see Fig 8). These two patterns can be found similarly in [29], there in particular for the bulk modulus under MBC, and for KUBC in [32].
The application of the regression algorithm in Eq (21) (with Z = E T ) leads to the values of Table 7. Greater f reduces V RVE,E T . The influence of the BC type on the material is found to be as predicted by [29]: For f = 10% under MBC, the condition of Eq (18) [29] isn't fulfilled any more. The precision of the values for this particular simulation point could be investigated further in future work by greater sample sets.

Prediction of fibre deformation mechanism
The fibre segment deformation inside the network sample by deflection w Ã (perpendicular to the fibre segment axis) is compared to the elongation ΔΛ (along the axis) by the median of their ratio under KUBC (see Fig 9). For greater simplicity, load cases i = 1 and 2 are summarised to T , i = 4 and 5 to the out-of-plane shear deformation shear-z . Load case i = 3 stands for tensile-z and i = 6 for shear-xy . Increased f has been shown to reduce the fibre segment length λ (see Table 1 and [36]) which reduces simultaneously in the present study the obtained median [w Ã /ΔΛ]. The obtained median ratio is in each case greater than 2, increasing for shear and out-of-plane deformation:

Conclusions and Future Work
The present study simulates the elastic mechanics of metallic fibre networks by two FE models (rigid inter-fibre joints or torsional spring connectors) and validates the results by experimentally obtained values from [36]. While [29,32] used solid volume meshes, beam theory is shown to offer a valid simplification method for the investigated material. BC types MBC and KUBC are simulated and a new model parameter h BC is introduced for the depth to which BC are prescribed into the material along n in on @V. MBC with fewer prescribed DOF provide a greater match to the experimental values of in-plane tensile testing than KUBC do. The results for rigid inter-fibre joints and for torsional spring connectors converge incrementally. During this process, the influence of the spring constant on the overall  structural stiffness decreases. This implies for the manufacturing process that the stiffness of inter-fibre joints becomes increasingly negligible as soon as a minimum joint strength has been achieved in the sintering process. In [29], an isotropic material behaviour was obtained for computer generated fibre networks without preferred fibre orientation direction. The material of the present study exhibits a dominantly in-plane fibre orientation [36] and an in-plane transversely isotropic mechanical behaviour. Depending on the applied BC type, a decreasing or increasing size effect is observed. Fibre segment deflection dominates the deformation inside the network over fibre segment elongation. These three findings imply for future work and future design studies that size and shape of the material are of importance, as well as the manufactured fibre orientation.
Of interest for future work is also the achievable strain magnitude in a matrix material located in the network's void phase and the material's suitability for mechanical bone growth stimulation [74].