## Abstract

Recently, micro-sphere-based methods derived from the angular integration approach have been used for excluding fibres under compression in the modelling of soft biological tissues. However, recent studies have revealed that many of the widely used numerical integration schemes over the unit sphere are inaccurate for large deformation problems even without excluding fibres under compression. Thus, in this study, we propose a discrete fibre dispersion model based on a systematic method for discretizing a unit hemisphere into a finite number of elementary areas, such as spherical triangles. Over each elementary area, we define a representative fibre direction and a discrete fibre density. Then, the strain energy of all the fibres distributed over each elementary area is approximated based on the deformation of the representative fibre direction weighted by the corresponding discrete fibre density. A summation of fibre contributions over all elementary areas then yields the resultant fibre strain energy. This treatment allows us to exclude fibres under compression in a discrete manner by evaluating the tension–compression status of the representative fibre directions only. We have implemented this model in a finite-element programme and illustrate it with three representative examples, including simple tension and simple shear of a unit cube, and non-homogeneous uniaxial extension of a rectangular strip. The results of all three examples are consistent and accurate compared with the previously developed continuous fibre dispersion model, and that is achieved with a substantial reduction of computational cost.

## 1. Introduction

Collagen fibres in soft biological tissues provide the overall stiffness and strength of the material. The latest imaging techniques, such as second-harmonic generation, have enabled detailed visualization of the underlying microscopic constitution of biological tissues such as arterial walls [1–3], carotid arteries [4], the myocardium [5,6], the pericardium [7], articular cartilage [8,9] and other tissues. The collagen fibres in these tissues may be dispersed randomly in space, in a certain pattern such as predominantly in a particular direction [10], as a rotationally symmetric dispersion about a mean direction, or as the recently observed non-symmetric dispersion in arterial walls [1,11], or other arrangements. Continuum constitutive laws that account for such underlying material microstructure have been proposed and employed extensively in the last few decades to model the mechanical response of these fibrous tissues (see [11–13]). In particular, constitutive laws that incorporate the three-dimensional (3D) fibre dispersion in fibrous tissues have attracted a lot of interest in the last decade. However, the precise description of the 3D fibre dispersion in a constitutive equation for the modelling of fibrous tissues poses formidable challenges even when considerable simplifications and idealizations are made.

Currently, there are two main approaches for modelling the dispersed fibre distributions in a constitutive equation, namely the ‘generalized structure tensor’ and the ‘angular integration’ (AI) approaches [11]. In the AI approach [14], the strain energy of a single collagen fibre is assumed to be a function of the fibre stretch. The fibre dispersion in the tissue is incorporated into the strain-energy function by an integration of the single fibre strain energy over all the fibre directions weighted with a continuous probability density function (PDF). As the terminology ‘AI’ is rather imprecise and does not explicitly mention fibre dispersion, in this paper instead of AI henceforth we use the terminology ‘continuous fibre dispersion’ (CFD). If the fibre dispersion is incorporated as a summation of a finite number of discrete fibre contributions, then we refer to this as the ‘discrete fibre dispersion’ (DFD) method.

The CFD approach has attracted a lot of interest, and since it was introduced in 1983 there have been numerous studies based on this approach; see [15] and references therein. Briefly, the fibre contribution *Ψ*_{f} to the strain-energy function of a tissue per unit reference volume is obtained by integrating the weighted strain energy *Ψ*_{n}(*λ*) for each fibre direction **N** over the unit sphere , i.e.
1.1where *λ* is the fibre stretch in the direction **N**, the PDF *ρ*(**N**) represents the probability density of the fibre in that direction in the reference configuration, and d*Ω* is the solid angle on the sphere. Note that we are considering one fibre family of the same type embedded in the matrix material. The original study [14] allowed for different fibres to have different properties, but here we consider all the fibres to have the same properties. If additional fibre families of the same or different type exist in the tissue, they can be included additively [16] with corresponding fibre volumetric ratios. The PDF in (1.1) must satisfy the normalization condition
1.2

The computational implementation of the CFD approach (1.1) requires two-dimensional integration over a unit sphere at each Gauss point during a finite-element analysis. In general, due to the complex natures of *ρ* and *Ψ*_{n}, analytical solutions of the integration in (1.1) only exist for some special cases. Frequently, it is evaluated by numerical methods such as
1.3where **N**_{n} and *w*_{n}, respectively, for *n* = 1, …, *m*, are integration points (orientations) and weights defined by the particular integration scheme over the unit sphere, *m* is the number of integration points, and *λ*_{n}, *n* = 1, …, *m*, are the stretches associated with the integration points. In a recent study [17], large errors in the stress–strain result have been observed with some of the commonly used numerical integration schemes over the unit sphere, such as the method in [18]. Verron [17] concluded that the errors observed could be ‘partially explained by the inability of standard methods to handle non-smooth functions even with a large number of integration points’. Out of the 20 numerical integration schemes over the unit sphere studied in [19], the best two integration schemes [18,20] have been found to be inaccurate for large deformation problems in a more recent study [21].

The only integration scheme found to be accurate is the FM900 scheme [22] which involves 900 integration points. However, the authors tested only the integration scheme under uniaxial and biaxial extensions. It is unknown whether this integration scheme would work for more complex loading conditions. Furthermore, in this scheme a few integration points have negative weights [22]. In the extreme case, if only fibres at these integration points are under tension, then the micro-sphere-based model would produce a negative stress even for fibres under tension. When the exclusion of compressed fibres is accounted for in finite-element analyses [23], the integration error is even larger with low-order integration schemes. Although the authors of [23] obtained a homogeneous stress distribution in the circumferential direction of the artery by using high-order integration schemes with 368 or 600 integration points, it is still unknown whether the tested high-order integration schemes are accurate enough for other large deformation problems.

Owing to the waviness and slenderness of the fibres, it is often assumed that they do not contribute to the strain energy of the tissue when loaded under compression [24]. To incorporate such tension–compression behaviour of the fibres in the strain-energy function, *Ψ*_{n}(λ_{n}) is often set to zero in (1.3) when the fibre stretch λ_{n} at a particular integration point **N**_{n} is less than one [23,25]. Setting *Ψ*_{n}(λ_{n}) to zero within a sub-domain of sphere renders some strain-energy functions or their derivatives discontinuous. However, existing numerical integration schemes are often proposed for continuous functions such as polynomials over the entire unit sphere, as in the widely used numerical integration schemes described in [18,20], which were proposed for polynomial functions of certain degrees. In addition, setting *Ψ*_{n}(λ_{n}) to zero within a sub-domain of the sphere is equivalent to using some of the integration points and weights for numerical integration of a function over only its complement. This treatment is questionable because the accuracy of the integration scheme may not be maintained for polynomial functions up to a certain degree locally over a sub-domain. Besides, if a realistic PDF measured from experiment is adopted, for example [3], the integrand could become more complex, and it could even be a discontinuous function.

The computational cost for the numerical evaluation of the fibre contribution over a subset of the unit sphere could be substantially reduced by using a recently proposed general invariant [26] or by using parallel computing platforms such as OpenMP [27]. Besides using parallel computing techniques on a high performance computing cluster, another possibility for modelling the tension–compression behaviour of fibres is to use the DFD approach, which could reduce computational time while maintaining accuracy, as shown in this study. Although it is not practical to count the actual number of fibres in a tissue, the concept of the DFD approach has already been applied in several areas. For example, a 3D DFD model with 30 fibre directions was proposed for the modelling of human skin [28], and a bivariate normal distribution was adopted to describe the fibre arrangement. With appropriate parameters, the model was able to represent various fibre arrangements such as aligned fibres or rotationally symmetric fibre dispersions. Although only 30 fibre directions were determined by using the parameters of the bivariate normal distribution, the model was able to capture the fibre undulation and fibre compression–tension behaviour of skin tissue under simple loading scenarios.

Similarly, a DFD model that consists of six weighted fibre bundles was proposed in [29], with all the fibres in each bundle oriented in the same direction. The volume fraction of fibres in each bundle was given by the weights. The six weights were then used as structure parameters in the model fitting. The case of six equal weights represents a 3D uniform fibre dispersion. The strain energy of each fibre bundle was assumed to be a function of the fibre stretch in the bundle direction, and that of all the fibres was determined by a weighted summation of the contributions from all six fibre bundles. Any fibre bundle under compression can be easily excluded from the fibre contribution to the strain-energy function. Although only six fibre bundles appeared in the constitutive equation, the model was able to approximately capture the mechanical response of rabbit skin, porcine skin, porcine aortic valve cusp and rat myocardial tissues. However, one limitation of this model, because of the discretization, is that the material response is not always isotropic even when the six weights are the same. To overcome this unphysical prediction, the authors improved the model by using a generalized strain invariant [30], and more recently by increasing the number of equally weighted fibre bundles [31] because the fitting of a large number of different weights is not practical.

In this study, we propose a systematic approach for determining the discrete fibre density of each fibre bundle directly from the fibre PDF in a straightforward way. Briefly, we first discretize the unit hemisphere into a finite number of spherical triangles and then compute a representative fibre direction at the centroid of each spherical triangle. The discrete fibre density for each representative fibre direction is then determined by numerical integration of the continuous fibre PDF over the corresponding spherical triangle. With this treatment, we can easily incorporate the tension–compression behaviour of the fibres over each spherical triangle. Through three numerical examples we demonstrate that the current model is able to accurately match predictions of the previously developed CFD model when compressed fibres are excluded, with a substantial reduction of computational cost.

This paper is structured as follows. In §2, we present the continuum mechanical framework for the proposed DFD model, including the discretization scheme for the unit sphere, the strain-energy function, and the Cauchy stress and elasticity tensors for 3D fibre dispersions. In §3, we introduce the von Mises distribution for the PDF specialized to a rotationally symmetric dispersion for illustration. Next, we provide a guideline for the implementation of the proposed discretization scheme in a finite-element programme. Additionally, specific forms of the strain-energy function associated with the fibres are provided. The theory introduced in §2 is then applied to several representative examples with the aim of demonstrating the efficacy and efficiency of the proposed model. Finally, §4 summarizes the proposed computational approach and suggests some future directions.

## 2. Discrete fibre dispersion model

In this section, we briefly present the continuum mechanical framework for the proposed DFD model including kinematics, the strain-energy function, and the corresponding Cauchy stress and elasticity tensors, which are introduced in a decoupled form.

### 2.1. Kinematics

The deformation map **x** = **χ**(**X**) transforms a material point **X** in the stress-free reference configuration into a spatial point **x** in the deformed configuration. The deformation gradient is defined as **F**(**X**) = ∂**χ**(**X**)/∂**X**, and its determinant represents the local volume ratio at point **X**, with *J* ≡ 1 for a strictly incompressible material. Following the multiplicative decomposition of the deformation gradient [32,33], we decouple **F** into a volumetric (dilatational) part *J*^{1/3}**I** and an isochoric (distortional) part , with . Based on **F** we define the right Cauchy–Green tensor as **C** = **F**^{T}**F** and its isochoric counterpart as , with the corresponding first invariants defined by
2.1respectively. Let **N** be a fixed vector in the reference configuration; then **C** : **N** ⊗ **N**, denoted *I*_{4}, represents the square of the stretch in the direction **N**. Its isochoric counterpart is denoted . Hence,
2.2Let us now introduce unit Cartesian basis vectors **E**_{1}, **E**_{2}, **E**_{3} and then express **N** in terms of spherical polar angles *Θ* and *Φ* relative to **E**_{1}, **E**_{2}, **E**_{3} such that
2.3

### 2.2. Strain-energy function

We assume that the 3D fibre dispersion inside the matrix material can be described by an integrable function *ρ*(**N**), which we now write as *ρ*(*Θ*, *Φ*), defined over the unit hemisphere
2.4Ideally, *ρ* should be determined by imaging analysis [3].

If the strain energy associated with an individual fibre direction **N** is described by a function *Ψ*_{n}(*I*_{4}), where *I*_{4} is defined in (2.2)_{1}, we require
2.5and, following the CFD approach [34], the strain-energy function of all fibres in the reference configuration can be written over the unit hemisphere as
2.6

Inspired by the discrete nature of fibres dispersed within the ground matrix, and the fact that the number of fibres is finite, we may treat the fibres in a discrete manner. Ideally, the dispersion of fibres and the number of fibres should be determined by experimental measurement. However, due to the large number of fibres, their actual number and orientations at a specific point in the tissue may not be possible to determine accurately. Thus, in this study, we first discretize the unit sphere into a finite number of elementary areas . An example of such a discretization with spherical triangles is shown in figure 1*b*. We then identify representative fibre angles (*Θ*_{n}, *Φ*_{n}) (associated with the centroids of the spherical triangles) for each elementary area and use these angles to represent all the fibres distributed within . Thus, the number of representative fibre angles is equal to the number of spherical triangles over the unit sphere, as can be seen in figure 1. Note that only half of the representative fibre directions are needed because of symmetry. If the area of is chosen to be very small, then the variation of the fibre directions within becomes negligible. In the extreme case when shrinks to a point, then the fibre direction is unique. The normalized number of fibres within each elementary area can then be determined from a set of *m* discrete fibre densities *ρ*_{n} defined by
2.7In fact, *ρ*(*Θ*, *Φ*) could be a discontinuous function over the hemisphere because the integration in (2.7) can be carried out over continuous sub-domains of the hemisphere. Similarly, the discrete fibre densities *ρ*_{n} satisfy the normalization condition over the unit hemisphere, i.e.
2.8which is the discrete counterpart of (1.2) on a hemisphere, where *m* denotes the number of representative fibre directions embedded in the matrix at a specific material point. The value of *ρ*_{n} depends on the area of and the fibre dispersion. The areas of the s on the sphere can vary or be nearly equal, as shown in figure 1*b*. If necessary, regions with higher fibre density, for example, the region near the mean fibre direction **M** in figure 1*b*, can be discretized with smaller . However, in this study, for the purpose of demonstration, we choose the triangular discretization shown in figure 1*b*. With a discretized hemisphere, we then re-define the strain-energy function (2.6) for all the fibre contributions as
2.9where *I*_{4n} = **C** : **N**_{n} ⊗ **N**_{n} and **N**_{n} is defined at the centroid of each spherical triangle via (2.3) with *Θ* = *Θ*_{n} and *Φ* = *Φ*_{n}; see the black arrows in figure 1*b*. Next, in order to exclude fibres under compression within a dispersion, we define *Ψ*_{n} as
2.10A representative plot of the strain-energy function associated with one fibre versus *I*_{4n} is shown in figure 2. As can be seen, the strain energy is only non-zero when *I*_{4n} > 1. Note that the first derivative of the strain-energy function with respect to *I*_{4n} is continuous and equal to zero at *I*_{4n} = 1, and hence, with the requirement (2.5), *f*(1) = *f*′(1) = 0.

For efficient computational implementation, we write the strain-energy function in a decoupled form, namely *Ψ* = *Ψ*_{vol} + *Ψ*_{iso}, where *Ψ*_{vol} is the volumetric strain-energy function and *Ψ*_{iso} is the isochoric part of the strain-energy function associated with one family of embedded fibres, given by
2.11where *Ψ*_{g} denotes the isochoric strain energy of the ground matrix, which is assumed to be isotropic and to depend only on . Then, from (2.9) and (2.11) we have
2.12where . For strictly incompressible materials, we have .

As our focus is on incompressible materials, the volumetric strain-energy function is used as a penalty function, and it is convenient to adopt a form for *Ψ*_{vol} given in FEAP [35], i.e.
2.13where *K* is a penalty parameter. The derivations of the volumetric parts of the stress and elasticity tensors are straightforward and have been well documented [12,36]. Hence, in the following, we only derive the isochoric parts of the stress and elasticity tensors.

### 2.3. Cauchy stress tensor

The so-called fictitious second Piola–Kirchhoff stress tensor is determined by differentiation of the isochoric strain-energy function with respect to . Thus, from (2.12) we obtain
2.14where **I** is the second-order unit tensor, , , and is given in (2.1)_{2}. Push-forward of yields the fictitious Cauchy stress tensor as
2.15where is the modified left Cauchy–Green tensor, and . The isochoric Cauchy stress tensor *σ*_{iso} is then determined as
2.16where is the fourth-order Eulerian projection tensor, and the symmetric fourth-order unit tensor is defined in component form by , where *δ*_{ad} is the Kronecker delta.

### 2.4. Elasticity tensor

The fourth-order fictitious elasticity tensor in the Lagrangian description is obtained by differentiation of with respect to followed by multiplication by a factor *J*^{−4/3}. This, with (2.14), gives
2.17where
2.18Note that for the neo-Hookean model, for which
2.19we have , where the constant *μ* (>0) is the shear modulus. We adopt this model here because it has been shown in [37] that it is sufficient to use a neo-Hookean model for the ground matrix.

A push-forward operation of with **F** yields the fictitious elasticity tensor in the Eulerian description, i.e.
2.20where the neo-Hookean model has been used. If , then an additional term should be included in (2.20)_{1}. Finally, with (2.20), we obtain the resulting isochoric part of the elasticity tensor in the Eulerian description, i.e.
2.21which is needed for the finite-element implementation together with the volumetric part [36].

## 3. Computational aspects and representative examples

### 3.1. Choices of fibre distribution and fibre model

We have implemented the proposed DFD model (2.12) in the general purpose finite-element analysis program FEAP [35] at the integration point level. Note that in the continuum mechanical framework described in §2, we have not specified any particular form of the fibre PDF *ρ*(*Θ*, *Φ*). Thus, our model is applicable to any type of fibre dispersion, symmetric or non-symmetric. Here, for illustration of the method, we choose the rotationally symmetric fibre dispersion described by the von Mises distribution
3.1where *b* is a concentration parameter describing how closely the fibres are distributed around the mean fibre direction **M**, erfi(*x*) =− i erf(i*x*) denotes the imaginary error function and erf(*x*) is the standard error function. This distribution will be used in all the numerical examples.

On substituting (3.1) into (2.7), we obtain a set of *m* discrete fibre densities *ρ*_{n} for the representative fibre directions **N**_{n}, *n* = 1, …, *m*. A general guideline for discretizing the hemisphere and for determining the *ρ*_{n} is given in the accompanying box (algorithm 1).

For the numerical examples, we also need to specify the single fibre strain-energy function . For example, one possible choice is the exponential form proposed in [24] that can capture the highly nonlinear behaviour of soft tissues, i.e.
3.2where *k*_{1} is a positive material parameter with the dimension of stress, and *k*_{2} is a positive dimensionless parameter. For comparison with existing results, we have also implemented the quadratic form (the standard fibre reinforcing model [40]) of , i.e.
3.3where *ν* is a non-negative material constant with the dimension of stress. It is easy to verify that both relations (3.2) and (3.3) satisfy *f*(1) = *f*^{′}(1) = 0, and for . On substituting either (3.2) or (3.3) into the isochoric Cauchy stress tensor (2.16) and the Eulerian fictitious elasticity tensor (2.20), we obtain the specific forms of the Cauchy stress and Eulerian elasticity tensors.

### 3.2. Representative numerical examples

To illustrate the performance of the proposed DFD model, we present three representative examples, specifically the homogeneous simple tension and simple shear of a unit cube with a 3D fibre dispersion, and the inhomogeneous extension of a rectangular strip with a mean fibre direction which does not coincide with the loading direction. For each of the three examples, we assume that the material is incompressible. To enforce the incompressibility condition, we use the augmented Lagrangian method [41] in FEAP. In each example, the geometry of the model is discretized with 8-node hexahedral mixed Q1/P0 elements, and each problem is solved by using the Newton–Raphson method. The proposed model is implemented at the Gauss point level. During the finite-element analysis, the frequency by which the stress and elasticity tensors are updated depends on the actual problem. The finite-element solutions of the first two examples are verified by Matlab or Mathematica solutions [42,43] based on analytical expressions obtained using the corresponding CFD model, as previously reported in [16]. The last example is verified by comparing with the numerical solution obtained by the CFD model, which is reported in [34].

### 3.2.1. Simple tension

Here, we consider a uniaxial tension test of an incompressible unit cube, as described in §3.1 of [16]. Briefly, we consider one family of fibres with a rotationally symmetric dispersion embedded in an isotropic matrix material. The mean fibre direction is taken to be aligned with the loading direction **E**_{3}, as depicted in figure 3. Thus, only a subset of fibres around the mean fibre direction over the hemisphere and the mean fibre direction itself are under tension. The fibres under tension form a right circular cone [44], and the fibres outside the cone are under compression. We apply a displacement boundary condition on the top face of the unit cube to impose a stretch of 1.2, as indicated in figure 3.

For this particular example, *I*_{4}(**N**) for any fibre direction **N** is given by
3.4where *λ* is the stretch in the loading direction. Note that *I*_{4}(**N**) is independent of *Φ* in this special case. The general form of the Cauchy stress tensor * σ* for this problem is given by Li

*et al*. [16] 3.5where for this special case we define the integration domain as

*Ω*= {(

*Θ*,

*Φ*) |

*Θ*∈ [0,

*π*/2],

*Φ*∈ [0, 2

*π*],

*I*

_{4}> 1},

*p*is the Lagrange multiplier,

**n**=

**F**

**N**, and, because of symmetry, the PDF

*ρ*(

*Θ*,

*Φ*) reduces to 3.6As given in [16], the uniaxial Cauchy stress

*σ*≡

*σ*

_{33}in the loading direction

**E**

_{3}is 3.7where

*α*and

*β*are defined over the domain

*Σ*= {

*Θ*∈ [0,

*π*/2] |

*I*

_{4}> 1} as 3.8The Cauchy stress–stretch result (3.7) was implemented in Matlab, and we obtained solutions of this problem with material parameters

*μ*= 1.64 kPa,

*k*

_{1}= 5.63 kPa and

*k*

_{2}= 14.25 [16]. The relationship between the Cauchy stress and the stretch in the loading direction is shown in figure 4 for

*b*= 0.01 and

*b*= 5 (solid curves). For comparison, we have also plotted the finite-element solutions (open circles) by using the proposed DFD model in FEAP with the same material parameters as in our previous study, and with

*m*= 640 discrete fibre directions, which is enough to obtain very accurate results (figure 4). Note that the implementation of the proposed model in FEAP is based on the continuum mechanical framework of §2. As can be seen in figure 4, a very good match between the Matlab and the finite-element solutions for different concentration parameters has been obtained. This indicates that the DFD model is able to predict the same result as the CFD model for this particular problem.

### Remark.

For reason of completeness we now address the numerical results obtained by using the micro-sphere-based approach. We rewrite the coefficients *α* and *β* in (3.7) over a subset of the whole unit sphere, i.e. *Ω* = {(*Θ*, *Φ*)|*Θ* ∈ [0, *π*], *Φ* ∈ [0, 2*π*], *I*_{4} > 1}. Thus,
3.9

We then adopt the numerical integration scheme FM900 of [22] with 900 integration points, as discussed in [19]. To exclude fibres under compression, we use only the angles (*Θ*_{i}, *Φ*_{i}) and weights *w*_{i} of the FM900 scheme if the fibre stretches in those directions are greater than one. Thus,

3.10

where *ω* represents the set of angles that lie inside *Ω*. We have carried out the numerical integration in Matlab and obtained the results for the Cauchy stress versus stretch with the concentration parameters *b* = 5 and *b* = 0.01. The relative error of the numerical result obtained by using the micro-sphere-based approach (3.10) with the FM900 integration scheme [22] is small compared with the Matlab solution based on the CFD model [16]. However, the numerical integration scheme with 21 points used in [18] is not able to reproduce the curves in figure 4. Similar results for the uniaxial and biaxial tension tests have also been observed for the case when the compressed fibres are not excluded [21]. Further numerical investigations are needed to verify whether this integration scheme is applicable in a finite-element analysis with complex loading conditions.

### 3.2.2. Simple shear

Similarly to the previous section, here we test the capability and efficiency of the proposed DFD model by subjecting the same unit cube to a simple shear deformation, as described in [16]. Briefly, with reference to figure 5, all the nodes on the bottom face of the cube in the (**E**_{1}, **E**_{2})-plane are constrained in all three translational degrees of freedom, and a horizontal displacement in the **E**_{1} direction is applied on the top face. We take the mean fibre direction **M** to be at 135° clockwise from the **E**_{3} direction in the (**E**_{1}, **E**_{3})-plane in the reference configuration, as illustrated on a cross-section of the cube in figure 5. This orientation was chosen so that the exclusion of fibres under compression has a significant influence on the resulting shear stress.

The Cauchy shear stress component *σ*_{13} in the (**E**_{1}, **E**_{3})-plane is given by Li *et al*. [16]
3.11where *α* and *γ* are defined by

3.12

and . For this particular example, the invariant *I*_{4}(**N**) has the explicit form
3.13We implemented the result (3.11) in Mathematica and obtained the solution for *σ*_{13} as a function of the amount of shear *c*. For this problem, we used the parameters *μ* = 7.64 kPa, *b* = 1.08, *k*_{1} = 996.6 kPa and *k*_{2} = 5.249 [16]. For comparison, we have also computed the case for *b* = 2. In figure 6, we have plotted the solutions from Li *et al*. [16] and the numerical results by using the proposed DFD model with the same material parameters as in our previous study and with *m* = 640 discrete fibre directions. As can be seen, the numerical results match very well with the corresponding results from Li *et al*. [16]. For higher values of *m*, we obtained the same stress response, but for smaller values of *m* the prediction of the stress response became less accurate.

### 3.3. Extension of a rectangular strip

In the previous two sections, we have verified the proposed DFD model with a unit cube under simple tension and simple shear. In this example, we consider the uniaxial extension of a rectangular strip with loading direction different from the mean fibre direction. For this problem, the deformation field is non-homogeneous, as distinct from the two previous examples. We choose this example because a solution of this problem obtained by using the corresponding CFD model has already been presented in [34]. Hence, we compare the results of the proposed DFD model with the existing result of the CFD model.

The geometry, the boundary conditions and the 3D fibre dispersion of the rectangular strip are described in [34]. However, here we align the longitudinal direction with the **E**_{3} axis instead of the **E**_{1} axis, as shown in figure 7. Briefly, a rectangular strip of 10 × 4 × 1 mm is discretized with 320 hexahedral elements. All nodes on the bottom face of the strip are constrained in the **E**_{3} direction. In addition, to prevent rigid body translation, we constrained the centre node of the bottom face in the **E**_{1} and **E**_{2} directions. Furthermore, the **E**_{2} degree of freedom at node A on the bottom face of the strip (figure 7*b*) is also constrained to prevent rigid body rotation about the **E**_{3} axis. A 3D rotationally symmetric fibre dispersion is assumed with the mean fibre direction aligned at 60° from the **E**_{3} direction counterclockwise in the (**E**_{1}, **E**_{3})-plane (figure 7*b*). A displacement boundary condition is applied on the top face of the strip to impose a stretch of λ = 1.4 in the **E**_{3} direction.

The numerical result for this problem obtained by using the CFD model with the quadratic fibre strain-energy function (3.3) has been provided in §4.2.1 of [34], and is reproduced here in figure 7*a*. In figure 7*b*–*e*, for comparison, the distributions of the Cauchy stress component *σ*_{33} in the deformed configuration obtained by using the DFD model with different values of *m* are also plotted. As can be seen, the numerical result by the DFD model approaches that of the CFD model as *m* increases. As shown in figure 7*b*, the DFD result with *m* = 4000 is essentially identical to the CFD result (figure 7*a*). This shows that we are able to approach the result of the CFD model with the proposed DFD model even when exclusion of compressed fibres is considered.

An important advantage of the DFD model is the substantial reduction in processing time. For example, the numerical simulation with the CFD model was completed in 3.3 h on a typical Windows computer with an Intel^{®} Core^{TM} i7-4770 Processor and 16 GB of memory. With the proposed DFD model, and a discretization density of *m* = 4000, the same problem was computed in only 53 s—a speed-up of 224! Certainly, a further reduction of the processing time can be achieved with a smaller value of *m*, but the accuracy would then be reduced.

## 4. Concluding remarks

In this study, we have proposed an efficient DFD model capable of excluding fibres under compression in the modelling of the nonlinear behaviour of fibrous tissues. We have introduced a systematic approach for discretizing the unit sphere into a finite number of spherical triangles and for computing the discrete fibre densities over each elementary spherical triangle. The discrete fibre densities associated with each spherical triangle are analogous to the weighting factors used in the micro-sphere-based approach but have a direct physical meaning. We were then able to conceive of and formulate a strain-energy function in terms of contributions from each fibre bundle distributed within the corresponding spherical triangle. This discrete treatment of fibre contributions, inspired by the discrete nature of fibres embedded in fibrous tissues, allows us to exclude any fibres under compression within a 3D dispersion in a rather straightforward way. The resulting nonlinear elastic constitutive model depends only on those fibre bundles under tension and the matrix material under arbitrary 3D finite deformations.

With a volumetric/isochoric split of the deformation gradient, we presented analytical expressions of the corresponding Cauchy stress and elasticity tensors in decoupled forms especially suitable for finite-element implementation. By using the augmented Lagrangian method for enforcing material incompressibility and a mixed finite-element formulation in FEAP, we have demonstrated the capability and efficiency of the proposed DFD model with three representative examples. For each of these examples, we have observed very good agreement between the discrete model and the previously developed continuous model [16,34]. The results indicate that the capability of the discrete model is equivalent to that of the continuous model when excluding fibres under compression but with a substantial computational speed-up.

For demonstration purposes, we discretized the unit sphere simply with spherical triangles. Future studies on the optimization of the discretization scheme may include non-uniform or locally refined discretizations. In addition, we have tested the proposed model with a rotationally symmetric fibre dispersion. It is straightforward to implement other types of fibre arrangements such as the recently observed non-symmetric fibre dispersion [11]. The only part that needs to be changed is the computation of the discrete fibre densities over each spherical triangle via numerical integration (algorithm 1). The remaining part of the implementation is the same.

Compared with the micro-sphere-based method, the proposed DFD model has several advantages.

— In the DFD model, dispersed fibres are treated in a discrete manner, as they would occur naturally. For example, figure 7

*e*shows the result of a rectangular strip with 40 discrete fibre bundles under uniaxial stretch. This result would then be exact for such a real fibre arrangement. On the other hand in the micro-sphere-based method, the dispersed fibres are considered to be continuously distributed. But their contribution is approximated by a numerical integration scheme, in which case the method is only valid for the functions for which the scheme was originally proposed. By contrast, the DFD method is independent of such numerical integration schemes, and any mechanical and failure properties of the fibres can be defined locally for each elementary area. If the representative fibre direction**N**_{n}is under compression, then all the fibres within can be excluded.— The discretization scheme used in the DFD method can be locally refined in a straightforward way as, for example, in [45]. When fibres are concentrated in several particular areas on a unit hemisphere, then the corresponding spherical triangles can be further divided into smaller ones to account for such locally concentrated fibre dispersions. However, if a micro-sphere-based method is used, it may not be possible to add additional integration points locally within a sub-domain of the unit hemisphere.

— Within each elementary area (spherical triangle), the discrete fibre density

*ρ*_{n}is evaluated by using a large number of*ρ*(*Θ*,*Φ*) values, whereas in the micro-sphere-based method*ρ*(*Θ*,*Φ*) is evaluated at a single integration point. Thus, the DFD model uses more*ρ*(*Θ*,*Φ*) information.— In the DFD model, the discrete fibre densities are computed by using a simple numerical integration scheme over the spherical triangles [38], and this integration is easy to perform because the variation of

*ρ*(*Θ*,*Φ*) over an elementary area is usually small. However, in the micro-sphere-based method the computation of the integration points over the unit sphere, and the associated weights of a particular integration scheme could be very complex or even impossible to perform by users.

A future comparison study of the two methods with locally concentrated fibres, as documented in [3], and various loading conditions will provide more insight into their relative performance.

As shown in the example of §3.3, the DFD solution approaches that of the CFD model as the discretization number *m* increases. With increasing *m*, the solution of the DFD model should theoretically approach the result of the CFD model for any problem. However, for larger *m*, the computational cost is higher. Thus, one aim when using the DFD model is to seek a balance between accuracy and computational efficiency. Certainly, we believe that the potential of this novel discrete model is far beyond its capability for excluding the compressed fibres that was demonstrated in the representative examples. The DFD model enables us to ‘fine tune’ the mechanical behaviour and failure properties of any fibre orientation at any point within the tissue. Future studies on the modelling of fibre recruitment and anisotropic fibre damage with the DFD model will demonstrate the further potential of the approach presented here.

## Authors' contributions

K.L. and G.A.H. designed the project and conceived the study; K.L. implemented the model in FEAP and conducted the numerical analyses. All authors interpreted the results. All authors contributed intellectually to the work presented, edited and approved the final version.

## Competing interests

We declare we have no competing interests.

## Appendix A. Mathematica code for discretization of a unit sphere

We have written a Mathematica code that was used for the discretization of the unit sphere into a finite number of spherical triangles, and we have tested the code in Mathematica v. 11.0 on a Windows machine. Thereby the number of spherical triangles can be increased by increasing the number ‘2’ in the definition of `gridFaces` in the following code.

`$PreRead=(# /.s_String /; StringMatchQ[s, NumberString] &&`

`Precision@ToExpression@s == MachinePrecision:> s <> “‘17.” &);`

`<< PolyhedronOperations‘`

`Quiet@Needs[“VectorAnalysis‘”]`

`gridFaces = Cases[Normal@Geodesate[PolyhedronData[“Icosahedron”], 2,`

`{0, 0, 0}, 1.0], _Polygon, Infinity];`

`Mappg = CoordinateTransformData[“Cartesian“ -> “Spherical”, “Mapping”];`

`Mg [vr_] = CoordinatesFromCartesian[vr, Spherical];`

`Grid[Table[`

`Flatten[{Mg[Partition[gridFaces[[i, 1]], 3][[1]][[1]]][[2 ;; 3]],`

`Mg[Partition[gridFaces[[i, 1]], 3][[1]] [[2]]][[2 ;; 3]],`

`Mg[Partition[gridFaces[[i, 1]], 3][[1]] [[3]]][[2 ;; 3]],`

` Normalize[RegionCentroid[ gridFaces[[i]] ]]}], {i,`

`Length[gridFaces]}] , ItemStyle -> “Text” ];`

`NumberForm[N[%, 16], {17, 8}, NumberFormat -> (#1 &)]`

For each spherical triangle, this program outputs the spherical coordinates of the three vertices in radians and the three Cartesian coordinates of its centroid in one row. Note that similar algorithms are also available for Matlab; see, for example, [46].

- Received October 14, 2017.
- Accepted January 5, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.