Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 11109.
Published online 2017 September 11. doi:  10.1038/s41598-017-07731-6
PMCID: PMC5593920

Mechanical bone growth stimulation by magnetic fibre networks obtained through a competent finite element technique


Fibre networks combined with a matrix material in their void phase make the design of novel and smart composite materials possible. Their application is of great interest in the field of advanced paper or as bioactive tissue engineering scaffolds. In the present study, we analyse the mechanical interaction between metallic fibre networks under magnetic actuation and a matrix material. Experimentally validated FE models are combined for that purpose in one joint simulation. High performance computing facilities are used. The resulting strain in the composite’s matrix is not uniform across the sample volume. Instead we show that boundary conditions and proximity to the fibre structure strongly influence the local strain magnitude. An analytical model of local strain magnitude is derived. The strain magnitude of 0.001 which is of particular interest for bone growth stimulation is achievable by this assembly. In light of these findings, the investigated composite structure is suitable for creating and for regulating contactless a stress field which is to be imposed on the matrix material. Topics for future research will be the advanced modelling of the biological components and the potential medical utilisation.


Materials incorporating a fibre structure have been an essential constituent in various fields of application for decades1,2. Today, a range of several distinguishable categories of fibre materials exists: polymeric non-woven fabrics3,4, actin networks and cytoskeletons in the field of biomaterials5,6, or also theoretical computer generated network materials79. The mechanical behaviour912, the thermal conductivity9,13, or the magnetic response14,15 are amongst the investigated physical properties.

The purpose of this present study is to investigate the suitability of a specific fibre-matrix composite material for its suitability as mechanically active tissue engineering scaffold. This fibre material is manufactured by sintering randomly stacked stainless steel fibres1619. Its basic properties have been investigated by several studies. It has been shown that beam theory offers an elegant simplification for numerical studies about its mechanics. Further it is known about this material that greater fibre volume fraction f reduces inside the material the average length λ of fibre segments between sintered inter-fibre bonds (Table 1). This has a strong effect on the local mechanical response of the material which is dominated by fibre deflection over fibre elongation12. The deformation under magnetic actuation of individual fibres has been investigated15. This study now extends the scope of the work to entire network geometries. Numerical studies offer a very useful tool in this context. The highly delicate geometry of this network material can be studied numerically at a great level of detail. Local behaviour of the material can be quantified. Those local results would equally be obtainable from experimental measurements only under a great expense of resources. Whenever available, the corresponding experimental measurements allow the calibration of the numerical models.

Table 1
Network samples.

Unlike previous work about this specific fibre material, the present study combines it with a matrix material to a fibre-matrix composite. Static, conventional applications of fibre-matrix composites have been proposed for composites of bioglass20,21 or cordierite22,23. In the field of biomedical engineering, metallic24,25 and since then also non-metallic26,27 fibre networks have been tested for their ability to improve bone cement mechanics. The exhibited properties document very clearly improvements from the composite structure for this orthopaedic application. In non-static or smart applications of fibre-matrix composites, material properties are amended during the application process. Magnetic nanoparticles combined with cellulose fibres are in the discussion for the design of data storage applications and magnetographic printing or filtering14,28. This work has been extended also to incorporate bacterial cellulose29,30. In a different context, magnetic fibre networks could be applied to design heat exchangers of variable drag31,32. The assembly of this present study has been discussed for its suitability as scaffold in tissue engineering33. In this case, the purpose of the fibres would be to deliver under magnetic actuation a mechanical stimulus for the enhanced growth of bone cells. Various other fibre network materials are also investigated for their application as tissue scaffolds29,34. The delivery of a mechanical stimulus from scaffolds for influencing cells and their development is an on-going research topic on tissue level and on organ level3537. For bone cells, it is known that cyclical loadings at 1 Hz of the strain magnitude 0.001 are beneficial for the tissue’s remodelling behaviour38. A beneficial effect of a stand-alone magnetic field on bone growth is equally known39,40.

Yet, the mechanical interaction between fibre phase and matrix phase in fibre-matrix composites is so far only insufficiently understood. The purpose of the present study is to make a contribution to this field and to investigate the mechanical response of a matrix inserted in the void phase of a fibre network under magnetic actuation. Previous work has already analysed for this specific material simplified single-fibre geometries under magnetic actuation15. Or global values have been predicted analytically for fibre assemblies41,42. This present study investigates complete fibre network geometries and analyses the matrix strain on local level. The results of this study are obtained by means of experimentally validated finite element (FE) models and lead to an analytical model for the local matrix strain in cube shaped fibre networks under magnetic actuation. The presented results were part of a dissertation project at the University of Cambridge19.

Mathematical notation

The following notation is used in this document. Scalars are given as x, vectors as x¯, and 2nd rank tensors as x¯¯. The vector product “×” and dot product “·” are applied. Relations which are greater-than and approximately equal are indicated as “[greater, similar]”. All fibre network nodes n of one sample i form the node set N i. The total number of nodes in that sample i is defined as Nˆi. Equality between two volumes is written “ [left and right double arrow ] ”, a volume containing a volume subset is given as “[subset or is implied by]”. The operator for geometry definition by the intersection of two bodies is “∩”. 3D spheres are specified for midpoint and radius as O (midpoint,R). Their respective boundary is written o (midpoint,R) = [partial differential]O (midpoint,R), geometric cube faces are written F.

Modelling Methods and Material

For the investigations of the present study, a FE model (Fig. 1) of linear elasticity is developed. It describes the matrix strain under magnetic actuation in interaction with the fibre material. The quantitative evaluation of local fibre network density and local strain fields follows the scheme of Eqs 10 to 13 and leads to the analytical model presented by this study (Eq. 16).

Figure 1
FE model: (a) Sample dimensions with cube face F x, (b) boundary conditions, (c) FE mesh.

Material samples and fibre network geometry extraction

FE modelling of the fibre network follows steps of an approach based on beam theory12. Three cube shaped fibre network samples are used in this present study (Fig. 1a and Table 1). This allows the investigation of three different values of fibre volume fraction f: 10, 15, and 20%. The influence of f on the mechanical response is of particular interest.

These networks are produced by N.V. Bekaert S.A. (Belgium) from American Iron And Steel Institute (AISI) 316 L (d  40 μm) or by Nikko Techno Ltd. (Japan) from AISI 444 steel fibres, assembled as fibre stacks, compressed, and sintered to fibre mats. During the sintering step, inter-fibre bonds form. The sample sections are cut by electronic discharge machining from the sintered mats and computed tomography (CT) scans are acquired by General Electric (Germany) for a resolution of res = 7.75 μm. Due to the available computing capacities, the sample cube volume is set to V = (775.00 μm)3. In the CT scans, this sample cube side length is the equivalent of 100 pixels.

A skeletonisation algorithm4345 is applied and returns as output the medial axis paths of the fibre bodies contained in the CT data. Resulting from the chosen CT scan resolution res = 7.75 μm, the location of all medial axes is defined inside the sample cubes on a [7.75 μm · 100]3-grid. Nodes defining network medial axis location are referred to as network nodes in the following. Their total number per sample Nˆi is known to increase for greater f 17,18 and is given in Table 1 for the three network samples.

FE meshing

While previous studies have simulated fibre networks for a void inter-fibre space9,12, the present study combines a fibre phase V f and a matrix phase V M in the sample volume V. The meshing of both phases and the defined interaction rules are discussed in the following (Fig. 1c). The mesh is implemented for the FE solver Abaqus 6.1346 (Tables 2 and and3).3). By their definition in the present study, V M equals V while V f is only a subset:

V ⇔ VM ⊂ Vf

Table 2
Implemented FE types.
Table 3
FE modelling parameters.

Eq. 1 simplifies the volume description by disregarding the marginal volume of fibres whose medial axis is located at a distance to the surface S of less than d/2.

Matrix phase

A mesh of regular hexahadra elements (C3D8R) fills V. These hexahedra represent V M. For simplifying the interaction with the fibre bodies, the hexahedra side length L hex is set to the CT scan resolution res = 7.75 μm = 1 pixel. The number of hexahedra in V follows as N hex = 1003, their volume as V hex = (7.75 μm)3.

The values of matrix stiffness E M are chosen as it would be applicable for three stages in bone tissue development: collagenous bone47, granulation tissue48, and immature bone49. The value for immature bone was chosen from the early phase of bone development. The value for collagenous bone is slightly greater than found in the literature for enhancing the FE solver convergence.

Beam assembly

The fibre medial axes are transferred in the present study into beam assemblies. Two Timoshenko beams5052 are implemented. Results for linear interpolation (B31) and quadratic interpolation (B32) are analysed and compared to predictions from previous work12. The beams are simulated for a simplified round cross section of diameter d = 40 μm (i.e. A f = (d/2)2 π), a Young’s modulus of E f = 200 GPa, and for a Poisson’s ratio of ν = 0.3. The sintered inter-fibre joints are implemented by torsional springs (CONN3D2). Their stiffness is simulated as K joint = sE f A f with s = 5 μm. This value of s has been found to a pproximate results from experimental measurements12,17. The number of network nodes in Table 1 refers to exclusively those nodes which define beam start points and end points. The implementation for the FE solver requires additional technical nodes. Those purely technical nodes provide additional integration points for B32 elements, indicate the normal relative to the beam axis, and one is needed for each CONN3D2 element.

Interaction rules

The implementation of the interaction between V fibre and V matrix (Fig. 1c) is based on a previously published model for osteosynthesis applications53. The regular hexahedra and the beam mesh share the 1003-node grid which is defined by the CT scan pixels. Three main simplifications are made in the modelling:

  • V f and V M overlap, a double section assignment exists for up to max(f) = 20% of V.
  • The interaction between both phases is reduced to the shared node grid.
  • The mechanical interaction is set to identical displacement u¯ at each node.

This means that more complex mechanical interaction as it appears in biological tissue between cells and scaffold is not considered at this point. These simplifications are made while, as in the original model, the following holds:

Ef ≫ EM

FE simulation and BC

In the case of a void inter-fibre space, the Cauchy stress tensor σ¯¯ is defined for a fibre network consisting of V f and the void matrix phase V void 9:


In this study, the relationship is modified for V M. Eq. 4 defines now the stress fields in the fibre-composite as:


This relationship might seem obvious. It is mentioned explicitly at this point because Eq. 4 is required for the evaluation of imposed matrix strain with regard to the magnitude of E M and with regard to the respective magnetic actuation.

Magnetic actuation and BC

In the simulations, a magnetic induction vector B¯ is imposed (Fig. 1b). The ferromagnetic mechanical response of the fibres is modelled as moment vector τ¯B. It is calculated for each beam element of length L f under the assumption of complete magnetisation parallel to the beam axis15:


In Eq. 9, τ¯M of all fibres in the respective sample is added into the equilibrium condition of moments. The actuation model is applied for M s = 1.6 MA/m. This value of M s has been experimentally validated for ferritic AISI 44415. The geometries of the network samples were obtained for austenitic AISI 316L. As part of this study’s modelling assumptions, the network is simulated for the hypothetical M s of AISI 444. It has been shown that the geometry of the AISI 316L networks can similarly also be produced by AISI 444 fibres17,18. The imposed magnetic induction vector B¯=[100]B is orientated along the x-axis. At the time of this study’s design, experiments were carried out at the University of Cambridge in which this type of in-plane actuation was tested.

The sample surface S of the sample volume V is defined by six quadratic cube faces, two of them perpendicular each to one of the three geometric axes. The cube face F x, marked in red in Fig. 1a, is used for the definition of the boundary conditions (BC). All three translational degrees of freedom (DOF) of beam elements located on the cube face F x are kinematically constrained up to a depth of h BC = 77.50 μm into the material. This value of h BC was applied as it provides a match to the experimental stiffness magnitude12,17.

Equilibrium conditions

The FE solver is run under the assumption of linear elasticity and for the equilibrium conditions of forces (Eq. 6) and moments (Eq. 8). It is required that div(σ¯¯)=0 is fulfilled. Body force per volume f¯ is set to 0 in the present study; gravitational forces being a typical example for f¯. Together with the surface traction vector t¯, it defines the equilibrium of forces. Surface traction itself is defined by σ¯¯ and the outer-normal vector: t¯=σ¯¯n¯out 55.


As the magnetic response is modelled as moment vector τ¯B, the general condition for equilibrium of forces can be simplified for the present study. Considering that f¯=0 it follows that:


Relative to the origin of point vector x¯, the general form of the equilibrium of moments55 is defined:


The simplification of Eq. 6 by f¯=0 can be repeated for Eq. 8. The sum of τ¯B which is imposed on the fibre volume V f is added. The modified equation representing the present study and considering the sum of imposed τ¯B follows as:


The Lagrangian reference frame56 is used for Eqs 6 to 9. The alternative Eulerian reference frame can be advantageous for the modelling of fluids57,58. During the study, the solver runs into up to eight singularities and one zero-pivot (Table 1). For the facilitation of the solving step, singularities and zero-pivots are constrained for their respective DOF. Due to the size of the beam mesh (>103) this number of additionally constrained nodes is negligible with regard to the presented mechanical analyses.

Network node counting and strain field averaging methods

Eqs 10 and 11 are applied in this present study for further analyses of local fibre network density. For that purpose, a centre point CP=[0L/2L/2]=[0387.50387.50] μm is defined on the constrained cube face F x (Fig. 2a). The function y (R) stands for the relative share of nodes in sample volume V which is located on the sphere boundary o (CP,R) of radius R.


Y(R)i=VO(CP,R)Ni/N^i=0Ry(R)dRi{10%,15%,20%}:samples of 10, 15, and 20% f

Figure 2
Averaging methods: (a) Averaging method g (R) around CP on constrained cube face F x (Eq. 12) and (b) f (R) around fibre network (Eq. 13).

Y (R) stands for the integral of y (R) over R and counts network nodes of V inside the sphere O (CP,R). y (R) and Y (R) lead to g (R) which is the first of two applied strain averaging methods. Both strain averaging methods, g (R) and f (R), quantify the local strain obtained in the matrix phase V M. Two equivalent strain measures, von-Mises strain measure ε v.M. and maximum principal strain ε Pmax 59,60, are calculated from the 3D strain fields of V M. While ε v.M. is generally recommended for ductile material, ε Pmax is of greater precision for brittle materials61. The material behaviour type of bone has been shown to depend on the deformation rate /dt 62.

g (R) quantifies the matrix strain on the boundary o (CP,R) where it intersects with V (Fig. 2a and Eq. 12). f (R) places in a first step a sphere O (n,R) on each matrix node n and averages the matrix strain within O (n,R). In the second step, the average is obtained across the entire node set N i (Fig. 2b and Eq. 13). g (R) quantifies the matrix strain depending on the distance to the constrained cube face. f (R) quantifies the matrix strain depending on the average distance to the fibre network structure.


equation M77


For the solving step, the FE code is run on the Cambridge High Performance Computing Cluster Darwin. The cluster provides a total of 9600 cores by 600 quad server Dell C6220 chassis. Two 2.60 GHz eight core Intel Sandy Bridge E5-2670 processors with sixteen cores in total form one node with 64 GB of RAM (4 GB per core)63.


The matrix phase of the meshed fibre-matrix composite (Fig. 1) is analysed for its deformation. Results about deformation patterns and deformation magnitude under magnetic actuation are presented. Distinctive patterns in the local matrix strain distribution are discovered. The results lead to an analytical model for the matrix strain magnitude in cube shaped fibre-matrix composites. A particular focus is the influence of the local fibre network density.

Local fibre network density

As explained above, the analysed fibre material is a statistical material. This requires for mechanical modelling knowledge about the distribution of fibre density and about the distribution of fibre orientation. In particular, this knowledge is of importance for the discussion of the local deformation patterns discovered in the present study. Previous analyses have shown that the main fibre orientation direction lies in the xy-plane; in the xy-plane itself, no preferred direction of orientation can be seen17. Figure 3 and Table 1 of the present study analyse now further the local distribution of material density across the three samples by quantification of the network node distribution. The density of network nodes relates directly to the fibre density in V.

Figure 3
Local network density: (a) to (c) network node distribution along Cartesian axes x, y, z, and (d)/(e) network node distribution on/in sphere O (CP,R).

The main result is that the variation of node density along the three Cartesian axes is far greater than the density variation relative to the sphere O (CP,R). The network nodes are counted per 2D pixel slice of a 253-grid along each of the three Cartesian axes. Table 1 shows the obtained numerical average, the median, and the standard deviation. A visible deviation exists in each of the three samples along each Cartesian axis. The maximum standard deviation relative to the average exists in min(f) = 10% along y and z. These two findings are very interesting as their cause lies in the compression step during manufacturing. Future network design might be able to exploit this further. However, this would make necessary a degree of control about production parameters greater than the one available in today’s process. The corresponding distributions are plotted in Fig. 3a to c together with each sample’s average. These three plots do not exhibit distinct local sample sections of outstanding greater/lesser material density. However, the extent of density variation along the three Cartesian axes becomes visible in these plots. The plot of Sample-10% shows clearly the greatest density variation relative to the sample average. Since the magnetic induction is imposed on the matrix material by the fibre network (Eq. 5), the rather uneven material density distribution along the three axes leads to local variations of obtained matrix strain (Figs 4 and and5b5b).

Figure 4
Deformation plot: ε Pmax of Sample-15% for beam B31 under magnetic actuation, green range on colour bar scaled to magnitude 0.001.
Figure 5
Matrix strain: (a) magnitude ε v.M. and ε Pmax depending on f, (b) spatial distribution of ε v.M. in immature bone, (c) statistical ε v.M. distribution.

The share of nodes located on/in O (CP,R) is plotted in Fig. 3d and e as function y (R) and Y (R) (Eqs 10 and 11). Both, y (R) and Y (R), exhibit a variation of the density distribution along R which is of far lesser magnitude than the one along the three Cartesian axes. This holds for each of the three samples. In Fig. 3d, y (R) is plotted. The three samples follow the same pattern of three distinct sections along R which can be seen in the distribution:

  • R  387.5 μm: y (R) increases exponentially. O (CP,R) corresponds to the smallest sphere drawn in Fig. 2a and is entirely contained in V.
  • 387.5 μm < R  775 μm: O (CP,R) exceeds the extensions of V at the cube’s sides and y (R) decreases, in a good approximation linearly.
  • 775 μm < R: O (CP,R) intersects with V only at the four cube corners opposite to cube face F x as drawn for the greatest sphere in Fig. 2a. y (R) drops sharply to 0.

In Fig. 3e, it can be seen that the value of Y (R) increases for greater R continuously at different rates. A node share of 100% is reached when the entire V is contained in O (CP,R), i.e R > 950 μm. Very important for the context of this study, the patterns are independent of f (i.e. it is obtained in each of the three samples). The findings of Fig. 3d and e are new for this specific material and are used in this present study for deriving the analytical model of the matrix strain magnitude (Fig. 6b and Eq. 16).

Figure 6
Matrix strain: (a) Distribution around fibre network and (b) around CP on cube face F x.

Matrix strain magnitude and distribution

The FE assembly is simulated for the three sample geometries. The obtained magnitude of matrix strain and its distribution inside the sample, as well as, its statistical distribution are discussed now in this section.

Deformation plot

Figure 4 contains the deformation plots of Sample-15% under magnetic actuation for all three matrix materials defined in Table 3. The green range on the colour bar is scaled to the magnitude 0.001 because of the particular interest with regard to applications including biomechanical bone growth stimulation38. Each of the nine plots contains locally this specific magnitude of matrix strain. The required magnitude of B is within the experimentally achievable range. Up to B = 1.00 can be seen as practical. This finding is of importance because it suggests the general suitability of this specific network material as scaffold for bone growth stimulation under magnetic actuation. Local maxima and minima are visible on S. For an increase from B = 0.25 to 1.00, the extension of the green areas increases, spreading further into the material. Greater E M reduces them respectively. The quantitative analysis of the matrix strain magnitude and distribution follows in Fig. 5.

Strain magnitude

Figure 5a plots the strain magnitude when averaged across the sample volume V for the previously shown interesting value of B = 1.00. The mechanical response is scaled linearly by B due to the modelling assumption of linear elasticity. The values of ε Pmax for brittle behaviour are marginally smaller than the values of ε v.M. for each value of f, and for each E M. This relationship is also predicted in the literature59 and holds independently of the Timoshenko beam (Eq. 15). Lesser E M increases the matrix strain non-linearly. Of interest is the result that the average matrix strain clearly exceeds 0.001 in the case of collagenous bone and granulation tissue. At this early stage of bone tissue development, lesser values of B are sufficient. For immature bone and B = 1.00, average matrix strain between 10−4 and 0.001 is obtained. A greater density of fibre network material increases the matrix strain for the given size of V. It is known for this network material that its mechanical response is influenced at this scale by a pronounced size effect12. Whether and how this size effect also influences the matrix strain under magnetic actuation is at this point one very worthwhile topic of future research.

The two Timoshenko beams lead to almost identical results for the matrix strain magnitude (Table 4). B32 is known to marginally reduce the network’s mechanical stiffness12:

EB31 ≳ EB32

Table 4
Timoshenko beam.

Conversely, under magnetic actuation the quadratic Timoshenko beam B32 leads to marginally greater matrix strain than linear B31 if B = const. This relation holds for both equivalent strain measures.

ifB=const,EM=const,f=const, it holds thatεPmaxεv.M., valid for both Timoshenko beams (B31 or B32)εB31εB32, valid for both equivalent strain measures (εv.M. or εPmax)

The obtained difference Δ is below the range of the magnitude of Δ is nearly independent of B. The obtained variation for B is in the range of the computing precision. The influence of f is not definite. For practical considerations in future prototype studies, the difference between both Timoshenko beams is negligible, especially with regard to other errors such as the one of the skeletonisation step. It is interesting that the difference increases for greater E M. One possible reason is that the stiffness requirement between matrix and fibres (Eq. 2) is better fulfilled for lesser E M. For numerical considerations, a further investigation of this model aspect through additional sample geometries will be highly worthwhile.

Spatial distribution

The spatial distribution of matrix strain (plotted on a 52-grid and averaged along the z-axis) shows a comparable pattern between the three samples including at the same time local variations (Fig. 5b). The main finding is that the matrix strain magnitude is not uniform across the samples under magnetic actuation. For uniaxial mechanical actuation, the samples show a very uniform deformation12. Under magnetic actuation in this present study, the minimum matrix strain can be found in the proximity of the kinematically constrained cube face F x at min(x). Matrix strain magnitude increases with proximity to the non-constrained cube faces. In each of the samples, a maximum strain peak is located at the grid cell of max(x) and max(y). Inside the sample V, maxima of matrix strain exist too in each case. In the case of Sample-15%, a region of greatest strain on the cube’s top face visible in the plots of Fig. 4 corresponds to the local maximum in Fig. 5b. Consideration of the results in Fig. 5b will be mandatory for future prototype designs. The strain magnitude and its location clearly depend on the applied BC and local variations exist. The design of BC imposed on the assembly will have to be considered in future studies. This is found to be a novel and reproducible feature of this specific material.

Statistical distribution

The statistical distribution function p (x) of ε v.M. is plotted in Fig. 5c and exhibits in each sample a single global p (x)-peak. This distribution pattern is obtained for each of the three samples. A change from B = 0.25 to 1.00 shifts this single global p (x)-peak towards greater values of ε v.M.. This shift demonstrates that for greater B the matrix strain increases uniformly inside the sample. Conversely, greater E M shifts the global p (x)-peak towards smaller values of ε v.M.. The magnitude of ε v.M.  0.001 required for bone growth stimulation38 can be found locally in each sample for B = 0.25 and 1. This finding is mandatory for the intended application as tissue engineering scaffold. A localization of strain magnitude in the sample volume follows in the next section below.

Matrix strain relative to fibre network and CP

The by far most relevant findings for the intended application purpose of this present study are shown in Fig. 6. Based on the strain distribution functions g (R) and f (R) (Eqs 12 and 13), it is possible to quantify and to analytically localize areas of greater or lesser matrix strain magnitude inside the material.

Figure 6a plots f (R) which is the average matrix strain in O (n,R) depending on R averaged for all network nodes n of the sample. The obtained results show that R has a strong influence on f (R). A specific pattern is discovered independently of f, E M, or B. In the close proximity to the fibre structure (R  0+), f (R) exhibits in each sample one global peak which exceeds the sample average by several magnitudes. For increasing R [set membership] [200 μm, 800 μm], f (R) drops below the sample average. In this range, f (R) is under the influence of the low-strain regions inside the samples and along F x (Fig. 5b). For values of R > 800 μm, the regions of greater matrix strain near the non-constrained cube faces influence f (R) stronger and it tends towards the sample average. This finding is particularly interesting for future research. One topic will be the mechanical interaction between fibres and the matrix at the interphase of these two. Biomechanical materials as used in the potential application of tissue engineering are characterized by non-linear mechanical response and additional parameters such as limited adhesion forces64. Those effects will be particularly important at the interface. Whether and how the predicted strain field changes remains to be seen.

The plots of Fig. 6b show g (R), the matrix strain obtained on a sphere o (n,R) according to the averaging method Eq. 12. The obtained matrix strain increases for greater R and a local maximum is located at around R = 600 μm. This pattern is obtained for each matrix material in each sample. The individual differences of network geometry and material density in the samples do not change this overall pattern. The absolute maximum of strain is obtained in the non-constrained sample corners at R  900 μm for Sample-10% and Sample-20%. In the case of Sample-15%, the sample maximum is reached at R  850 μm. In each of the three samples, the obtained matrix strain drops afterwards. One difference between the samples is that this drop is much more pronounced in the case of Sample-15%. Low magnitude of matrix strain in the non-constrained corners of this sample is visible in the deformation plots (Fig. 4). It is worth mentioning that at this value of R only a fraction of the total 106 matrix hexahedra are located. The intersection of o (n,R) with V for R > 850 μm in the sample corners affects only ≈2% of all matrix hexahedra. This causes the obtained differences which are statistically not representative.

Based on the results of Fig. 6b, it is possible to derive for the first time a simple analytical model for the prediction of local matrix strain inside the sample’s V. Eq. 16 approximates the obtained values by a power function:


For R > 860 μm, the fit between Eq. 16 and the data declines sharply. This is why the regression analysis is run for R  860.25 μm. This means that 9,724 of 106 hexahedra are not included in Eq. 16, ≈1%. The matrix strain of the other 99% is modelled by Eq. 16. Table 5 contains the regression values for a, b, and the coefficient of determination r 2. The values are given for ε v.M., ε Pmax, and for both Timoshenko beams. r 2 varies between 0.82 and 0.94 which indicates a very good fit. The values of a and b are very similar in each case for the two Timoshenko beams. A clear trend of a and b for changes of f can’t be identified at this stage. Whether this is achievable for greater V is one highly interesting topic for future work. A clear result of this study is that the strain magnitude exhibits a specific pattern around CP and that it can be described by this simple analytical model. Differences can be seen between the three samples. The overall pattern repeats.

Table 5
Regression analysis.


We have presented above the results for the deformation of a fibre-matrix composite under magnetic actuation. The local matrix strain has been studied and quantified in depth for the first time. The results show for the investigated assembly distinct deformation patterns under magnetic actuation and conclusions about the influence of model parameters can be drawn. With regard to the general material response, it must be taken from now on into consideration that the matrix deformation is not uniform across the samples under magnetic actuation. This is in contrast to previous work which has been able to show that this particular network material exhibits under uniaxial mechanical actuation a very uniform deformation12.

One additional conclusion from this study is that BC have a direct influence on the resulting deformation field under magnetic actuation. For the cube shaped samples in this study, the local matrix strain reaches its minimum in the proximity of the constrained cube face. Local maxima exist inside the sample volume V. Maxima are located in the free corners of the sample cubes. In future studies, the applied BC will have to be explicitly considered. One major material parameter is the fibre volume fraction f. It can be concluded for this model parameter that the magnitude of matrix strain increases for greater f. Whether this trend will also be obtained for other sample volume V is one very interesting topic for future research. A pronounced size effect of the mechanical response is known for this network material12.

The statistical distribution contains in each sample one single global peak which shifts depending on matrix stiffness and magnetic induction. The conclusion is that controlling the strain magnitude of this incidence peak is mandatory for controlling the average and overall response in possible future applications. At the same time, also local sections of greater or lesser strain magnitude exist inside the matrix. The results for two derived distribution functions (Eqs 12 and 13) allow to gain novel insights about the localization of matrix strain magnitude inside the material. First, it can be shown that the matrix strain increases in each sample in the proximity of the fibre structure. Second, the distance to the kinematically constrained cube face has a strong influence on matrix strain magnitude and can be used for deriving a simple analytical model of matrix strain magnitude (Eq. 16). This model is valid for 99% of the matrix volume. The achieved fit for the regression analysis is a very good match with r 2 = 0.82 or greater. A remaining challenge are the free, unconstrained corners of the cube.

With regard to the purpose of this study to investigate the suitability of the fibre-matrix composite for possible applications as bioactive scaffold in tissue engineering, one very important result is obtained. From experimental work, it has been known since 1985 that a cyclical strain of 0.001 of the frequency 1 Hz is beneficial for the stimulation of bone growth38. This strain magnitude must be considered as design target for the scaffold application. The results in this study are showing that this magnitude is achievable by this investigated fibre-network composite assembly. In conclusion, it is possible to say that the principal suitability of this material for bone growth stimulation under magnetic actuation is confirmed. Future work will have to relax the simplifications made for the matrix model. For the application as tissue engineering scaffold, understanding of local strain at the interface between cells and steel is mandatory. The advanced version of this current model could integrate biological aspects such as the detaching of focal adhesion points build by bone cells on the fibre structure64. Also a dynamic matrix model including the active migration of bone cells through the fibre network65 after the cell seeding would be a worthwhile aspect and of great interest for understanding bone development. The current model still relies on a static mesh.

In future work, one further effect must not be neglected in the evaluation of bone growth. It is known that magnetic fields alone already can stimulate bone growth. It is applied in the treatment of avascular necrosis or pseudarthrosis39. Designs for hip joint prostheses make use of that effect40. If benefits are achieved by the proposed assembly it will have to be made clear what the actual cause of the improvement is. The strain field imposed by the scaffold or the direct reaction of the tissue to a simple magnetic field? In the ideal case, it will be possible to raise synergies from both effects.


The authors wish to thank Dr. Arul Britto from CUED for his advice about model implementation in Abaqus and Dr. Stuart Rankin from the Cambridge University High Performance Computing Division for his advice about the Darwin Cluster. Prof Rainer Bader and Dr. Jan Wieding are acknowledged for the very useful discussions about their steel-bone interaction model. This research was supported by the European Research Council Grant No. 240446.

Author Contributions

Author Contributions

W.A.B. designed the experiments, W.A.B. conducted the experiment, W.A.B. analysed the results, W.A.B. wrote the manuscript.


Competing Interests

The author declares that they have no competing interests.


Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


1. Ducheyne P, Aernoudt E, De Meester P. The mechanical behaviour of porous austenitic stainless steel fibre structures. Journal of Materials Science. 1978;13:2650–2658. doi: 10.1007/BF02402752. [Cross Ref]
2. Cox H. The elasticity and strength of paper and other fibrous materials. British Journal of Applied Physics. 1951;3:72–79. doi: 10.1088/0508-3443/3/3/302. [Cross Ref]
3. Mannarino MM, Rutledge GC. Mechanical and tribological properties of electrospun PA 6(3)T fiber mats. Polymer. 2012;53:3017–3025. doi: 10.1016/j.polymer.2012.04.039. [Cross Ref]
4. Farukh F, et al. Mechanical analysis of bi-component-fibre nonwovens: Finite-element strategy. Composites Part B: Engineering. 2015;68:327–335. doi: 10.1016/j.compositesb.2014.09.003. [Cross Ref]
5. Lieleg O, Schmoller KM, Claessens MMAE, Bausch AR. Cytoskeletal polymer networks: viscoelastic properties are determined by the microscopic interaction potential of cross-links. Biophysical journal. 2009;96:4725–4732. doi: 10.1016/j.bpj.2009.03.038. [PubMed] [Cross Ref]
6. Yamaoka H, Matsushita S, Shimada Y, Adachi T. Multiscale modeling and mechanics of filamentous actin cytoskeleton. Biomechanics and modeling in mechanobiology. 2012;11:291–302. doi: 10.1007/s10237-011-0317-z. [PubMed] [Cross Ref]
7. Head DA, MacKintosh FC, Levine AJ. Nonuniversality of elastic exponents in random bond-bending networks. Physical Review E. 2003;68:025101(R). doi: 10.1103/PhysRevE.68.025101. [PubMed] [Cross Ref]
8. Shahsavari AS, Picu RC. Size effect on mechanical behavior of random fiber networks. International Journal of Solids and Structures. 2013;50:3332–3338. doi: 10.1016/j.ijsolstr.2013.06.004. [Cross Ref]
9. Dirrenberger J, Forest S, Jeulin D. Towards gigantic RVE sizes for 3D stochastic fibrous networks. International Journal of Solids and Structures. 2014;51:359–376. doi: 10.1016/j.ijsolstr.2013.10.011. [Cross Ref]
10. Picu RC. Mechanics of random fiber networks - a review. Soft Matter. 2011;7:6768–6785. doi: 10.1039/c1sm05022b. [Cross Ref]
11. Rodney D, Gadot B, Martinez OR, du Roscoat SR, Orgéas L. Reversible dilatancy in entangled single-wire materials. Nature Materials. 2016;15:72–77. doi: 10.1038/nmat4429. [PubMed] [Cross Ref]
12. Bosbach WA. The Elastic Behaviour of Sintered Metallic Fibre Networks: A Finite Element Study by Beam Theory. Plos One. 2015;10:e0143011. doi: 10.1371/journal.pone.0143011. [PMC free article] [PubMed] [Cross Ref]
13. Wang, M., Kang, Q. & Pan, N. Thermal conductivity enhancement of carbon fiber composites. Applied Thermal Engineering29, 418–421.
14. Small AC, Johnston JH. Novel hybrid materials of magnetic nanoparticles and cellulose fibers. Journal of Colloid and Interface Science. 2009;331:122–126. doi: 10.1016/j.jcis.2008.11.038. [PubMed] [Cross Ref]
15. Markaki, A. E. & Justin, A. W. A magneto-active scaffold for stimulation of bone growth. Materials Science & Technology30, 1590–1598.
16. Jin MZ, Chen CQ, Lu TJ. The mechanical behavior of porous metal fiber sintered sheets. Journal of the Mechanics and Physics of Solids. 2013;61:161–174. doi: 10.1016/j.jmps.2012.08.006. [Cross Ref]
17. Neelakantan S, Bosbach W, Woodhouse J, Markaki A. Characterization and deformation response of orthotropic fibre networks with auxetic out-of-plane behaviour. Acta Materialia. 2014;66:326–339. doi: 10.1016/j.actamat.2013.11.020. [Cross Ref]
18. Spear RL, et al. Physical and Biological Characterization of Ferromagnetic Fiber Networks: Effect of Fibrin Deposition on Short-Term In Vitro Responses of Human Osteoblasts. Tissue Engineering Part A. 2015;21:463–474. doi: 10.1089/ten.tea.2014.0211. [PMC free article] [PubMed] [Cross Ref]
19. Bosbach, W. A. The mechanical and magnetic behaviour of sintered fibre networks and their suitability for a therapeutic, biomedical application. Ph.D. thesis, University of Cambridge (Date of access: 30/08/2017) (2015).
20. Ducheyne P, Hench L. The processing and static mechanical properties of metal fibre reinforced bioglass. Journal of Materials Science. 1982;17:595–606. doi: 10.1007/BF00591494. [Cross Ref]
21. Boccaccini AR, Ovenstone J, Trusty PA. Fabrication of woven metal fibre reinforced glass matrix composites. Applied Composite Materials. 1997;4:145–155.
22. Kaya C, Kaya F, Boccaccini AR. Fabrication of Stainless-Steel-Fiber-Reinforced Cordierite-Matrix Composites of Tubular Shape Using Electrophoretic Deposition. Journal of the American Ceramic Society. 2002;85:2575–2577. doi: 10.1111/j.1151-2916.2002.tb00499.x. [Cross Ref]
23. Kaya C, Kaya F, Boccaccini AR. Electrophoretic deposition infiltration of 2-D metal fibre-reinforced cordierite matrix composites of tubular shape. Journal of Materials Science. 2002;37:4145–4153. doi: 10.1023/A:1020087819697. [Cross Ref]
24. Taylor D, Clarke F, McCormack B, Sheehan J. Reinforcement of Bone Cement Using Metal Meshes. Journal of Engineering in Medicine. 1989;203:49–53. doi: 10.1243/PIME_PROC_1989_203_007_01. [PubMed] [Cross Ref]
25. Kotha SP, Li C, Schmid SR, Mason JJ. Fracture toughness of steel-fiber-reinforced bone cement. Journal of biomedical materials research. Part A. 2004;70A:514–521. doi: 10.1002/jbm.a.30107. [PubMed] [Cross Ref]
26. Yang J-m, Huang P-y, Yang M-c, Lo SK. Effect of MMA-g-UHMWPE Grafted Fiber on Mechanical Properties of Acrylic Bone Cement. Journal of biomedical materials research. Part A. 1997;38:361–369. doi: 10.1002/(SICI)1097-4636(199724)38:4<361::AID-JBM9>3.0.CO;2-M. [PubMed] [Cross Ref]
27. Puska MA, Närhi TO, Aho AJ, Yli-Urpo A, Vallittu PK. Flexural properties of crosslinked and oligomer-modified glass-fibre reinforced acrylic bone cement. Journal of Materials Science: Materials in Medicine. 2004;15:1037–1043. [PubMed]
28. Marchessault RH, Rioux P, Raymond L. Magnetic cellulose fibres and paper: preparation, processing and properties. Polymer. 1992;33:4024–4028. doi: 10.1016/0032-3861(92)90600-2. [Cross Ref]
29. Vitta, S., Drillon, M. & Derory, A. Magnetically responsive bacterial cellulose: Synthesis and magnetic studies. Journal of Applied Physics108 (2010).
30. Gao X, Sozumert E, Shi Z, Yang G, Silberschmidt VV. Assessing stiffness of nanofibres in bacterial cellulose hydrogels: Numerical-experimental framework. Materials Science and Engineering: C. 2017;77:9–18. doi: 10.1016/j.msec.2017.03.231. [PubMed] [Cross Ref]
31. Tadrist L, Miscevic M, Rahli O, Topin F. About the use of fibrous materials in compact heat exchangers. Experimental Thermal and Fluid Science. 2004;28:193–199. doi: 10.1016/S0894-1777(03)00039-6. [Cross Ref]
32. Golosnoy IO, Cockburn A, Clyne TW. Optimisation of Metallic Fibre Network Materials for Compact Heat Exchangers. Advanced Engineering Materials. 2008;10:210–218. doi: 10.1002/adem.200800021. [Cross Ref]
33. Bosbach, W. A. Finite element modelling of bone growth stimulation and its suitability for therapeutic, biomedical applications (Date of access: 11/05/2017). In von-Behring-Röntgen-Symposium, 1–2, (Giessen (Germany), 2015).
34. Li WJ, Laurencin CT, Caterson EJ, Tuan RS, Ko FK. Electrospun nanofibrous structure: A novel scaffold for tissue engineering. Journal of Biomedical Materials Research. 2002;60:613–621. doi: 10.1002/jbm.10167. [PubMed] [Cross Ref]
35. Mammoto T, Ingber DE. Mechanical control of tissue and organ development. Development. 2010;137:1407–1420. doi: 10.1242/dev.024166. [PubMed] [Cross Ref]
36. Badylak SF, Taylor D, Uygun K. Whole-Organ Tissue Engineering: Decellularization and Recellularization of Three-Dimensional Matrix Scaffolds. Annual Review of Biomedical Engineering. 2011;13:27–53. doi: 10.1146/annurev-bioeng-071910-124743. [PubMed] [Cross Ref]
37. Marx V. Organs from the lab. Nature. 2015;522:373–377. doi: 10.1038/522373a. [PubMed] [Cross Ref]
38. Rubin CT, Lanyon LE. Regulation of bone mass by mechanical strain magnitude. Calcified tissue international. 1985;37:411–417. doi: 10.1007/BF02553711. [PubMed] [Cross Ref]
39. Grunert P, et al. Establishment of a novel in vitro setup for electric and magnetic stimulation of human osteoblasts. Cell Biochem Biophys. 2014;70:805–817. doi: 10.1007/s12013-014-9984-6. [PubMed] [Cross Ref]
40. Schmidt C, Zimmermann U, Van Rienen U. Modeling of an Optimized Electrostimulative Hip Revision System under Consideration of Uncertainty in the Conductivity of Bone Tissue. IEEE Journal of Biomedical and Health Informatics. 2015;19:1321–1330. doi: 10.1109/JBHI.2015.2423705. [PubMed] [Cross Ref]
41. Clyne TW, Markaki AE, Tan JC. Mechanical and magnetic properties of metal fibre networks, with and without a polymeric matrix. Composites Science and Technology. 2005;65:2492–2499. doi: 10.1016/j.compscitech.2005.05.037. [Cross Ref]
42. Markaki AE, Clyne TW. Magneto-mechanical stimulation of bone growth in a bonded array of ferromagnetic fibres. Biomaterials. 2004;25:4805–4815. doi: 10.1016/j.biomaterials.2003.11.041. [PubMed] [Cross Ref]
43. Lee TC, Kashyap RL, Chong-Nam C. Building Skeleton Models via 3-D Medial Surface/Axis Thinning Algorithms. Graphical Models and Image Processing. 1994;56:462–478. doi: 10.1006/cgip.1994.1042. [Cross Ref]
44. Lindquist, W. B. 3DMA General Users Manual (Stony Brook (USA), 1999).
45. Yang, H. A Geometric and Statistical Analysis of Fibrous Materials from Three-Dimensional High Resolution Images. Ph.D. thesis, State University of New York at Stony Brook (2001).
46. Dassault Systèmes. Abaqus 6.13 Online Documentation (Date of access: 11/05/2017), (2013).
47. Engler AJ, Sen S, Sweeney HL, Discher DE. Matrix elasticity directs stem cell lineage specification. Cell. 2006;126:677–89. doi: 10.1016/j.cell.2006.06.044. [PubMed] [Cross Ref]
48. Wang S, Liu GR, Hoang KC, Guo Y. Identifiable range of osseointegration of dental implants through resonance frequency analysis. Medical engineering & physics. 2010;32:1094–106. doi: 10.1016/j.medengphy.2010.07.015. [PubMed] [Cross Ref]
49. Torzilli PA, Takebe K, Burstein AH, Heiple KG. Structural properties of immature canine bone. Journal of biomechanical engineering. 1981;103:232–238. doi: 10.1115/1.3138286. [PubMed] [Cross Ref]
50. Timoshenko S. On the Correction for Shear of the Differential Equation for Transverse Vibrations of Prismatic Bars. Philosophical Magazine Series 6. 1921;41:744–746. doi: 10.1080/14786442108636264. [Cross Ref]
51. Timoshenko S. On the transverse vibrations of bars of uniform cross-section. Philosophical Magazine Series 6. 1922;43:125–131. doi: 10.1080/14786442208633855. [Cross Ref]
52. Dassault Systèmes. Abaqus Analysis User’s Guide 29.3.3 Choosing a beam element (Date of access: 11/05/2017), (2013).
53. Wieding J, Souffrant R, Fritsche A, Mittelmeier W, Bader R. Finite Element Analysis of Osteosynthesis Screw Fixation in the Bone Stock: An Appropriate Method for Automatic Screw Modelling. PloS ONE. 2012;7:1–10. doi: 10.1371/journal.pone.0033776. [PMC free article] [PubMed] [Cross Ref]
54. National Institute for Subatomic Physics. Materials Technical Specification (Date of access: 28/10/2014) (2003).
55. Dassault Systèmes. Abaqus Theory Guide 1.5.1 Equilibrium and virtual work (Date of access: 11/05/2017), (2013).
56. Logan, J. D. Applied Mathematics 6 edn (Lincoln, Nebraska (USA) 2006).
57. Atzberger PJ. Stochastic Eulerian Lagrangian methods for fluid-structure interactions with thermal fluctuations. Journal of Computational Physics. 2011;230:2821–2837. doi: 10.1016/ [Cross Ref]
58. Azarnykh D, Litvinov S, Adams NA. Numerical methods for the weakly compressible Generalized Langevin Model in Eulerian reference frame. Journal of Computational Physics. 2016;314:93–106. doi: 10.1016/ [Cross Ref]
59. Hosford, W. Mechanical behavior of materials 2nd edn (New York (USA), 2010).
60. Alting, L. Manufacturing engineering processes 2nd ed (Lyngby (Denmark), 1994).
61. Gross, D., Hauger, W., Schröder, J., Wall, W. A. & Bonet, J. Engineering Mechanics 2: Mechanics of Materials (Heidelberg (Germany), 2011).
62. Kirchner H. Ductility and brittleness of bone. International Journal of Fracture. 2006;139:509–516. doi: 10.1007/s10704-006-0050-2. [Cross Ref]
63. Cambridge University. High Performance Computing Service - The Darwin cluster (Date of access: 14/03/2016), (2016).
64. Fritsche A, et al. Measuring bone cell adhesion on implant surfaces using a spinning disc device. Materialwissenschaft und Werkstofftechnik. 2010;41:83–88. doi: 10.1002/mawe.200900547. [Cross Ref]
65. Murphy CM, Haugh MG, O’Brien FJ. The effect of mean pore size on cell attachment, proliferation and migration in collagen-glycosaminoglycan scaffolds for bone tissue engineering. Biomaterials. 2010;31:461–466. doi: 10.1016/j.biomaterials.2009.09.063. [PubMed] [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group