Theoretical Solutions of Dynamic Responses of Cancellous Bone

Shaohua Wang¹, Jianli Jiang², Xiaobing Lu³

¹The Second People’s Hospital of Yichang, Yichang, China.
²The First People’s Hospital of Yichang, Yichang, China.
³The Institute of Mechanics, Chinese Academy of Sciences, Beijing, China.

DOI: 10.4236/jbm.2019.712013                

Abstract

Human bone may be damaged by impact in the cases of traffic accidents and ship impact. The impact responses of cancellous bone were analyzed based on the two-phase media theory. A direct analytical method is introduced for solving this type of problems. First, flow function and potential function were introduced to decouple the controlling equations. Then direction solving method was used to obtain the solution. The solution is determined by the parameters of a (related with wave speed) and b (related with damping), as well as the boundary conditions. These two parameters a and b determine the propagation speed of the responses along the bone and the attenuation rate. It is shown that the responses: deformation, stress and pressure of the corpus medullae caused by loading, propagate toward the other end when the impact is acted on one end of the bone. The responses are discontinuous during propagate. The discontinuous surface moves with a constant speed. The responses at the cross section increase gradually from the bottom to the top because of the distribution of the loading at the boundary. The solutions can be used as the basis for certification of numerical simulation as well as the design of impact prevention of bone.

Keywords

Cancellous Bone, Impact Load, Responses, Decoupling Method

1. Introduction

In the cases of traffic accidents, ship impact and blast, human bone is easy to be hurt by the impact load [1] [2]. The human tibia and femur are the positions being often damaged. It is practical to study the dynamic responses of the human bone under impact load for either the prevention of bone damage or recovery.

The human tibia and femur are a kind of long tubular bone composed of cancellous bone with high percentage porosity. The cancellous bone is composed of the skeleton and the corpus medullae. Recently, the mechanical properties and responses were concerned because of the requirement of fracture treatment and artificial limb design [3] [4] [5]. The permeability, strength and stiffness of the cancellous bone were studied experimentally. These studies provide the basic data for the following research.

The cancellous bone is composed of a kind of cellular porous material. The corpus medullae are filled fully in the porosity. The mechanical properties of the cancellous bone are closely related with the density [6] [7].

The position and the form of damage by impact are related with the human posture during bearing load. For example, the lower limbs, especially the calcaneus and the femur, are easy to be hurt when a person is standing up. The damage degree is determined by the load level and the mechanical properties of the bone [8]. Under impact, the cancellous bone is compressed and the corpus medullae flow out of the porosity totally or partially. The stress and displacement increase. The pore pressure changes accordingly by the load level and boundary conditions. The displacement rebounds totally or partially and the corpus medullae flow back into porosities once the impact is unloaded [6] [7]. The region of large stress is located near the load end of bone such as the neck of femur under impact while it is located at the backbone under static load [9].

The damage of the cancellous bone has also been studied. Abdel-Rahman et al. [10] showed that the impact can cause the occurrence of the micro fractures at the joint. The study of Dakin et al. [11] showed that the pelvic joint becomes flabby after healing up if one person bears impact load. Yang et al. [12] showed that the strain of ligament and the motion trajectory of the joint will change after impact. The shape of the sacral fracture, propagation of stress wave and strain responses were studied by Quan et al. [13].

However, it is less considered that the cancellous bone should the taken as two-phase media in the former study on the damage of bone during impact. In fact, the properties and behavior of the two media: skeleton and corpus medullae, are very different. The loading acted on the cancellous bone causes the deformation of bone trabecula, which leads to the pressure changes and flow of corpus medullae. Conversely, the flow of corpus medullae affects the deformation of bone trabecula. The relative motion between corpus medullae and bone trabecula forms the resistance [14]. Therefore, it is limited to take the bone as a kind of one-phase media during the analysis of the responses under impact.

In this paper, the theoretical analysis on the dynamic responses of cancellous bone under impact load is carried out to provide analytical solutions as the basis for certification of numerical simulation and design of impact prevention of bone. The solutions can also be the references of the impact experiments of bone. The behavior of cancellous bone is described as poro-elastic media. By introducing the potential and flow functions, the controlling equations are simplified and decoupled. The analytical solutions are obtained by using the direct method. The development of dynamic responses is discussed.

2. Formulation of the Problem

As two-phase media, dynamic load will certainly result in the changes of pressure of corpus medullae, stress and deformation of the cancellous bone. Assumption of a plane strain problem is adopted for simplicity considering the aim of this paper is to present an analytical method for solving this kind of problem. The impact load (p shown in Figure 1) is applied on one end of the cancellous bone. The load is divided into pressure of the corpus medullae and stress of the skeleton. At the other end, the pressure is zero, indicating the corpus medullae can freely flow in or out, while the stress is equal to that at load end in amplitude but inverse in direction which is caused by the person’s weight and limitations. The other surfaces are free of deformation but permeability is forbidden. The change of porosity is small and its gradient can be neglected. Density of each phase is constant. Hooke’s law is adopted to describe the behavior of the skeleton at elastic stage (Figure 1) [6].

3. Controlling Equations

Momentum conservation of the skeleton and the corpus medullae

The skeleton and the corpus medullae affect each other by the deformation, pore pressure and resistance between them. Biot is the first one to propose the mathematical model for two-phase media and the model have been widely used [15]. But the couple mass coefficient is difficult to determine in his model. Thus most researchers adopted the modified model. The controlling equations are similar to Biot’s model in this paper by neglecting the couple mass coefficient but considering the changes of porosity [16].

The skeleton and the corpus medullae should satisfy the mass conversation respectively. Since there are no source items, the mass conservation equations can be written as

$\{\begin{array}{l}\frac{\partial \epsilon \rho }{\partial t}+\frac{\partial \epsilon \rho v_i}{\partial x_i}=0\\ \frac{\partial \left(1-\epsilon \right){\rho }_s}{\partial t}+\frac{{\partial }^2\left(1-\epsilon \right){\rho }_s{u}_i}{\partial t\partial x_i}=0\end{array}$ (1)

Figure 1. Sketch of the problem.

Considering the assumption of constant densities, the following equation can be obtained by combining the two equations in Equation (1)

$\frac{\partial \epsilon v_i}{\partial x_i}+\frac{{\partial }^2\left(1-\epsilon \right)u_i}{\partial t\partial x_i}=0$ (2)

The momentum conservation equations are as follows

$\{\begin{array}{l}\frac{\partial {\sigma }_{ij}}{\partial x_j}-\frac{\partial p}{\partial x_i}-\epsilon \rho g_i-\left(1-\epsilon \right){\rho }_s{g}_i=\left(1-\epsilon \right){\rho }_s\frac{{\partial }^2{u}_i}{\partial t^2}+\epsilon \rho \frac{\partial v_i}{\partial t}\\ -\frac{\partial \epsilon p}{\partial x_i}+p\frac{\partial \epsilon }{\partial x_i}-\frac{{\epsilon }^2\mu }{K}\left(v_i-\frac{\partial u_i}{\partial t}\right)-\epsilon \rho g_i=\epsilon \rho \frac{\partial v_i}{\partial t}\end{array}$ (3)

in which K is the physical permeability, ε is the porosity, u is the displacement of skeleton, v is the velocity of corpus medullae, p is the pressure of corpus medullae, σ_ij is the stress of the skeleton, ρ is the density of corpus medullae, ρ_s is the density of skeleton, g is the gravity acceleration, μ is the viscosity of corpus medullae, i and j are tensor notations and denotes the three directions. To remove the gravity terms we shall replace henceforth ${\sigma }_{ij}$ by ${\sigma }_{ij}-\left(1-\epsilon \right)\left({\rho }_s-\rho \right)gy{\delta }_{ij}$ and p by $p+\rho gy$. Then the above equations become homogeneous

$\{\begin{array}{l}\frac{\partial {\sigma }_{ij}}{\partial x_j}-\frac{\partial p}{\partial x_i}=\left(1-\epsilon \right){\rho }_s\frac{{\partial }^2{u}_i}{\partial t^2}+\epsilon \rho \frac{\partial v_i}{\partial t}\\ -\frac{\partial p}{\partial x_i}-\frac{\epsilon \mu }{K}\left(v_i-\frac{\partial u_i}{\partial t}\right)=\rho \frac{\partial v_i}{\partial t}\end{array}$ (4)

From Equation (4), p can be eliminated, so

$\frac{\partial {\sigma }_{ij}}{\partial x_j}+\frac{\epsilon \mu }{K}\left(v_i-\frac{\partial u_i}{\partial t}\right)=\left(1-\epsilon \right)\left({\rho }_s\frac{{\partial }^2{u}_i}{\partial t^2}-\rho \frac{\partial v_i}{\partial t}\right)$ (5)

Constitutive equations―Hooke’s law

In this study, the skeleton is taken as poroelastic media, e.f., Hooke’s law is suited for describing the mechanical behavior:

${\sigma }_{ij}=G\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+{\delta }_{ij}\frac{2G\nu }{1-2\nu }\frac{\partial u_k}{\partial x_k}$ (6)

Decoupling of equations

Obviously, the above Equations (1)-(5) are difficult to solve directly because the variables are coupled. In other words, the variables are contained in each partial differential equation, which causes the solving process very much complex. Only to decouple these equations into those containing one variable in each equation, the solving process can be simplified. Therefore, the following potential and “flow” functions are adopted:

$\begin{array}{l}\bar{u}=Grad{\phi }_s+Curl{\bar{\psi }}_s\\ \bar{v}=Grad\phi +Curl\bar{\psi }\\ \text{with}\\ Div{\bar{\psi }}_s=0\\ Div\bar{\psi }=0\end{array}$ (7)

Substitute Equation (7) into the mass conservation Equation (2), the following equation can be obtained

$\epsilon \phi +\left(1-\epsilon \right)\frac{\partial {\phi }_s}{\partial t}=0$ (8)

From the second of Equation (4), the relation between the pore pressure and the potential and flow functions can be obtained as follows considering Equation (8)

$p=-\rho \frac{\partial \phi }{\partial t}-\frac{\epsilon \mu }{K}\left(\phi -\frac{\partial {\phi }_s}{\partial t}\right)$ (9)

and

$\frac{\epsilon \mu }{K}\left(\bar{\psi }-\frac{\partial {\bar{\psi }}_s}{\partial t}\right)=-\rho \frac{\partial \bar{\psi }}{\partial t}$ (10)

Similarly taking the divergence and curl of Equation (5) respectively, we can obtain

${\lambda }_m{\nabla }^2{\phi }_s=\left(1-\epsilon \right)\left(1+\frac{1-\epsilon }{\epsilon }\frac{\rho }{{\rho }_s}\right)\frac{{\partial }^2{\phi }_s}{\partial t^2}+\frac{\mu }{{\rho }_{s}K}\frac{\partial {\phi }_s}{\partial t}$ (11)

and

$\frac{G}{{\rho }_s}{\nabla }^2{\bar{\psi }}_s=\left(1-\epsilon +\frac{\epsilon \rho }{{\rho }_s}\right)\frac{{\partial }^2{\bar{\psi }}_s}{\partial t^2}$ (12)

in which ${\lambda }_m=\frac{2G\left(1-\nu \right)}{{\rho }_s\left(1-2\nu \right)}$.

Now, the above equations have been decoupled as Equations (8)-(12). From Equations (11) and (12), ${\phi }_s$ and ${\bar{\psi }}_s$ can be solved solely. Then instituting ${\phi }_s$ and ${\bar{\psi }}_s$ into Equations (8)-(10), $\phi$, p, $\bar{\psi }$ can be obtained. The solving process is presented in the following section.

4. Solution and Discussion

To obtain the solutions, the above equations are further changed. Instituting Equation (8) into Equation (9) to eliminate $\phi$ :

${\nabla }^{2}p=\frac{\rho \left(1-\epsilon \right)}{\epsilon }\frac{{\partial }^2}{\partial t^2}{\nabla }^2{\phi }_s+\frac{\mu }{K}\frac{\partial }{\partial t}{\nabla }^2{\phi }_s$ (13)

Instituting Equation (11) into Equation (13), the controlling equation of p can be obtained:

${\nabla }^{2}p=a\frac{{\partial }^{2}p}{\partial t^2}+b\frac{\partial p}{\partial t}$ (14)

in which $a=\left(1-\epsilon \right)\left({\rho }_s+\rho \frac{1-\epsilon }{\epsilon }\right)\frac{1-2\nu }{2G\left(1-\nu \right)}$, $b=\frac{\mu }{K}\frac{1-2\nu }{2G\left(1-\nu \right)}$

Assuming pressure of corpus medullae caused by the external load at the boundary x = 0 is p₀y (indicating that the pressure linearly distributed along the section), at the other end of x = l, the pore pressure is a constant (here set as zero). At the boundary y = 0 and y = h, $\partial p/\partial y=0$, e.g. the seepage is forbidden. Then the controlling equation of pore pressure and the boundary and initial conditions can be written as follows:

$\{\begin{array}{l}{\nabla }^{2}p=a\frac{{\partial }^{2}p}{\partial t^2}+\frac{\mu }{K}\frac{\partial p}{\partial t}\\ {p|}_{x=0}=p_{0}y,{p|}_{x=l}=0\\ {\frac{\partial p}{\partial y}|}_{y=0}=0,{\frac{\partial p}{\partial y}|}_{y=h}=0\\ {p|}_{t=0}=0,{\frac{\partial p}{\partial t}|}_{t=0}=0\end{array}$ (15)

p can be solved as [17] :

1) if $c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]=0$

$p=A\left(1+\frac{b}{2}\right){\text{e}}^{-\frac{bt}{2}}$ (16)

2) if $c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]>0$, $c>0$

$p=A\left(\frac{b+c}{2c}{\text{e}}^{\frac{-b+c}{2}t}+\frac{c-b}{2c}{\text{e}}^{\frac{-b-c}{2}t}\right)$ (17)

3) if $-c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]<0$, $c>0$

$p=A\left(\cos ct+\frac{b}{2c}\sin ct\right){\text{e}}^{-\frac{bt}{2}}$ (18)

in which

$A=-\frac{p_{0}h}{2\pi }\underset{m=1}{\overset{m\to \infty }{{\displaystyle \sum }}}\frac{1}{m}\sin \frac{m\pi x}{l}+\frac{8p_{0}h}{{\pi }^3}\underset{n=1}{\overset{n\to \infty }{{\displaystyle \sum }}}\underset{m=1}{\overset{m\to \infty }{{\displaystyle \sum }}}\frac{1}{{\left(2n+1\right)}^2}\frac{1}{m}\sin \frac{m\pi x}{l}\cos \frac{n\pi y}{h}$

By Equation (13) and initial condition ${{\phi }_s|}_{t=0}=0$, ${\phi }_s$ can be obtained as:

1) if $c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]=0$

${\phi }_s=-{\displaystyle {\sum }_{m=1}^{m\to \infty }\frac{A}{\frac{\rho \left(1-\epsilon \right)b^2}{4\epsilon }-\frac{b\mu }{2K}}\left(1+\frac{b}{2}\right){\text{e}}^{-\frac{bt}{2}}}$ (19)

2) if $c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]>0$, $c>0$

$\begin{array}{c}{\phi }_s={\displaystyle {\sum }_{m=1}^{m\to \infty }A(\frac{\frac{-b-c}{2c}}{\frac{\rho \left(1-\epsilon \right){\left(c-b\right)}^2}{4\epsilon }-\frac{\left(c-b\right)\mu }{2K}}{\text{e}}^{\frac{-b+c}{2}t}}\\ +\frac{\frac{-b+c}{2b}}{\frac{\rho \left(1-\epsilon \right){\left(c-b\right)}^2}{4\epsilon }-\frac{\left(c-b\right)\mu }{2K}}{\text{e}}^{\frac{-b-c}{2}t})\end{array}$ (20)

3) if $-c^2=b^2-4\left[{\left(\frac{m\pi }{l}\right)}^2+{\left(\frac{n\pi }{h}\right)}^2\right]<0$, $c>0$

${\phi }_s={\displaystyle {\sum }_{m=1}^{m\to \infty }-A\left(A_1\cos ct-A_2\sin ct\right)}$ (21)

in which $A_1=\frac{\rho \left(1-\epsilon \right)}{\epsilon }\left(\frac{3b^2}{4}-c^2\right)-\frac{b\mu }{K},A_2=\frac{\rho \left(1-\epsilon \right)}{\epsilon }\left(-\frac{b^3}{8c}+bc\right)+\frac{\mu }{K}\left(c+\frac{b^2}{4c}\right)$.

By Equation (12), ${\bar{\psi }}_s$ can be obtained in the following way.

The controlling equation and the boundary and initial conditions are rewritten as:

$\begin{array}{l}\frac{G}{{\rho }_s}{\nabla }^2{\bar{\psi }}_s=\left(1-\epsilon +\frac{\rho \epsilon }{{\rho }_s}\right)\frac{{\partial }^2{\bar{\psi }}_s}{\partial t^2}\\ {\sigma }_x=E{\left[\frac{{\partial }^2{\phi }_s}{\partial x^2}-\nu \frac{{\partial }^2{\phi }_s}{\partial x^2}+\left(1-\nu \right)\frac{{\partial }^2{\bar{\psi }}_s}{\partial x\partial y}\right]|}_{x=0}=E\left(1-\nu \right){\frac{{\partial }^2{\bar{\psi }}_s}{\partial x\partial y}|}_{x=0}=p_0\\ E\left(1-\nu \right){\frac{{\partial }^2{\bar{\psi }}_s}{\partial x\partial y}|}_{x=l}=p_0\\ 2{\frac{{\partial }^2{\phi }_s}{\partial x\partial y}+\frac{{\partial }^2{\bar{\psi }}_s}{\partial y^2}-\frac{{\partial }^2{\bar{\psi }}_s}{\partial x^2}|}_{y=0}={\frac{{\partial }^2{\bar{\psi }}_s}{\partial y^2}-\frac{{\partial }^2{\bar{\psi }}_s}{\partial x^2}|}_{y=0}=0\\ {\frac{{\partial }^2{\bar{\psi }}_s}{\partial y^2}-\frac{{\partial }^2{\bar{\psi }}_s}{\partial x^2}|}_{y=h}=0\\ {{\bar{\psi }}_s|}_{t=0}=0,{\frac{\partial {\bar{\psi }}_s}{\partial t}|}_{t=0}=0\end{array}$ (22)

The boundary conditions indicate that the axial stresses of the skeleton at the two ends are equal which are caused by the external load. The shear stresses at the other sides are equal to zero, e.g. the boundaries are free. The initial values are zero.

${\bar{\psi }}_s$ can then be solved as:

$\begin{array}{c}{\bar{\psi }}_s=-\frac{p_{0}hl}{E\left(1-\nu \right)\pi }\underset{n=1}{\overset{n\to \infty }{{\displaystyle \sum }}}\frac{{\left(-1\right)}^{n+1}}{n}\left[1-{\displaystyle {\sum }_{k=0}^{k\to \infty }\frac{8}{{\pi }^2{\left(2k+1\right)}^2}\cos \frac{\left(2k+1\right)\pi x}{l}}\right]\\ \times \sin \frac{n\pi y}{l}\cos dt\end{array}$ (23)

in which $d=\sqrt{D/C}$, $D=\frac{G}{{\rho }_s}{\left(\frac{n\pi }{l}\right)}^2+{\left(\frac{m\pi }{h}\right)}^2$, $C=1-\epsilon +\frac{\rho \epsilon }{{\rho }_s}$.

Based on the above solutions of p, ${\phi }_s$ and ${\bar{\psi }}_s$, it is easy to obtain the displacements and stresses by using of the following equations

$u_{sx}=\frac{\partial {\phi }_s}{\partial x}+\frac{\partial {\bar{\psi }}_s}{\partial y}$

$u_{sy}=\frac{\partial {\phi }_s}{\partial y}-\frac{\partial {\bar{\psi }}_s}{\partial x}$

${\sigma }_x=\lambda {\nabla }^2{\phi }_s+2\mu \left(\frac{{\partial }^2{\phi }_s}{\partial x^2}+\frac{{\partial }^2{\bar{\psi }}_s}{\partial x\partial y}\right)$

${\sigma }_y=\lambda {\nabla }^2{\phi }_s+2\mu \left(\frac{{\partial }^2{\phi }_s}{\partial y^2}-\frac{{\partial }^2{\bar{\psi }}_s}{\partial x\partial y}\right)$

To reduce the article length, the detailed expressions of the displacements and stresses are not shown here.

The responses of the cancellous bone can be investigated by just instituting the values of parameters into the above analytical solutions. The basic parameters are referenced from the literature [6] and shown as follows: The elastic modulus E = 2.5 × 10⁵ N/m², the Poisson’s ratio υ = 0.30, porosity ε = 0.60, density of skeleton ρ_s = 2.0 × 10³ kg/m³, density of corpus medullae ρ_f = 930 kg/m³, permeability k = 4.0 × 10⁻² m/s, and the load amplitude 20 N - 2000 N. The permeability and elastic modulus are changed based on the above values to investigate the effects of parameters.

Figures 2-4 show the responses (stress, displacement and pore pressure) of the cancellous bone along the length at different time. It can be seen that the pressure of the corpus medullae increases near the load end and there is a discontinuous surface and then decreases with distance (Figure 4). The reason is the load causes the increase of pressure, while the other end is permeable freely so the pressure is zero. The stress and displacement propagate also with a discontinuous surface respectively. Before the discontinuous surface, the changes are obvious while after the discontinuous surface the responses are small. At the cross section, the responses increase gradually from the bottom to the top (Figure 5). The reason is that the loading at the boundary is linearly distributed along the cross section.

Figure 6 and Figure 7 show the dynamic responses with distance and parameters of a and b at the time t = 0.18 s. It is shown that the dynamic responses

Figure 2. Development of stress along the length.

Figure 3. Development of displacement along the length.

Figure 4. Development of pore pressure along the length.

Figure 5. Development of displacement along cross section.

(a) (b)

Figure 6. Effects of parameters a. (a) Pressure of corpus medullae; (b) Stress of skeleton.

(a) (b)

Figure 7. Effects of parameters b. (a) Pressure of corpus medullae; (b) Stress of skeleton.

increase with the decrease of the parameters a. It can be seen from Equation (14) that 1/a is the square of the wave speed. Thus the decrease of a means the disturbance caused by the loading propagates speed increase and, the propagation distance increases at the same time. Accordingly, the attenuation of disturbance in the bone increases. In other words, the responses (displacement, strain pressure and stress) decrease with decrease of a.

The increase of b means the increase of damping, so the responses decrease with time and distance. It can be seen from Figure 7 that the responses decrease more than 10 times with the value of b changing from 13.5 m/s² to 1.35 × 10⁻³ m/s². However, the propagation speed keeps as a constant, which is determined by only the value of a.

5 Conclusions

Dynamic responses of the cancellous bone under impact load were analyzed theoretically. Analytical solutions were presented by using a decoupling method. First, based on liquid-solid media theory, two dimensional two-phase controlling equations were obtained. Flow function and potential function were then introduced to decouple the controlling equations and the solutions were solved.

The solutions show that the responses (pressure of the corpus medullae, strain, stress and deformation) increase near the load end and then decrease with distance. There is a discontinuous surface. The dynamic responses increase with the decrease of the parameters a. The increase of b causes the dynamic responses decrease with time and distance. The responses at the cross section increase gradually from the bottom to the top similar to the distribution of the loading at the boundary.

According to the obtained solutions, the end near the loading, such as the lower end of tibia and femur, is the priority position for protection. In experiments, the loading form and distribution at the loading end and the boundary condition at the other end should be strictly controlled because these factors greatly affect the dynamic responses of the cancellous bone. More sensors should be placed near the loading end to measure fine information.

Nomenclature

ρ_w: density of fluid

ρ_s: density of solid

$u_i$ : displacements of solid in three directions

$v_i$ : velocities of fluid in three directions

ε: porosity

p: fluid pressure

${\sigma }_{ij}$ : effective stress

g: gravitational acceleration

G: shear modulus

ν: Poisson ratio

${\phi }_s$ : potential functions of solid

${\psi }_s$ : flow functions of solid

$\phi$ : Potential functions of water

$\psi$ : flow functions of water

K: the physical permeability ( $={\rho }_{w}g/\left(\mu k\right)$ )

k: the Darcy’s permeability

g: the gravity acceleration

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Xu, S.C., Huang, Y., Fu, L.H., et al. (2016) Comparison and Analysis of Occupant Tibia Injury Criteria under Explosion Shock Loading. Journal of Automotive Safety and Energy, 7, 371-376 (In Chinese).
[2]	Huang, J.S. and Zhou, J.P. (2008) Shock Injury Criteria and Tolerance of Personnel Ship Board. Beijing Biomedical Engineering, 27, 300-304 (In Chinese).
[3]	Li, D.Y. and Chen, H.B. (2001) Test of Fluid Penetrability through Cancellous Bone and Its Theoretical Prediction. Journal of Chongqing Normal University (Natural Science Edition), 18, 19-22 (In Chinese).
[4]	Ou, Y.J., Yang, G.T., Zhao, L.M., Zhao, Y.G., Zhu, Q.G. and Zhong, S.Z. (1996) Impact Response of Human Lumbar Vertebral Trabecular Bone. Journal of Biomedical Engineering, 13, 29-33 (In Chinese).
[5]	Linde, F., Norgaad, I., Hvid, A., Odgaard, A. and Soballe, K. (1991) Mechanical Properties of Trabecular Bone Depend on Strain Rate. Journal of Biomechanics, 3, 296-301. https://doi.org/10.1016/0021-9290(91)90305-7
[6]	Li, D.Y., Yao, G.W., Liu, Z.F. and Zhang, X.W. (2000) Numerical Analysis of Impact Response of Cancellous Bone in Human Tibia. Journal of Chongqing University (Natural Science Edition), 23, 82-86 (In Chinese).
[7]	Liu, Z.F., Chen, H.B., Li, D.Y. and Wang, Z.G. (2000) One-Dimensional Impact Dynamic Response of Cancellous Bone. Journal of Biomedical Engineering, 17, 266-269 (In Chinese).
[8]	Huang, J.S., Hua, H.X. and Zhou, J.P. (2008) Lower Extremities Axial Compression Injury Characteristic of Personnel Shipboard by Simulating Ship Shock Motions. Chinese Journal of Biomedical Engineering, 27, 410-415 (In Chinese).
[9]	Jiang, H.B. and Ge, S.R. (2006) Numerical Modeling of Impact Character in Human Femur. Chinese Journal of Clinical Rehabilitation, 10, 114-117 (In Chinese).
[10]	Abdel-Rahman, E.M. and Hefzy, M.S. (1998) Three-Dimensional Dynamic Behaviour of the Human Knee Joint under Impact Loading. Medical Engineering & Physics, 20, 276-290. https://doi.org/10.1115/1.1372321
[11]	Dakin, G.J., Arbelaez, R.A., Molz, F.J., Alonso, J.M. and Eberhardt, A.W. (2001) Elastic and Viscoelastic Properties of the Human Pubic Symphysis Joint: Effects of Lateral Impact Joint Loading. Journal of Biomechanical Engineering, 123, 218-226. https://doi.org/10.1115/1.1372321
[12]	Yang, J.K. and Kajzer, J. (1997) Computer Simulation of Impact Response of the Human Knee Joint in Car-Pedestrian Accidents. Doktorsavhandlingar vid Chalmers Tekniska Hogskola, 1320, 15.
[13]	Quan, R.F., Yang, D.S. and Wang, Y.J. (2005) Experimental Study of the SacrumFracture under the Dynamics Clashing Load. Chinese Journal of Clinical Anatomy, 23, 188-192.
[14]	Li, D.Y. and Chen, H.B. (2001) Numerical Analysis on the Viscoelastic Characteristic of Cancellous Bone. Journal of Chongqing University (Natural Science Edition), 24, 91-94.
[15]	Biot, M.A. (1956) Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid: II: High Frequency Range. The Journal of the Acoustical Society of America, 28, 168-191. https://doi.org/10.1121/1.1908241
[16]	Lu, X.B., Tan, Q.M., Cheng, C.M., et al. (2004) Liquefaction and Displacement of Saturated Sand under Vertical Vibration. Acta Mechanica Sinica, 20, 96-105. https://doi.org/10.1007/BF02493578
[17]	Liang, K.M. (2010) Mathematical Methods for Physics. Fourth Edition, High Education Press, Beijing.