A finite element implementation for biphasic contact of hydrated porous media under finite deformation and sliding.

The study of biphasic soft tissue contact is fundamental to understand the biomechanical behavior of human diarthrodial joints. However, to date, only few biphasic finite element contact analyses for three-dimensional physiological geometries under finite deformation have been developed. The objective of this article is to develop a hyperelastic biphasic contact implementation for finite deformation and sliding problem. An augmented Lagrangian method was used to enforce the continuity of contact traction and fluid pressure across the contact interface. The finite element implementation was based on a general purpose software, COMSOL Multiphysics. The accuracy of the implementation is verified using example problems, for which solutions are available by alternative analyses. The implementation was proven to be robust and able to handle finite deformation and sliding.

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.

Biphasic finite element contact analysis of the knee joint using an augmented Lagrangian method.

Biphasic contact analysis is essential to obtain a more complete understanding of soft tissue biomechanics; however, only a limited number of studies have addressed these types of problems. In this paper, a theoretically consistent biphasic finite element solution of the 2D axisymmetric human knee was developed, and an augmented Lagrangian method was used to enforce the biphasic continuity across the contact interface. The interaction between the fluid and solid matrices of the soft tissues of the knee joint, the stress and strain distributions within the meniscus, and the changes in stress and strain distributions in the articular cartilage of the femur and tibia after complete meniscectomy were investigated. It was found that (i) the fluid phase carries more than 60% of the load, which reinforces the need for the biphasic model for knee biomechanics; (ii) the inner third and outer two-thirds of the meniscus had different strain distributions; and (iii) the distributions of both maximum shear stress and maximum principal strain in articular cartilage changed after complete meniscectomy - with peak values increasing by over 350%.

An augmented Lagrangian method for sliding contact of soft tissue.

Despite the importance of sliding contact in diarthrodial joints, only a limited number of studies have addressed this type of problem, with the result that the mechanical behavior of articular cartilage in daily life remains poorly understood. In this paper, a finite element formulation is developed for the sliding contact of biphasic soft tissues. The augmented Lagrangian method is used to enforce the continuity of contact traction and fluid pressure across the contact interface. The resulting method is implemented in the commercial software COMSOL Multiphysics. The accuracy of the new implementation is verified using an example problem of sliding contact between a rigid, impermeable indenter and a cartilage layer for which analytical solutions have been obtained. The new implementation's capability to handle a complex loading regime is verified by modeling plowing tests of the temporomandibular joint (TMJ) disc.

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.

A linear, biphasic model incorporating a brinkman term to describe the mechanics of cell-seeded collagen hydrogels.

Protein-based hydrogels are commonly used as in vitro models of native tissues because they can mimic specific aspects of the three-dimensional extracellular matrix present in vivo. Bulk mechanical stimulation is often applied to these gels to determine the response of embedded cells to biomechanical factors such as stress and strain. This study develops and applies a linear, biphasic formulation of hydrogel mechanics that includes a Brinkman term to account for viscous effects. The model is used to predict fluid pressure, relative velocity, and estimated shear stress exerted on cells seeded within a cyclically strained collagen hydrogel with and without imposed cross flow. The model was validated using a confined compression creep test of a cardiac fibroblast-seeded collagen type I hydrogel, and the effect of the added Brinkman term was assessed. The model indicated that the effects of strain and interstitial fluid flow are strongly interdependent in the collagen hydrogel. Our results suggest that the contribution of the Brinkman term is greater in protein hydrogels than in native tissues, and that studies that apply cyclic strain to cell-seeded hydrogels should account for the induced interstitial fluid flow. This study, therefore, has relevance to the increasing number of studies that examine cellular responses to mechanical stresses using in vitro hydrogel models.

The impact of biomechanics in tissue engineering and regenerative medicine.

Biomechanical factors profoundly influence the processes of tissue growth, development, maintenance, degeneration, and repair. Regenerative strategies to restore damaged or diseased tissues in vivo and create living tissue replacements in vitro have recently begun to harness advances in understanding of how cells and tissues sense and adapt to their mechanical environment. It is clear that biomechanical considerations will be fundamental to the successful development of clinical therapies based on principles of tissue engineering and regenerative medicine for a broad range of musculoskeletal, cardiovascular, craniofacial, skin, urinary, and neural tissues. Biomechanical stimuli may in fact hold the key to producing regenerated tissues with high strength and endurance. However, many challenges remain, particularly for tissues that function within complex and demanding mechanical environments in vivo. This paper reviews the present role and potential impact of experimental and computational biomechanics in engineering functional tissues using several illustrative examples of past successes and future grand challenges.

A two-dimensional computational model of lymph transport across primary lymphatic valves.

This study utilizes a finite element model to characterize the transendothelial transport through overlapping endothelial cells in primary lymphatics during the uptake of interstitial fluid. The computational model is built upon the analytical model of these junctions created by Mendoza and Schmid-Schonbein (2003, "A Model for Mechanics of Primary Lymphatic Valves," J. Biomed. Eng., 125, pp. 407-414). The goal of the present study is to investigate how adding more sophisticated and physiologically representative biomechanics affects the model's prediction of fluid uptake. These changes include incorporating a porous domain to represent interstitial space, accounting for finite deformation of the deflecting endothelial cell, and utilizing an arbitrary Lagrangian-Eulerian algorithm to account for interacting and nonlinear mechanics of the junctions. First, the present model is compared with the analytical model in order to understand its effects on parameters such as cell deflection, pressure distribution, and velocity profile of the fluid entering the lumen. Without accounting for the porous nature of the interstitium, the computational model predicts greater cell deflection and consequently higher lymph velocities and flow rates than the analytical model. However, incorporating the porous domain attenuates the cell deflection and flow rate to values below that predicted by the analytical model for a given transmural pressure. Second, the present model incorporates recent experimental data for parameters such as lymph viscosity, transmural pressure measurements, and others to evaluate the ability of these junctions to act as unidirectional valves. The volume of flow through the valve is calculated to be 0.114 nL/microm per cycle for a transmural pressure varying between 8.0 mm Hg and -1.0 mm Hg at 0.4 Hz. Though experimental data for the absorption of lymph through these endothelial junctions are scarce, several measurements of lymph velocity and flow rates are cited to validate the present 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 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 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.

A penetration-based finite element method for hyperelastic 3D biphasic tissues in contact: Part 1--Derivation of contact boundary conditions.

In this study, we extend the penetration method, previously introduced to simulate contact of linear hydrated tissues in an efficient manner with the finite element method, to problems of nonlinear biphasic tissues in contact. This paper presents the derivation of contact boundary conditions for a biphasic tissue with hyperelastic solid phase using experimental kinematics data. Validation of the method for calculating these boundary conditions is demonstrated using a canonical biphasic contact problem. The method is then demonstrated on a shoulder joint model with contacting humerus and glenoid tissues. In both the canonical and shoulder examples, the resulting boundary conditions are found to satisfy the kinetic continuity requirements of biphasic contact. These boundary conditions represent input to a three-dimensional nonlinear biphasic finite element analysis; details of that finite element analysis will be presented in a manuscript to follow.

Biphasic finite element simulation of the TMJ disc from in vivo kinematic and geometric measurements.

Understanding the biomechanical nature of the degeneration of the temporomandibular joint requires a coupling between experimental measurements and numerical simulation. In this study, geometry measured from MRI, and motion obtained from a specially designed optoelectronic system are fed into a three-dimensional biphasic finite element analysis to generate the spatial and temporal mechanical response of the disc. This study demonstrates how this coupling can effectively predict the biomechanical response of the temporomandibular joint disc to physiological loading. For small jaw opening movements, asymmetries in the load of the disc are found, with especially higher shear stresses in the lateral portion.

Mixed and Penalty Finite Element Models for the Nonlinear Behavior of Biphasic Soft Tissues in Finite Deformation: Part II - Nonlinear Examples.

This two-part paper addresses finite element-based computational models for the three-dimensional (3-D) nonlinear analysis of soft hydrated tissues, such as articular cartilage in diarthrodial joints, under physiologically relevant loading conditions. A biphasic continuum description is used to represent the soft tissue as a two-phase mixture of incompressible inviscid fluid and a hyperelastic, transversely isotropic solid. Alternate mixed-penalty and velocity-pressure finite element formulations are used to solve the nonlinear biphasic governing equations, including the effects of strain-dependent permeability and a hyperelastic solid phase under finite deformation. The resulting first-order, nonlinear system of equations is discretized in time using an implicit finite difference scheme, and solved using the Newton-Raphson method. Details of the formulations were presented in Part I [1]. In Part II, the two formulations are used to develop two-dimensional (2-D) quadrilateral and triangular elements and three-dimensional (3-D) hexahedral and tetrahedral elements. Numerical examples, including those representative of soft tissue material testing and simple human joints, are used to validate the formulations and to illustrate their applications. A focus of this work is the comparison of the alternate formulations for nonlinear problems. While it is demonstrated that both formulations produce a range of converging elements, the velocity-pressure formulation is found to be more efficient computationally.

Mixed and Penalty Finite Element Models for the Nonlinear Behavior of Biphasic Soft Tissues in Finite Deformation: Part I - Alternate Formulations.

This paper addresses finite element-based computational models for the three-dimensional, (3-D) nonlinear analysis of soft hydrated tissues, such as the articular cartilage in diarthrodial joints, under physiologically relevant loading conditions. A biphasic continuum description is used to represent the soft tissue as a two-phase mixture of incompressible, inviscid fluid and a hyperelastic solid. Alternate mixed-penalty and velocity-pressure finite element formulations are used to solve the nonlinear biphasic governing equations, including the effects of a strain-dependent permeability and a hyperelastic solid phase under finite deformation. The resulting first-order nonlinear system of equations is discretized in time using an implicit finite difference scheme, and solved using the Newton-Raphson method. Using a discrete divergence operator, an equivalence is shown between the mixed-penalty method and a penalty method previously derived by Suh et al. [1]. In Part II [2], the mixed-penalty and velocity-pressure formulations are used to develop two-dimensional (2-D) quadrilateral and triangular elements and 3-D hexahedral and tetrahedral elements. Numerical examples, including those representative of soft tissue material testing and simple human joints, are used to validate the formulations and to illustrate their applications. A focus of this work is the comparison of alternate formulations for nonlinear problems. While it is demonstrated that both formulations produce a range of converging elements, the velocity-pressure formulation is found to be more efficient computationally.