A Study of preconditioned Krylov subspace methods with reordering for linear systems from a biphasic v–p finite element formulation

A study was conducted on combinations of preconditioned iterative methods with matrix reordering to solve the linear systems arising from a biphasic velocity–pressure (v–p) finite element formulation used to simulate soft hydrated tissues in the human musculoskeletal system. Krylov subspace methods were tested due to the symmetric indefiniteness of our systems, specifically the generalized minimal residual (GMRES), transpose-free quasi-minimal residual (TFQMR), and biconjugate gradient stabilized (BiCGSTAB) methods. Standard graph reordering techniques were used with incomplete LU (ILU) preconditioning. Performance of the methods was compared on the basis of convergence rate, computing time, and memory requirements. Our results indicate that performance is affected more significantly by the choice of reordering scheme than by the choice of Krylov method. Overall, BiCGSTAB with one-way dissection (OWD) reordering performed best for a test problem representative of a physiological tissue layer. The preferred methods were then used to simulate the contact of the humeral head and glenoid tissue layers in glenohumeral joint of the shoulder, using a penetration-based method to approximate contact. The distribution of pressure and stress fields within the tissues shows significant through-thickness effects and demonstrates the importance of simulating soft hydrated tissues with a biphasic model.     

A Lagrange multiplier mixed finite element formulation for three-dimensional contact of biphasic tissues

A three-dimensional (3D) contact finite element formulation has been developed for biological soft tissue-to-tissue contact analysis. The linear biphasic theory of Mow, Holmes, and Lai (1984, J. Biomech., 17(5), pp. 377-394) based on continuum mixture theory, is adopted to describe the hydrated soft tissue as a continuum of solid and fluid phases. Four contact continuity conditions derived for biphasic mixtures by Hou et al. (1989, ASME J. Biomech. Eng., 111(1), pp. 78-87) are introduced on the assumed contact surface, and a weighted residual method has been used to derive a mixed velocity-pressure finite element contact formulation. The Lagrange multiplier method is used to enforce two of the four contact continuity conditions, while the other two conditions are introduced directly into the weighted residual statement. Alternate formulations are possible, which differ in the choice of continuity conditions that are enforced with Lagrange multipliers. Primary attention is focused on a formulation that enforces the normal solid traction and relative fluid flow continuity conditions on the contact surface using Lagrange multipliers. An alternate approach, in which the multipliers enforce normal solid traction and pressure continuity conditions, is also discussed. The contact nonlinearity is treated with an iterative algorithm, where the assumed area is either extended or reduced based on the validity of the solution relative to contact conditions. The resulting first-order system of equations is solved in time using the generalized finite difference scheme. The formulation is validated by a series of increasingly complex canonical problems, including the confined and unconfined compression, the Hertz contact problem, and two biphasic indentation tests. As a clinical demonstration of the capability of the contact analysis, the gleno-humeral joint contact of human shoulders is analyzed using an idealized 3D geometry. In the joint, both glenoid and humeral head cartilage experience maximum tensile and compressive stresses are at the cartilage-bone interface, away from the center of the contact area.     

A penetration-based finite element method for hyperelastic 3D biphasic tissues in contact. Part II: finite element simulations

The penetration method allows for the efficient finite element simulation of contact between soft hydrated biphasic tissues in diarthrodial joints. Efficiency of the method is achieved by separating the intrinsically nonlinear contact problem into a pair of linked biphasic finite element analyses, in which an approximate, spatially and temporally varying contact traction is applied to each of the contacting tissues. In Part I of this study, we extended the penetration method to contact involving nonlinear biphasic tissue layers, and demonstrated how to derive the approximate contact traction boundary conditions. The traction derivation involves time and space dependent natural boundary conditions, and requires special numerical treatment. This paper (Part II) describes how we obtain an efficient nonlinear finite element procedure to solve for the biphasic response of the individual contacting layers. In particular, alternate linearization of the nonlinear weak form, as well as both velocity-pressure, v-p, and displacement-pressure, u-p, mixed formulations are considered. We conclude that the u-p approach, with linearization of both the material law and the deformation gradients, performs best for the problem at hand. The nonlinear biphasic contact solution will be demonstrated for the motion of the glenohumeral joint of the human shoulder joint.     

Hydrostatic pressurization and depletion of trapped lubricant pool during creep contact of a rippled indenter against a biphasic articular cartilage layer

This study presents an analysis of the contact of a rippled rigid impermeable indenter against a cartilage layer, which represents a first simulation of the contact of rough cartilage surfaces with lubricant entrapment. Cartilage was modeled with the biphasic theory for hydrated soft tissues, to account for fluid flow into or out of the lubricant pool. The findings of this study demonstrate that under contact creep, the trapped lubricant pool gets depleted within a time period on the order of seconds or minutes as a result of lubricant flow into the articular cartilage. Prior to depletion, hydrostatic fluid load support across the contact interface may be enhanced by the presence of the trapped lubricant pool, depending on the initial geometry of the lubricant pool. According to friction models based on the biphasic nature of the tissue, this enhancement in fluid load support produces a smaller minimum friction coefficient than would otherwise be predicted without a lubricant pool. The results of this study support the hypothesis that trapped lubricant decreases the initial friction coefficient following load application, independently of squeeze-film lubrication effects.     

Effects of friction on the unconfined compressive response of articular cartilage: a finite element analysis

A finite element analysis is used to study a previously unresolved issue of the effects of platen-specimen friction on the response of the unconfined compression test; effects of platen permeability are also determined. The finite element formulation is based on the linear KLM biphasic model for articular cartilage and other hydrated soft tissues. A Galerkin weighted residual method is applied to both the solid phase and the fluid phase, and the continuity equation for the intrinsically incompressible binary mixture is introduced via a penalty method. The solid phase displacements and fluid phase velocities are interpolated for each element in terms of unknown nodal values, producing a system of first order differential equations which are solved using a standard numerical finite difference technique. An axisymmetric element of quadrilateral cross-section is developed and applied to the mechanical test problem of a cylindrical specimen of soft tissue in unconfined compression. These studies show that interfacial friction plays a major role in the unconfined compression response of articular cartilage specimens with small thickness to diameter ratios.     

A finite element analysis of the indentation stress-relaxation response of linear biphasic articular cartilage.

The indentation problem of a thin layer of hydrated soft tissue such as cartilage or meniscus by a circular plane-ended indenter is investigated. The tissue is represented by a biphasic continuum model consisting of a solid phase (collagen and proteoglycan) and a fluid phase (interstitial water). A finite element formulation of the linear biphasic continuum equations is used to solve an axisymmetric approximation of the indentation problem. We consider stress-relaxation problems for which analytic solution is intractable; where the indenter is impermeable (solid) and/or when the interface between the indenter and tissue is perfectly adhesive. Thicknesses corresponding to a thin and thick specimen are considered to examine the effects of tissue thickness. The different flow, pressure, stress and strain fields which are predicted within the tissue, over time periods typically used in the mechanical testing of soft tissues, will be presented. Results are compared with the case of a porous free-draining indenter with a perfectly lubricated tissue-indenter interface, for which an analytic solution is available, to show the effects of friction at the tissue-indenter interface, and the effects of an impermeable indenter. While these effects are present for both thin and thick tissues, they are shown to be more significant for the thin tissue. We also examine the effects of the stiffness of the subchondral bone on the response of the soft tissue and demonstrate that the subchondral bone substrate can be modeled as a rigid, impermeable boundary. The effects of a curved tissue-subchondral bone interface, and the early time response are also studied. For physiologically reasonable levels of curvature, we will show that the curved tissue-subchrondal bone interface has negligible influence on the tissue response away from the interface. In addition, the short-time stress-relaxation responses of the tissue (e.g., at times less than 1s) demonstrate the essential role of the fluid phase in supporting the load applied to the tissue, and by extrapolation to shorter times characteristics of normal joint motion, suggest the essential role of a biphasic model in representing soft tissue behavior in joint response.     

A Novel Method for Single Sample Multi-Axial Nanoindentation of Hydrated Heterogeneous Tissues Based on Testing Great White Shark Jaws

Nanomechanical testing methods that are suitable for a range of hydrated tissues are crucial for understanding biological systems. Nanoindentation of tissues can provide valuable insights into biology, tissue engineering and biomimetic design. However, testing hydrated biological samples still remains a significant challenge. Shark jaw cartilage is an ideal substrate for developing a method to test hydrated tissues because it is a unique heterogeneous composite of both mineralized (hard) and non-mineralized (soft) layers and possesses a jaw geometry that is challenging to test mechanically. The aim of this study is to develop a novel method for obtaining multidirectional nanomechanical properties for both layers of jaw cartilage from a single sample, taken from the great white shark (Carcharodon carcharias). A method for obtaining multidirectional data from a single sample is necessary for examining tissue mechanics in this shark because it is a protected species and hence samples may be difficult to obtain. Results show that this method maintains hydration of samples that would otherwise rapidly dehydrate. Our study is the first analysis of nanomechanical properties of great white shark jaw cartilage. Variation in nanomechanical properties were detected in different orthogonal directions for both layers of jaw cartilage in this species. The data further suggest that the mineralized layer of shark jaw cartilage is less stiff than previously posited. Our method allows multidirectional nanomechanical properties to be obtained from a single, small, hydrated heterogeneous sample. Our technique is therefore suitable for use when specimens are rare, valuable or limited in quantity, such as samples obtained from endangered species or pathological tissues. We also outline a method for tip-to-optic calibration that facilitates nanoindentation of soft biological tissues. Our technique may help address the critical need for a nanomechanical testing method that is applicable to a variety of hydrated biological materials whether soft or hard.     

An evaluation of three-dimensional diarthrodial joint contact using penetration data and the finite element method

We have developed an approximate method for simulating the three-dimensional contact of soft biphasic tissues in diarthrodial joints under physiological loading. Input to the method includes: (i) kinematic information describing an in vitro joint articulation, measured while the cartilage is deformed under physiological loads, (ii) geometric properties for the relaxed (undeformed) cartilage layers, obtained for the analyses in this study via stereophotogrammetry, and (iii) material parameters for the biphasic constitutive relations used to represent cartilage. Solid models of the relaxed tissue layers are assembled in physiological positions, resulting in a mathematical overlap of the cartilage layers. The overlap distribution is quantified and converted via the biphasic governing equations into applied traction boundary conditions for both the solid and fluid phases for each of the contacting layers. Linear, biphasic, three-dimensional, finite element analysis is performed using the contact boundary conditions derived for each of the contacting layers. The method is found to produce results consistent with the continuity requirements of biphasic contact. Comparison with results from independent, biphasic contact analyses of axisymmetric problems shows that the method slightly underestimates the contact area, leading to an overestimation of the total traction, but yields a good approximation to elastic stress and solid phase displacement.     

A transversely isotropic biphasic model for unconfined compression of growth plate and chondroepiphysis.

Using the biphasic theory for hydrated soft tissues (Mow et al., 1980) and a transversely isotropic elastic model for the solid matrix, an analytical solution is presented for the unconfined compression of cylindrical disks of growth plate tissues compressed between two rigid platens with a frictionless interface. The axisymmetric case where the plane of transverse isotropy is perpendicular to the cylindrical axis is studied, and the stress-relaxation response to imposed step and ramp displacements is solved. This solution is then used to analyze experimental data from unconfined compression stress-relaxation tests performed on specimens from bovine distal ulnar growth plate and chondroepiphysis to determine the biphasic material parameters. The transversely isotropic biphasic model provides an excellent agreement between theory and experimental results, better than was previously achieved with an isotropic model, and can explain the observed experimental behavior in unconfined compression of these tissues.     

Technical Note: Modelling Soft Tissue Using Biphasic Theory - A Word of Caution

In recent years the biphasic approach initiated by Mow and coworkers, has been very popular in modelling soft, hydrated, cartilage tissues as well as other soft tissues, such as the brain. This work points out that due to the inherent inability of biphasic models in their present form to account for stress-strain rate dependence resulting from the viscoelasticity of the solid phase, the applicability of these models is limited to the loading conditions producing large relative velocities of phases.     

Modelling Soft Tissue Using Biphasic Theory — A Word of Caution

In recent years the biphasic approach initiated by Mow and coworkers, has been very popular in modelling soft, hydrated, cartilage tissues as well as other soft tissues, such as the brain. This work points out that due to the inherent inability of biphasic models in their present form to account for stress-strain rate dependence resulting from the viscoelasticity of the solid phase, the applicability of these models is limited to the loading conditions producing large relative velocities of phases.     

Three-dimensional finite element simulations of the mechanical response of the fingertip to static and dynamic compressions

The analysis of the mechanics of the contact interactions of fingers/handle and the stress/strain distributions in the soft tissues in the fingertip is essential to optimize design of tools to reduce many occupation-related hand disorders. In the present study, a three-dimensional (3D) finite element (FE) model for the fingertip is proposed to simulate the nonlinear and time-dependent responses of a fingertip to static and dynamic loadings. The proposed FE model incorporates the essential anatomical structures of a finger: skin layers (outer and inner skins), subcutaneous tissue, bone and nail. The soft tissues (inner skin and subcutaneous tissue) are considered to be nonlinearly viscoelastic, while the hard tissues (outer skin, bone and nail) are considered to be linearly elastic. The proposed model has been used to simulate two loading scenarios: (a) the contact interactions between the fingertip and a flat surface and (b) the indentation of the fingerpad via a sharp wedge. For case (a), the predicted force/displacement relationships and time-dependent force responses are compared with the published experimental data; for case (b), the skin surface deflection profiles were predicted and compared with the published experimental observations. Furthermore, for both cases, the time-dependent stress/strain distributions within the tissues of the fingertip were calculated. The good agreement between the model predictions and the experimental observations indicates that the present model is capable of predicting realistic time-dependent force/displacement responses and stress/strain distributions in the soft tissues for dynamic loading conditions.     

Three-dimensional finite element simulations of the mechanical response of the fingertip to static and dynamic compressions

The analysis of the mechanics of the contact interactions of fingers/handle and the stress/strain distributions in the soft tissues in the fingertip is essential to optimize design of tools to reduce many occupation-related hand disorders. In the present study, a three-dimensional (3D) finite element (FE) model for the fingertip is proposed to simulate the nonlinear and time-dependent responses of a fingertip to static and dynamic loadings. The proposed FE model incorporates the essential anatomical structures of a finger: skin layers (outer and inner skins), subcutaneous tissue, bone and nail. The soft tissues (inner skin and subcutaneous tissue) are considered to be nonlinearly viscoelastic, while the hard tissues (outer skin, bone and nail) are considered to be linearly elastic. The proposed model has been used to simulate two loading scenarios: (a) the contact interactions between the fingertip and a flat surface and (b) the indentation of the fingerpad via a sharp wedge. For case (a), the predicted force/displacement relationships and time-dependent force responses are compared with the published experimental data; for case (b), the skin surface deflection profiles were predicted and compared with the published experimental observations. Furthermore, for both cases, the time-dependent stress/strain distributions within the tissues of the fingertip were calculated. The good agreement between the model predictions and the experimental observations indicates that the present model is capable of predicting realistic time-dependent force/displacement responses and stress/strain distributions in the soft tissues for dynamic loading conditions.     

Biphasic finite element modeling of hydrated soft tissue contact using an augmented Lagrangian method

A study of biphasic soft tissues contact is fundamental to understanding the biomechanical behavior of human diarthrodial joints. To date, biphasic-biphasic contact has been developed for idealized geometries and not been accessible for more general geometries. In this paper a finite element formulation is developed for contact of biphasic tissues. The augmented Lagrangian method is used to enforce the continuity of contact traction and fluid pressure across the contact interface, and the resulting method is implemented in the commercial software COMSOL Multiphysics. The accuracy of the implementation is verified using 2D axisymmetric problems, including indentation with a flat-ended indenter, indentation with spherical-ended indenter, and contact of glenohumeral cartilage layers. The biphasic finite element contact formulation and its implementation are shown to be robust and able to handle physiologically relevant problems.     

An augmented Lagrangian finite element formulation for 3D contact of biphasic tissues

Biphasic contact analysis is essential to obtain a complete understanding of soft tissue biomechanics, and the importance of physiological structure on the joint biomechanics has long been recognised; however, up to date, there are no successful developments of biphasic finite element contact analysis for three-dimensional (3D) geometries of physiological joints. The aim of this study was to develop a finite element formulation for biphasic contact of 3D physiological joints. The augmented Lagrangian method was used to enforce the continuity of contact traction and fluid pressure across the contact interface. The biphasic contact method was implemented in the commercial software COMSOL Multiphysics 4.2^® (COMSOL, Inc., Burlington, MA). The accuracy of the implementation was verified using 3D biphasic contact problems, including indentation with a flat-ended indenter and contact of glenohumeral cartilage layers. The ability of the method to model multibody biphasic contact of physiological joints was proved by a 3D knee model. The 3D biphasic finite element contact method developed in this study can be used to study the biphasic behaviours of the physiological joints.     

A biphasic model for micro-indentation of a hydrogel-based contact lens

The stiffness and hydraulic permeability of soft contact lenses may influence its clinical performance, e.g., on-eye movement, fitting, and wettability, and may be related to the occurrence of complications; e.g., lesions. It is therefore important to determine these properties in the design of comfortable contact lenses. Micro-indentation provides a nondestructive means of measuring mechanical properties of soft, hydrated contact lenses. However, certain geometrical and material considerations must be taken into account when analyzing output force-displacement (F-D) data. Rather than solely having a solid response, mechanical behavior of hydrogel contact lenses can be described as the coupled interaction between fluid transport through pores and solid matrix deformation. In addition, indentation of thin membranes ( approximately 100 microm) requires special consideration of boundary conditions at lens surfaces and at the indenter contact region. In this study, a biphasic finite element model was developed to simulate the micro-indentation of a hydrogel contact lens. The model accounts for a curved, thin hydrogel membrane supported on an impermeable mold. A time-varying boundary condition was implemented to model the contact interface between the impermeable spherical indenter and the lens. Parametric studies varying the indentation velocities and hydraulic permeability show F-D curves have a sensitive region outside of which the force response reaches asymptotic limits governed by either the solid matrix (slow indentation velocity, large permeability) or the fluid transport (high indentation velocity, low permeability). Using these results, biphasic properties (Young's modulus and hydraulic permeability) were estimated by fitting model results to F-D curves obtained at multiple indentation velocities (1.2 and 20 microm/s). Fitting to micro-indentation tests of Etafilcon A resulted in an estimated permeability range of 1.0 x 10(-15) to 5.0 x 10(-15) m(4)N s and Young's modulus range of 130 to 170 kPa.