Multiscale modeling of cancellous bone considering full coupling of mechanical, electric and magnetic effects.

Modeling of cancellous bone has important applications in the detection and treatment of fatigue fractures and diseases like osteoporosis. In this paper, we present a fully coupled multiscale approach considering mechanical, electric and magnetic effects by using the multiscale finite element method and a two-phase material model on the microscale. We show numerical results for both scales, including calculations for a femur bone, comparing a healthy bone to ones affected by different stages of osteoporosis. Here, the magnetic field strength resulting from a small mechanical impact decreases drastically for later stages of the disease, confirming experimental research.

Introduction

In the present contribution, we develop a multiscale model for cancellous bone taking mechanical, electric and magnetic effects into account. An important application of this model is the early detection of osteoporosis. This bone disease reduces the mass density of the bone, making it thinner and weaker, increasing the likelihood of fractures. Sonography is used as a cheap, fast and non-invasive early detection technique for osteoporosis (Kaufman et al. 2008). Material modeling and numerical simulations are helpful tools in order to understand and evaluate experimental measurements and enable medical diagnostics based on this method.

Bone is a composite material with impressive properties, drawing the interest of researchers of many different fields. As a material, it is very strong and stiff and has a high fracture toughness, while also maintaining a light weight (Hamed et al. 2010). Thus in recent decades, a lot of different approaches to investigate and simulate the material behavior of bone have appeared. Many analytical solutions are based on Biot’s famous theory (Biot 1956a, b). Examples include (Buchanan and Gilbert 2007; Chen et al. 2018; Steeb 2010). Here, cortical bone is modeled as a solid, while bone marrow is assumed to be a fluid. The acoustic properties of bone material are then used to obtain mechanical material parameters of bone and the parameters of Biot’s model. Additionally, the results are compared with the findings of experiments.

In contrast to the analytical solutions, many numerical approaches exist in the scope of bone modeling as well. The finite difference method was used in Kaufman et al. (2008) to obtain numerical results of ultrasound propagation in bone. Applications of the finite element method (FEM) on the topic of bone modeling include the simulation of mechanical properties of bone (Gardner et al. 2000; Miller et al. 2002) and the simulation of osteogenic effects (Wang et al. 2017). In Christen et al. (2010), patient-specific FEM simulations are proposed in order to estimate the likelihood of osteoporotic fractures.

Since the bone microstructure is very complex and heterogenous, material modeling should take place on different scales. Currently used single-scale models are criticized in Christen et al. (2010) as oversimplified and multiscale approaches are proposed instead. In Hamed et al. (2010), the mechanical properties of bone are modeled on five different length scales from the nanoscale to the macroscale. Multiscale approaches can also be combined with numerical methods. The finite element square method () extends the standard FEM approach by applying the multiscale concept and solving the differential equation systems on two scales via the FEM. An overview of the method can be found in Schröder (2000); Schröder and Hackl (2013). Basic works on this method include for example (Willis 1981; Suquet 1987; Castañeda and Suquet 1997) and applications to different materials can be found for example in Ilić and Hackl (2004); Miehe et al. (2002). An application of the within the scope of bone modeling can be found in Ural and Mischinski (2013); Podshivalov et al. (2011); Pahr and Zysset (2008), proposing different models to capture the microstructure of bone, allowing to investigate mechanical effects. In Ilic et al. (2010) and Klinge et al. (2013), macroscopic material parameters were recovered by simulations on the microscale. The results obtained were subsequently used for macroscale simulations of wave propagation.

So far, all presented contributions focus only on the mechanical effects of bone. However, cortical bone possesses the properties of a piezoelectric solid. After the discovery of this effect (Fukada and Yasuda 1957; Shamos et al. 1963), research considering these coupled physical effects has started. A review on computer modeling of bone piezoelectricity can be found in Mohammadkhah et al. (2019). There, applications are discussed as well. Since electric and magnetic effects are coupled physically via the Maxwell equations, it may be necessary to include magnetic effects as well. In Güzelsu and Saha (1981), bone was modeled as a hollow cylinder and analytical solutions of the coupled equations of all three effects were studied. The results were then compared to in vitro experimental measurements.

In this work, we present a fully coupled multiscale approach for modeling cancellous bone considering mechanical, electric and magnetic effects and using two scales, the macro- and microscale. At the microscale, we assume a heterogenous material consisting of two phases, cortical bone and bone marrow. Cortical bone is modeled as piezoelectric, insulating solid, bone marrow as viscoelastic, conducting solid. Electric and magnetic effects are coupled via the Maxwell equations. Based on energy methods in mechanics, we establish a thermodynamically consistent material model and derive the weak and strong form of the corresponding boundary value problem. We apply the FEM to solve the problem numerically. For multiscale analysis, we resort to the method. To apply this method, we constructed a periodic representative volume element (RVE) and discuss the transition between scales.

The article is structured as follows: in Sect. 2 we discuss the material structure of cancellous bone and the method. Then, we introduce the microscopic material model and derive the weak and strong form of the corresponding variational problem. Additionally, we cover the macroscale boundary value problem. In Sect. 3 we present the FEM implementation of the model and show details regarding scale transition and programming. In Sect. 4 we present numerical results, starting with microscale calculations, on to multiscale simulations for a cylindrical body and finally a true to scale model of a human femur bone. To close this article, we draw a short conclusion and give an outlook to future research envisioned in Sect. 5.

Material model

Structure and properties of cancellous bone

Our work focuses on the description of the internal structure of cancellous (spongy) bone, which consists of small beams or shells of interconnected cortical bone and interstitial bone marrow. Cortical bone is mainly composed of elastic collagen fibers, which act as charge carriers. When applying a shear stress, these collagen fibers slip past each other, thus producing a piezoelectric effect. This was first measured in Fukada and Yasuda (1957) and later validated in Shamos et al. (1963). This means that, whenever a mechanical strain is present in the bone, an electric field is generated due to the piezoelectric effect. A time-dependent fluctuation of the electric field then creates a magnetic field due to Ampère’s circuital law, coupling mechanical, electric and magnetic effects all together.

An important application of bone modeling is the early detection of osteoporosis, a bone disease, which manifests itself in the reduction of the cortical bone phase, thus reducing the strength of the bone and increasing the likelihood of fractures. Compared to a healthy bone, the volume fraction of cortical bone for an affected bone can be reduced from 30 to (Steeb 2010; Ilic et al. 2010). Figure 1 shows a comparison depending on the osteoporosis stage and illustrates the heterogeneity of the material. During the course of osteoporosis, the cortical bone (represented brighter) reduces and is replaced by bone marrow (represented in dark). Thus, we will employ different RVEs for the simulations. Here, the cortical bone phase is represented in gray, while the bone marrow phase is drawn in transparent red color.

Fig. 1

Bone phases depending on osteoporosis stage (cf. Laboratoires 2019) and corresponding RVEs

Early detection of osteoporosis can be done via sonography: ultrasonic waves enter the bone and due to the described effects create a magnetic field, which can be measured (Güzelsu and Saha 1981) and—depending on the results—conclusions on the health status of the investigated bone can be drawn. In this contribution, we introduce a material model including all the described effects. It is important to note, that there are two different forms of coupling: while the piezoelectric coupling is captured via a suitable material model, the Maxwell coupling is of physical (electrodynamical) nature.

Concept of the method

To include micro-heterogeneities directly, an extremely fine resolution of the problem would be necessary, resulting in a very high computation cost for the simulations. Alternatively, the method is a homogenization technique, which captures the structure of micro-heterogeneities by introducing a second—smaller—scale to the problem. If the material is statistically regular on the smaller scale, it can be modeled by a corresponding RVE (Schröder 2000; Schröder and Hackl 2013). In this paper, we denote the larger scale as the macroscale and the smaller scale as the microscale. To obtain accurate results, the quotient of the characteristic lengths between micro- and macroscale should tend to zero, so the RVE has to be much smaller than the simulated macroscopic body. Figure 2 illustrates this procedure: instead of using a material model on the macroscale, the state variables are linked to the microscale, where the RVE problem is solved. The microscale calculations yield average flux quantities and consistent tangent matrices, which then can be used for the solution of the macroscale problem, replacing a macroscopic material model.

Fig. 2

Transition between macro- and microscale. State variables enter as boundary conditions of the RVE problem. Flux quantities at the macroscale are calculated by averaging the RVE quantities

We denote spatial coordinates on the macroscale by and on the microscale by . Quantities denoted as are affiliated to the macroscale. The transition between the scales regarding energy conservation and numerical treatment is discussed in Sect. 3.2.

Variational formulation of the microscale problem

The domain , representing the RVE of the micro-problem, is split into a cortical bone part and a bone marrow part . For any quantity, the indices and are used to denote the affiliation to each phase. If no index is present, the quantity or equation is valid for both phases. We employ the following thermodynamic energy functional at the microscale:The functional contains the energy densities and of both phases, a volume constraint , dissipation and gauge functionals ( and ) and the potential of the generalized external forces . The main variables of the problem are then the mechanical displacements , the electric scalar potential and the magnetic vector potential , yielding seven unknown variables for the three-dimensional model. The state variables are the mechanical strain , the electric field and the magnetic flux density , calculated asThis way, two of the four Maxwell equations are already satisfied:For the mechanical strain, we use Voigt’s notation (Mehrabadi and Cowin 1990) asThen, the energy densities for both phases areconsisting of quadratic energies for mechanical, electric and magnetic effects, resulting in a linear problem. We include a piezoelectric energy term for the cortical bone phase. For the bone marrow phase, an inelastic strain is introduced. Here, is the mechanical stiffness tensor, is the permittivity tensor, is the inverse permeability tensor and is the piezoelectric tensor. While it is possible to switch between state and flux variables via a Legendre transformation, the present formulation proves as the most suitable for our model, as it allows an easy inclusion of the Maxwell coupling and the electric dissipation. For linear problems, the transformation would change an extremal into a saddle point problem, thus excluding solvers, that require positive definiteness of the system matrix as a precondition. The constraint function readsenforcing volume conservation of the inelastic deformation. Here, is a Lagrange multiplier. The dissipation function isThus, governs the evolution of the inelastic strain and the energy loss due to conduction. The latter satisfies Ohm’s law (Eq. (7), right). Both parts of the dissipation only occur in the bone marrow phase. Here, the viscosity parameter , the electric conductivity tensor , with the identity tensor , the electric conductivity and the electric current density are introduced. The gauge function isand ensures, that a unique solution for the magnetic vector potential is obtained by penalizing its divergence, effectively requiring, that vanishes and thus improving the numerical stability (Semenov et al. 2006). The penalty parameter is a numerical parameter used to control the gauge term. Finally, the potential of generalized external forces isHere, and are the mechanical volume and surface forces, and are the electric volume and surface charges and and are the volume and surface currents.

By calculating the derivative of the energy density with respect to the state variables, we find the following constitutive equations for both phases:For the bone marrow, the additional constitutive equations areintroducing the flux quantities mechanical stress , electric displacement and magnetic field strength . For the cortical bone phase, the viscosity parameter and the electric conductivity tensor vanish. The material tensors satisfy

Weak and strong form of the microscale problem

To calculate the weak and strong form of the problem, the energy functional has to become stationary with respect to the main variables and internal variables, leading toThe stationary condition of the first variation of the energy functional reads thenThe variation of the generalized external forces isUsing the introduced energy densities, constraint, dissipation and gauge functions, Eqs. (5), (6), (7), (8), and inserting the constitutive equations Eqs. (10), (14) simplifies toHere, the identity vector is denoted as . We find the evolution equation of the inelastic strain:To calculate the Lagrange multiplier, the trace is applied to Eq. (17):The second term in Eq. (18) must vanish because of the introduced volume constraint. This leads to the final evolution equationwith denoting the deviatoric part of the mechanical stress . The time integration of the evolution equation is discussed in Sect. 3.

To calculate the strong form of the problem, the remaining variational equation is used:This form is later used to insert a FEM ansatz. We apply partial integration to each term. Details can be found in Appendix A. We obtainrecovering the mechanical equilibrium condition, the two remaining Maxwell equations and boundary conditions, including the gauge. Here, is the normal vector pointing outwards. Additionally, we receive the jump conditions between the phaseson the interface and the evolution equation of the inelastic strain Eq. (19) in . Here denotes the difference between the phases. It should be noted that the strong form is valid for both phases, but the calculation of the flux variables and the inelastic strain evolution depends on the specific material parameters and thus in which phase the calculation is done.

Macroscale problem

For the macroscale, the following boundary value problem in the domain has to be solved: find the set , such thatwith the state variablesand the calculation of the fluxes depending on the microscale calculationsWe transform the strong form into the weak form by multiplying with test functions of the main variables and again using partial integration:Here, the variation of the macroscopic generalized external forces isThis form is again used in the next section to formulate the FEM.

Numerical implementation

Finite element method

To solve the boundary value problems on both scales, we insert a standard finite element approach (Zienkiewicz et al. 2005) into the weak form of the problem for all main variables. In this section, we derive the resulting system for the microscale. It should be noted that the same system has to be solved for the macroscale, but each quantity has to be replaced by its macro-average quantity . The inelastic strain is only present on the microscale and vanishes on the macroscale. Its calculation is not done via the FEM, but directly by using the evolution equation Eq. (19) on the integration point level. Details regarding the calculation of macro-fluxes and consistent material tensors are given in the next subsection. Here, we denote nodal FEM values by . For the evolution equation of the inelastic strain on the micro-scale, we apply an explicit Euler scheme, yielding:Here, is the time increment between two time steps. The standard FEM approach for the remaining system isapproximating the main variable and their variations by shape functions times the nodal values of the functions . For the state variables and the gauge, this approach yieldsHere, the operator matrices areInserting these equations into the reduced weak form of the micro-problem Eq. (20) and by using the arbitrariness of the test functions, we find the final equation system in matrix form as follows (a detailed derivation is given in Appendix B):with the residual vector and the generalized force and displacement vectors together with the mass, damping and stiffness matrices as follows:The material tensors depend again on the phase. We calculate the mechanical stiffness tangent matrix by introducing a time discretization as:For the bone marrow phase, the calculation depends on the inelastic strain :withThe mechanical stiffness tangent matrix for the bone marrow phase is thenIn order to solve the resulting second-order differential equation system, a suitable time integration scheme is necessary. Here we use a JWH--scheme introduced in Kadapa et al. (2017), where also details regarding advantages and implementation of this method can be found. For the time integration, the time increment and the additional numerical parameter are needed. By combining the method with a regular Newton–Raphson scheme, we transform the matrix system of Eq. (32) towith the index denoting the iteration and the generalized tangent matrixwhich is the Jacobian of the system. Here is the increment of the solution vector and , and are numerical parameters depending on (Kadapa et al. 2017). The residual is calculated from either initial conditions for the first iteration of the first time step or else from the previous increment (Kadapa et al. 2017). The resulting tangent matrix is neither symmetric nor positive definite, limiting the choices for a suitable solver of the linear system.

Transition between the scales

To connect the macro- and microscale in , it is important to discuss the transition between the scales. The Hill-Mandel conditions (Hill 1963, 1972; Schröder 2009; Schröder et al. 2016; Labusch et al. 2019; Karimi et al. 2019) have to be fulfilled, guaranteeing energy conservation during the scale transition. Thus, the virtual work on the macroscale has to be equal to the virtual work on the microscale:For the macro-to-micro-transition, these conditions can be fulfilled by three different types of boundary conditions on the microscale: Dirichlet, Neumann and periodic boundary conditions (Ilic et al. 2010; Schröder 2000; Schröder and Hackl 2013). Here we chose periodic boundary conditions, as they are the only type of boundary condition, where the results on the microscale are independent from the relative geometry of the RVE (Schröder 2000; Schröder and Hackl 2013). Additionally, as the RVE is periodic in space, this type of boundary condition is the most suitable. In the program, the periodic boundary conditions were applied by fixing all degrees of freedom at all corner nodes, preventing rigid body motions, and linking all degrees of freedom at opposite faces of the RVE, ensuring the periodicity. The micro-state variables consist then of two parts: a term resulting from the microscopic main variables (denoted by ), whose fluctuations are calculated, and a term contributed by the macroscale:This way, we calculate the flux variables on the microscale. For the micro-to-macro transition, the volume average of these flux quantities is sent back to the macroscale:In this model, energy dissipation is considered in two ways. For the electric current , the average is calculated and included in the scale transition, resulting in no energy loss during the scale transition. For the inelastic strain , the complete state in every point and for every RVE is saved. Thus, the dissipation occurs only on the microscale and the energy conservation is fullfilled, as the virtual work send to the microscale is equal to the virtual work send back added to the energy dissipation on the microscale. With the flux variables available on the macroscale, it is now possible to obtain the macro-residual for the Newton–Raphson method and the calculation of consistent macro-tangent moduli remains, which are needed for the iteration. The definitions of those moduli readThe calculation can be done by applying a small numerical perturbation to each entry of the corresponding state variablewith the i-th unit vector , and then calculating each entry of the macroscopic tangent tensors by evaluating the perturbated fluxes by means of the RVE asSince for our model the same RVE is used everywhere and the nonlinearity from the inelastic strain is very small, this calculation has to be done only once for all RVEs and all time steps, making this approach very efficient. Together with the calculated macrostate variables, this allows to solve the macroscopic FE problem.

Implementation

For the simulations, we implemented a computer program in the language julia (The julia programming language 2021), using mainly the packages juafem (Carlsson and Ekre 2021) and coherentstructures (de Diego et al. 2021). As the microscale calculations are not dependent on each other, we have parallelized the macroscale element routine, increasing the speed of the computations drastically. As the inelastic strain is only present in the microscale, we used HDF5 files to store the complete state of the inelastic strain for every RVE for the previous and current time step. Thus, for the inelastic strain evolution no information is lost. In order to solve the linear systems, we used the BiCGStab(l) method of the package krylovmethods (Krylovmethods 2021), as it is stable, fast even without preconditioning the problem and can be used for any matrix type. Regarding the structure of the program, Fig. 3 shows the procedure.

Fig. 3

Program flow of the multiscale simulations

Simulation results

Parameters and material tensors

In this subsection, we discuss the numerical and material parameters employed. Unless explicitly stated otherwise, the parameters from this subsection are used in all simulations. Regarding the numerical parameters, we use the same parameters for both scales. The time integration parameter is , the Newton–Raphson tolerance is and the gauge penalty parameter is . The load and numerical time step increment depending on the model are shown in Table 1.

Table 1

Load and numerical time step increment for the different models

Model	Load	Time step
Microscale cube	ε¯_yz=1·10^-5	Δ_t=1·10^-3 s
Cylinder	u_max=2·10^-6 m	Δ_t=1·10^-2 s
Femur bone	u_max=2·10^-6 m	Δ_t=1·10^-2 s

The used default material parameters are shown in Table 2. Young’s modulus and Poisson’s ratio for both phases can be found in Steeb (2010). The piezoelectric coefficient can be found in Fukada and Yasuda (1957). For the magnetic properties, bone is considered as a nonmagnetizable material, thus having the same permeability as the vacuum (Güzelsu and Saha 1981). All other parameters are of rather academical nature and influence the results only marginally. The resulting material tensors read

Table 2

Default material parameters

Material parameter		Cortical bone		Bone marrow
Young’s modulus	E	22.0	GPa	2.0	GPa
Poisson’s ratio	ν	0.32	–	0.3	–
Permittivity	ξ_1	8.85·10^-12	F/m	8.85·10^-12	F/m
Permeability	μ_c	1.257·10^-6	H/m	1.257·10^-6	H/m
Piezoelectric coefficient	e_15	3.0·10^-3	As/m^2	0	As/m^2
Electric conductivity	κ_1	0	S/m	1.0·10^4	S/m
Viscosity parameter	μ_v	0	s/GPa	0.5·Δ_t	s/GPa

We assume linear isotropic material everywhere, excluding the piezoelectric tensor which is preferential in the z-axis due to the longitudinal orientation of the collagen fibers. It should be noted, that due to the form of the piezoelectric tensor, the material model as a whole is non-isotropic.

For the generation of the meshes, we used the program gmsh (Geuzaine and Remacle 2021). We did the visualization of the results with paraview (Paraview 2021) and julia (The julia programming language 2021).

Microscale model

In this subsection, we restrict ourselves to microscale simulations. In order to compare periodic RVEs for different stages of osteoporosis, we introduce the lengths parameters a and b (Fig. 4), which allow us to control the volume fractions of the phases. By using this convention, the total volume of the RVE is . We only use RVEs with the same total volume of , which is a suitable size for the microscale calculations (Ilic et al. 2010), making it easy to compare different RVEs. Thus, the choice of a and b is restricted by . The volume fraction of cortical bone for our RVE is .

Fig. 4

Periodic RVE with cortical bone phase (gray) and bone marrow phase (transparent red) and lengths parameters

In our first example, we use a healthy bone RVE with the parameters and , resulting in . We compare different mesh resolutions. The first RVE consists of two elements in each phase block, resulting in six elements for each spatial direction. The second RVE consists of four elements in each block, resulting in twelve elements for each spatial direction. Here, all degrees of freedom for all corner nodes are restricted to zero and all opposite nodes are linked, to guarantee periodicity. Figure 5 shows the results of the simulations.

Fig. 5

Microscale simulation results of a coarse and fine mesh (left and right respectively) for all flux quantities. Top left: mechanical stress , top right: mechanical stress in the xz-plane with , bottom left: magnitude of the electric displacement field , bottom right: magnitude of the magnetic field strength

Both simulations show quadratic convergence behavior and periodic results. For all quantities, the results between the two different used meshes look nearly identical confirming mesh independence of the results. This is not only fulfilled on the surface of the model, but also in the inner parts, as the slice (top right) shows. It should be noted that since the method uses volume averaging, the coarse mesh with only six elements in each spatial direction is sufficient enough to create accurate results for the multiscale method and is mostly used in the remaining examples of this paper.

The calculation of the magnetic field strength is susceptible for numerical errors. These errors can occur, when there are sharp edges in the mesh or between the phase transitions, which can amplify the results. To investigate this issue, we constructed smoother RVEs by using tetrahedron elements and different mesh resolutions. Figure 6 shows the used meshes. We found that despite the smoother approach, the numerical results of the RVEs with tetrahedron elements are much worse compared to the RVEs with hexahedron elements, showing worse convergence behavior and overestimating the magnetic field strength. One reason for this result could be that through the mesh refinement, many additional corners are introduced, which in total amplify the magnetic field strength more than a low number of very sharp corners. Moreover, smaller element size leads to amplified singularities of the corresponding fields at corners, which otherwise are regularized by the employed shape functions. Similarly, Fig. 5 shows an increase in the magnetic field strength for the finer mesh resolution of the RVEs with hexahedron elements. In conclusion, the coarse RVE with hexahedron elements shows the best numerical performance despite the low mesh resolution and the included sharp edges.

Fig. 6

Used meshes for the complete RVE (left) and only the cortical bone phase (right). a coarse hexahedron mesh, b fine hexahedron mesh, c coarse tetrahedron mesh, d fine tetrahedron mesh

To compare the model behavior for different stages of osteoporosis, we created RVEs with different volume fractions of cortical bone. Table 3 shows the choice of the lengths parameters and the resulting volume fractions. The macroscopic mechanical stiffness tensor was now evaluated for all RVEs by applying a small numerical perturbation as discussed in Sect. 3.2. We then calculate the effective Young’s modulus asFigure 7 shows a plot of the macroscopic Young’s modulus against the volume fraction of cortical bone. Here, we observe a drastical reduction of the macroscopic Young’s modulus with decreasing cortical bone fraction. Compared to a healthy bone (), the effective Young’s modulus of the degenerated bone () decreases to (from to ). Similar results can be found in Ilic et al. (2010).

Table 3

Lengths parameter of the different RVEs

no.	a [mm]	b [mm]	ρ_b
1	0.43	0.14	5.3%
2	0.40	0.20	10.4%
3	0.38	0.24	14.5%
4	0.36	0.28	19.1%
5	0.34	0.32	24.2%
6	0.32	0.36	29.5%

Fig. 7

Effective Young’s modulus against cortical bone volume fraction for different RVEs

Cylinder model

In this section, we show results for a cylinder model, which has a length of and a diameter of . The mesh and the displacement boundary conditions are shown in Fig. 8. The mesh consists of 1767 nodes and 1440 hexahedral elements. The left and right face is fixed, resulting in the boundary conditions on the faces. Additionally, in the inner part of the left face () depicted in Fig. 9, the cylinder is assumed to be grounded, resulting in and . We apply a time-dependent mechanical displacement in x-direction to the middle part of the cylinder and calculate 100 time steps. Figure 10 shows the amplitude of the displacement function a versus the time t.

Fig. 8

Cylinder mesh and displacement boundary conditions (red: all directions restricted, orange: only the x-direction restricted, blue-gray: no directions restricted)

Fig. 9

Cylinder front in the xy-plane for with grounded nodes in red

Fig. 10

Amplitude of the displacement function a against the time step t

First, we examine the simulation results for the healthy bone (RVE 6, ). Here, we observe quadratic convergence behavior for the macroscale as well. Figures 11 and 12 show the magnitude of the average electric displacement field and the magnitude of the average magnetic field strength , respectively, plotted against time t. The history of the average electric displacement field mimics the displacement boundary condition. Thus, the electric displacement field is caused mainly by the piezoelectric effect of the cortical bone material phase. In contrast, the magnitude of the average magnetic field strength increases until time , where the maximum is reached. Then, the magnitude decreases again and at the end of the simulation, only a small amount of the magnetic field is present. We conclude, that the magnetic field is caused mainly by the time change of the electric displacement field as described by the Maxwell equations.

Fig. 11

Magnitude of the average electric displacement field , plotted against the time t

Fig. 12

Magnitude of the average magnetic field strength , plotted against the time t

To compare the different stages of osteoporosis, we use different RVEs (Table 3). The simulation results are shown in Figs. 13, 14, 15, 16. Here, the number of the specific RVE increases from top to bottom.

Fig. 13

Simulation results for RVE 1 (top) to 6 (bottom): stress ,

Fig. 14

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the electric displacement field ,

Fig. 15

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the magnetic field strength ,

Fig. 16

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the electric current density ,

As an additional example for the cylinder model, we performed a parameter study for the electric conductivity parameter , aiming to understand the interaction between the time derivative of the electric displacement field and the electric current density in the Maxwell equation. Figures 17 and 18 show the results for RVE 1 and .

Fig. 17

Simulation results for RVE 1 for the magnetic field strength with (top), (in the middle) and (bottom),

Fig. 18

Simulation results for RVE 1 for the electric current density with (top), (in the middle) and (bottom),

For all quantities, we observe an increase for RVEs with higher volume fractions of cortical bone. Additionally, the difference between the RVEs is greater, the lower the volume fraction of cortical bone is. While the difference is barely noticeable between RVE 5 and 6, the change of all quantities excluding the stress is distinct between RVEs 1 and 2. Qualitatively, we notice similar results between the different RVEs.

Regarding the parameter study of the electric conductivity, we observe nearly identical results for the magnetic field strength for the first two choices of , but a significant increase for . Similarly the electric current density increases proportionally to the increase of the material parameter. Thus, for the first two choices of , nearly no magnetic field and electric current is visible.

Real bones can be highly anisotropic. To investigate possible effects on the simulation results, we constructed an anisotropic RVE, which is longer (Fig. 19) and therefore also is divided into ten instead of six elements in z-direction. We used the parameters and , resulting in a total RVE volume and a volume fraction of cortical bone . We compare our calculations to the isotropic RVE 6, which has a similar volume fraction of cortical bone. The results are shown in Figs. 20 and 21. We obtain similar results for both RVE geometries. The calculated stress is slightly higher for the anisotropic RVE. The magnetic field strength is about increased for the anisotropic RVE compared to the cubic RVE.

Fig. 19

Anisotropic RVE with cortical bone phase (gray) and bone marrow phase (transparent red) and lengths parameters

Fig. 20

Simulation results for RVE 6 (top) and the anisotropic RVE (bottom) for the stress ,

Fig. 21

Simulation results for RVE 6 (top) and the anisotropic RVE (bottom) for the magnetic field strength ,

True to scale bone model

We examine a true to scale model of a human femur bone from Lifescience (2021) and slightly modify it by using the software blender (Blender 2021), improving the mesh. Again the model has a length of about . The mesh and the displacement boundary conditions are shown in Fig. 22. The mesh consists of 1660 nodes and 4944 tetrahedral elements. The grounded nodes are shown in Fig. 23. Again, we apply the mechanical displacement depicted in Fig. 10 to the middle section and calculate 100 time steps. To compare different stages of osteoporosis, we use again different RVEs (Table 3). Figures 24, 25, 26, 27, 28, 29 show the results.

Fig. 22

Femur bone mesh and displacement boundary conditions (red: all directions restricted, orange: only the x-direction restricted, blue-gray: no directions restricted)

Fig. 23

Femur bone front with grounded nodes in red

Fig. 24

Simulation results for RVE 1 (top) to 6 (bottom): stress ,

Fig. 25

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the electric displacement field ,

Fig. 26

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the magnetic field strength ,

Fig. 27

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the electric current density ,

Fig. 28

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the magnetic field strength , slice,

Fig. 29

Simulation results for RVE 1 (top) to 6 (bottom): magnitude of the electric current density , slice,

Again, the simulations show qualitatively similar results, but a significant increase for all quantities the higher the cortical bone volume fraction is. Compared to the cylinder model, we receive slightly higher numerical values, which lie in the same magnitudes. The reason for this is most likely the used mesh, which has sharper corners due to the geometry of bone. Additionally, tetrahedron elements usually perform worse compared to hexahedron elements. The difference between the RVEs is smaller the higher the volume fraction of cortical bone is. Thus, both the functionality of the bone and the results of the sonography are only slightly affected at earlier stages of osteoporosis, but significantly at later ones. This confirms the disease as being often imperceptible for many subjects at earlier stages. This is especially important regarding the magnetic field strength , as it is the quantity measured at sonography-aided early detection. To further examine the results, we calculate the average and maximum magnetic field strength at time step for the different RVEs. The results are shown in Figs. 30 and 31.

Fig. 30

Average magnetic field strength for the different RVEs at

Fig. 31

Maximum magnetic field strength for the different RVEs at

Here, for both quantities a similar behavior can be observed. While there is nearly no reduction between the two RVEs with the highest volume fraction of cortical bone, the difference between the single RVEs increases for lower volume fractions of cortical bone. The average magnetic field strength reduces for the ill bone () to compared to the healthy bone (), from to . The maximum magnetic field strength for the healthy bone is , while the maximum for the degenerated bone is only . This equals a reduction to . These results show the order of magnitude to be expected for the results of experimental research. For advanced stages of osteoporosis, sonography should measure a magnetic field strength, whose magnitude is only about one third compared to a healthy bone.

Conclusion and outlook

In this contribution, we present a fully coupled multiscale model for cancellous bone considering mechanical, electric and magnetic effects. We model bone as a two-phase material with the cortical bone phase assumed as a piezoelectric, insulating solid and the bone marrow phase described as a viscoelastic, conducting solid. Electric and magnetic effects are coupled via the Maxwell equations. Based on energy methods in mechanics, we establish a thermodynamically consistent material model and derive the weak and strong form of the microscale boundary value problem.

In order to solve the macroscale problem, we create an RVE and apply the FEM to solve the problem numerically. For the time integration of the FEM, we use a JWH--scheme (Kadapa et al. 2017). The numerical simulations on the microscale show mesh independence and quadratic convergence. For finer mesh resolutions or smoother geometries of the phases, the model tends to overestimate the magnetic field strength. Additionally, we show that the effective Young’s modulus of the RVE depends strongly on the volume fraction of the different phases. Here, we find a reduction by for the degenerated bone () compared to the healthy bone (), achieving similar results as in Ilic et al. (2010).

For the multiscale calculations, we use and apply periodic boundary conditions and volume averaging for the transition between the scales. We apply a time-dependent displacement boundary condition. The macroscopic cylinder model again shows quadratic convergence. To compare different stages of osteoporosis with a healthy bone, we create six different RVEs with different volume fractions of cortical bone phase and run calculations for all RVEs. The simulations show a strong reduction of all quantities with decreasing volume fraction of cortical bone phase. The differences between the healthy bone RVE () and a slightly degenerated bone () are very small, while the differences in the later stages of the illness, ( compared to ), increase drastically. To examine the interaction between the time derivative of the electric displacement field and the electric current density in the Maxwell equation, we perform a parameter study regarding the electric conductivity parameter . Here, the results show a significant increase of the electric current density and the magnetic field strength with increasing . To investigate the effect of anisotropy on the model, we compared our cubic RVE with an anisotropic cuboid RVE. Depending on the used RVE geometry, the results can vary slightly.

As a final example, we apply our model to a true to scale model of a human femur bone. Here, the results show again a similar behavior for all quantities. Between the two RVEs with the highest volume fraction of cortical bone phase, nearly no reduction of the magnetic field strength can be observed. With decreasing , the differences grow increasingly larger. Compared to the healthy bone (), the bone with late stage osteoporosis () shows a drastic reduction of the magnetic field strength by nearly two thirds. These results show, in which order of magnitude differences between healthy and degenerated bones can be expected, when performing experimental research and sonography for the purpose of early detection of osteoporosis.

For future research, we aim to solve the inverse problem by using an Artificial Neural Network to predict simulation outputs for random microstructures. Here, the network should recover the distribution of cortical bone phase in the macroscopic model from the magnetic field data, thus diagnosing the state of the bone. Additionally, wave propagation in cancellous bone will be investigated in more detail. The comparison of experimental with simulation results could provide further insights. Accurate material parameters could be obtained from the experiments, which then could be used for the simulations. To make precise predictions for experimental setups, it is of great importance to address possible numerical problems of the simulations. The used RVEs have a very coarse mesh resolution and contain sharp edges. While our investigations so far show, that the coarse RVE with hexahedron elements performs best, it is still relevant to investigate in detail how the magnetic field strength is affected for different, smoother RVEs, which model the microstructure of bone in a more realistic way. Another important aspect is to investigate the microscale behavior for RVEs which differ in size and structure of the phases. Depending on the geometry of the used RVE, the simulation results can vary. Thus, for the future we plan to investigate this effect in detail. The usage of different function spaces could improve the results. Finally, our macroscopic models could be extended to include a surrounding medium like air or water, allowing proper decay of the magnetic field.