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

Abstract

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.

Introduction

Materials incorporating a fibre structure have been an essential constituent in various fields of application for decades^1,2. Today, a range of several distinguishable categories of fibre materials exists: polymeric non-woven fabrics^3,4, actin networks and cytoskeletons in the field of biomaterials^5,6, or also theoretical computer generated network materials^7,8,9. The mechanical behaviour^9,10,11,12, the thermal conductivity^9,13, or the magnetic response^14,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 fibres^16,17,18,19. 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 elongation¹². The deformation under magnetic actuation of individual fibres has been investigated¹⁵. 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 bioglass^20,21 or cordierite^22,23. In the field of biomedical engineering, metallic^24,25 and since then also non-metallic^26,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 filtering^14,28. This work has been extended also to incorporate bacterial cellulose^29,30. In a different context, magnetic fibre networks could be applied to design heat exchangers of variable drag^31,32. The assembly of this present study has been discussed for its suitability as scaffold in tissue engineering³³. 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 scaffolds^29,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 level^35,36,37. 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 behaviour³⁸. A beneficial effect of a stand-alone magnetic field on bone growth is equally known^39,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 actuation¹⁵. Or global values have been predicted analytically for fibre assemblies^41,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 Cambridge¹⁹.

Mathematical notation

The following notation is used in this document. Scalars are given as x, vectors as $\bar{x}$, and 2^nd rank tensors as $\bar{\bar{x}}$. The vector product “×” and dot product “·” are applied. Relations which are greater-than and approximately equal are indicated as “≳”. All fibre network nodes n of one sample i form the node set N ⁱ. The total number of nodes in that sample i is defined as ${\hat{N}}^{i}$. Equality between two volumes is written “ ⇔ ”, a volume containing a volume subset is given as “⊂”. 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) = ∂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 theory¹². 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)³. In the CT scans, this sample cube side length is the equivalent of 100 pixels.

A skeletonisation algorithm^43,44,45 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]³-grid. Nodes defining network medial axis location are referred to as network nodes in the following. Their total number per sample ${\hat{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 space^9,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.13⁴⁶ (Tables 2 and 3). By their definition in the present study, V _M equals V while V _f is only a subset:

$$V\iff {V}_{M}\subset {V}_{f}$$

(1)

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  = 100³, their volume as V _hex  = (7.75 μm)³.

The values of matrix stiffness E _M are chosen as it would be applicable for three stages in bone tissue development: collagenous bone⁴⁷, granulation tissue⁴⁸, and immature bone⁴⁹. 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 beams^50,51,52 are implemented. Results for linear interpolation (B31) and quadratic interpolation (B32) are analysed and compared to predictions from previous work¹². The beams are simulated for a simplified round cross section of diameter d = 40 μm (i.e. A _f  = (d/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 measurements^12,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 applications⁵³. The regular hexahedra and the beam mesh share the 100³-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 $\bar{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:

$${E}_{f}\gg {E}_{M}$$

(2)

FE simulation and BC

In the case of a void inter-fibre space, the Cauchy stress tensor $\bar{\bar{\sigma }}$ is defined for a fibre network consisting of V _f and the void matrix phase V _void ⁹:

$$\bar{\bar{\sigma }}\ne 0\,\forall \bar{x}\in {V}_{f},\bar{\bar{\sigma }}=0\,\forall \bar{x}\in {V}_{void}$$

(3)

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

$$\bar{\bar{\sigma }}\ne 0\,\forall \bar{x}\in {V}_{f},\bar{\bar{\sigma }}\ne 0\,\forall \bar{x}\in {V}_{M}$$

(4)

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 $\bar{B}$ is imposed (Fig. 1b). The ferromagnetic mechanical response of the fibres is modelled as moment vector ${\bar{\tau }}_{B}$. It is calculated for each beam element of length L _f under the assumption of complete magnetisation parallel to the beam axis¹⁵:

$${\bar{\tau }}_{B}={A}_{f}{L}_{f}(\overline{{M}_{s}}\times \overline{B})$$

(5)

In Eq. 9, ${\overline{\tau }}_{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 444¹⁵. 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 fibres^17,18. The imposed magnetic induction vector $\bar{B}=[\begin{array}{c}1\\ 0\\ 0\end{array}]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 magnitude^12,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(\bar{\bar{\sigma }})=0$ is fulfilled. Body force per volume $\overline{f}$ is set to 0 in the present study; gravitational forces being a typical example for $\overline{f}$. Together with the surface traction vector $\bar{t}$, it defines the equilibrium of forces. Surface traction itself is defined by $\bar{\bar{\sigma }}$ and the outer-normal vector: $\bar{t}=\bar{\bar{\sigma }}\,\cdot \,{\bar{n}}_{out}$ ⁵⁵.

$${\int }_{S}\bar{t}\,dS+{\int }_{V}\overline{f}\,dV=0$$

(6)

As the magnetic response is modelled as moment vector ${\overline{\tau }}_{B}$, the general condition for equilibrium of forces can be simplified for the present study. Considering that $\overline{f}\mathrm{=}\,0$ it follows that:

$${\rm{Eq}}{\rm{.}}\,{\rm{6}}\mathop{\Rightarrow }\limits^{\bar{f}\mathrm{=}\,0}{\int }_{{F}_{x}}\bar{t}\,dS=0$$

(7)

Relative to the origin of point vector $\bar{x}$, the general form of the equilibrium of moments⁵⁵ is defined:

$${\int }_{S}(\bar{x}\times \bar{t})\,dS+{\int }_{V}(\bar{x}\times \bar{f})\,dV=0$$

(8)

The simplification of Eq. 6 by $\bar{f}=0$ can be repeated for Eq. 8. The sum of ${\bar{\tau }}_{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 ${\bar{\tau }}_{B}$ follows as:

$${\rm{Eq}}\,{\rm{.8}}\mathop{\Rightarrow }\limits^{\bar{f}=\,\mathrm{0,}\,+{\bar{\tau }}_{B}}\quad {\int }_{{F}_{x}}(\bar{x}\times \bar{t})\,dS\mathop{\underbrace{+\sum _{{V}_{fibre}}{\overline{\tau }}_{B}}}\limits_{{\rm{magnetic}}\,\,{\rm{actuation}}}=0$$

(9)

The Lagrangian reference frame⁵⁶ is used for Eqs 6 to 9. The alternative Eulerian reference frame can be advantageous for the modelling of fluids^57,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 (>10³) 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=[\begin{array}{c}0\\ L/2\\ L/2\end{array}]=[\begin{array}{c}0\\ 387.50\\ 387.50\end{array}]$ μ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}=\sum _{V\cap {o}_{(CP,R)}}{N}^{i}/{\hat{N}}^{i}$$

(10)

$$\begin{array}{ccc}{Y}_{(R)}^{i} & = & {\sum }_{V\cap {O}_{(CP,R)}}{N}^{i}/{\hat{N}}^{i}\\ & = & {\int }_{0}^{R}{y}_{(R)}\,dR\\ i & \in & \{10{\rm{ \% }},15{\rm{ \% }},20{\rm{ \% }}\}:\text{samples of 10, 15, and 20\%}\,f\end{array}$$

(11)

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 materials⁶¹. The material behaviour type of bone has been shown to depend on the deformation rate dε/dt ⁶².

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 ⁱ (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.

$${g}_{(R)}^{i,j}=\frac{1}{S}{\int }_{V\cap {o}_{(CP,R)}}{\varepsilon }_{j}\,dS$$

(12)

$$\begin{array}{cc}{f}_{(R)}^{i,j} & =\mathop{\underbrace{\frac{1}{{\hat{N}}^{i}}\quad \sum _{{N}^{i}}}}\limits_{2.\,{\rm{a}}{\rm{v}}{\rm{e}}{\rm{r}}{\rm{a}}{\rm{g}}{\rm{e}}\,{\rm{f}}{\rm{o}}{\rm{r}}\,{\rm{n}}{\rm{o}}{\rm{d}}{\rm{e}}\,{\rm{s}}{\rm{e}}{\rm{t}}\,{N}^{i}}\mathop{\overbrace{[\frac{1}{\frac{4}{3}\pi {R}^{3}}\,\,{\int }_{{O}_{(n,R)}}{\varepsilon }_{j}dV]}}\limits^{1.\,{\rm{a}}{\rm{v}}{\rm{e}}{\rm{r}}{\rm{a}}{\rm{g}}{\rm{e}}\,{\varepsilon }_{j}\,{\rm{i}}{\rm{n}}\,O\,{\rm{o}}{\rm{n}}\,{\rm{n}}{\rm{o}}{\rm{d}}{\rm{e}}\,n}\\ i\, & \in \{10{\rm{ \% }},15{\rm{ \% }},20{\rm{ \% }}\}:\text{samples of 10, 15, and 20\%}\,f\\ {\varepsilon }_{j}\, & \in \{{\varepsilon }_{v.M.},\,{\varepsilon }_{Pmax}\}\end{array}$$

(13)

Hardware

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)⁶³.

Results

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 seen¹⁷. 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 25³-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 5b).

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 stimulation³⁸. 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 literature⁵⁹ 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⁻⁴ 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 effect¹². 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 stiffness¹²:

$${E}_{{\rm{B}}31}\gtrsim {E}_{{\rm{B}}32}$$

(14)

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.

$$\begin{array}{cc}\text{if}\,\,B=const,\,\,\, & {E}_{M}=const,\,\,\,f=const\text{, it holds that}\\ {\varepsilon }_{Pmax} & \,\lesssim \,{\varepsilon }_{v.M.}\,\,\,\text{, valid for both Timoshenko beams (B31 or B32)}\\ {\varepsilon }_{\text{B31}} & \,\lesssim \,{\varepsilon }_{\text{B32}}\,\,\,\text{, valid for both equivalent strain measures (}{\varepsilon }_{v.M.}\,\text{or}\,{\varepsilon }_{Pmax}\text{)}\end{array}$$

(15)

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 5²-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 deformation¹². 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 stimulation³⁸ 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 ∈ [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 forces⁶⁴. 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 10⁶ 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:

$$\begin{array}{c}{\varepsilon }_{(R)}=aB{e}^{bR}\\ a\,[{{\rm{T}}}^{-1}],\,b\,[{{\rm{\mu }}{\rm{m}}}^{-1}]\in {{\mathbb{R}}}^{+}\end{array}$$

(16)

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 10⁶ 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 ². The values are given for ε _v.M., ε _Pmax , and for both Timoshenko beams. r ² 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.

Discussion

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 deformation¹².

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 material¹².

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 ² = 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 growth³⁸. 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 structure⁶⁴. Also a dynamic matrix model including the active migration of bone cells through the fibre network⁶⁵ 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 pseudarthrosis³⁹. Designs for hip joint prostheses make use of that effect⁴⁰. 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.