- Open Access
Simulating the swelling and deformation behaviour in soft tissues using a convective thermal analogy
BioMedical Engineering OnLinevolume 1, Article number: 8 (2002)
It is generally accepted that cartilage adaptation and degeneration are mechanically mediated. Investigating the swelling behaviour of cartilage is important because the stress and strain state of cartilage is associated with the swelling and deformation behaviour. It is well accepted that the swelling of soft tissues is associated with mechanical, chemical, and electrical events.
The purpose of the present study was to implement the triphasic theory into a commercial finite element tool (ABAQUS) to solve practical problems in cartilage mechanics. Because of the mathematical identity between thermal and mass diffusion processes, the triphasic model was transferred into a convective thermal diffusion process in the commercial finite element software. The problem was solved using an iterative procedure.
The proposed approach was validated using the one-dimensional numerical solutions and the experimental results of confined compression of articular cartilage described in the literature. The time-history of the force response of a cartilage specimen in confined compression, which was subjected to swelling caused by a sudden change of saline concentration, was predicted using the proposed approach and compared with the published experimental data.
The advantage of the proposed thermal analogy technique over previous studies is that it accounts for the convective diffusion of ion concentrations and the Donnan osmotic pressure in the interstitial fluid.
The extracellular matrix (ECM) of articular cartilage is negatively charged and hydrated [1, 2]. Under normal physiological conditions, the fixed negative charge in ECM is electrically neutralized by mobile cations and anions contained in the synovial fluid. The difference in ion concentration between the ECM and synovial fluid introduces a fluid pressure difference between them, the so-called Donnan osmotic pressure . At equilibrium, the swelling pressure, which consists of the Donnan osmotic pressure and a charge-to-charge expansion stress within the tissue, is counteracted by structural elements in the solid matrix.
It is well accepted that the swelling of soft tissues is associated with mechanical, chemical, and electrical events. The triphasic theory , which was generalized by Gu et al.  to include multiple ion species, is the only realistic model currently available to describe the mechanics of articular cartilage swelling. Gu et al.  and Gu  used the triphasic theory to analyse the coupling between fluid and ion transport, streaming potential, and deformation of cartilage during compression tests. Bryant and McDonnell  applied the triphasic theory to analyse corneal swelling under steady-state conditions in which the effects of fluid flow were neglected. Recently, Hon et al.  developed a rigorous numerical algorithm of the triphasic model for two-dimensional problems based on radial basis functions. However, all numerical/analytical approaches mentioned above are problem-oriented, and cannot be used for solving realistic biological problems, which are usually nonlinear, three-dimensional, irregular and complex in geometry, and undergo large deformations.
In order to apply the triphasic theory to biological systems, or applied biomechanical problems, the triphasic constitutive equations must be implemented into a finite element model (FEM). Simon et al.  proposed a poroelastic finite element formulation that included transport and swelling effects in soft tissue. This model might be considered a special case of the triphasic model by Lai et al. . Simon et al.  developed a one-dimensional FEM that was used to demonstrate the properties of triphasic materials. However, it could not be used for practical problems in two- or three-dimensional space. Sun et al.  developed sophisticated FEM formulations that included the coupling between mechanical, electrical and chemical events in articular cartilage. However, the approach proposed by  still in a conceptual stage and needs substantial research and development until it can be used for solving realistic biological problems.
Because the process of diffusion is mathematically equivalent to the process of thermal transfer, the transport of ions in cartilage may be analysed using a thermal analogy. Myers et al.  simulated swelling of articular cartilage using a thermal analogous technique. Kaufmann et al.  analysed poroelastic problems with swelling deformation using a commercially available FEM software package (ABAQUS). In these last two studies, the diffusion process of ion concentrations was analysed by neglecting convection, i.e., by assuming that the velocity of the tissue and the diffusion of ions are decoupled. For most substances, the diffusion process is strongly coupled with the velocity of the media. Therefore, the coupling between the velocity of the tissue and the diffusion of ions may have non-negligible effects in physiological problems.
Understanding the mechanisms of swelling and deformation is important when investigating the mechanics and adaptative behaviour of articular cartilage. The triphasic theory [4, 5] can be used to describe swelling and transport of ions in articular cartilage. However, a commercial finite element tool based on the triphasic theory, which bioengineers can use to analyse swelling and deformation behaviour for real-life biophysical problems, does not exist. The purpose of this study was to implement the triphasic theory [4, 5] into a commercially available finite element software package (ABAQUS), and to validate the proposed approach using the one-dimensional numerical solutions by Simon et al. , and the experimental results of confined articular cartilage compression obtained by Eisenberg and Grodzinsky . In order to solve the problem using ABAQUS, the triphasic model was converted into a thermo-mechanical interaction problem with convective heat transfer. It was assumed that the volume occupied by ions is small and that the effect of electric current on the flow of ions can be ignored.
where I is the unit tensor; σs, σf, and σtrepresent the stresses in the solid phase, in the interstitial fluid, and in the total tissue; λ and μ are the first and the second Lamé constants; p and ε are fluid pressure and strain of the solid phase, respectively. The stress in the solid phase (σs) is decomposed here into three parts: a component associated with the fluid pressure (-Φs p I), a chemical expansion stress (T c I), and an elastic stress associated with the deformation of the matrix. Φsand Φfare the solid and fluid volume fractions, respectively. Assuming that the volumetric fraction of ions, Φ- and Φ+, is negligible, as suggested by Gu et al.  and Sun et al. , one can write the equation of the mixture in the form:
Φs+ Φf+ Φ- + Φ+ ≈ Φs+ Φf= 1 (2)
The ECM is negatively charged, with a fixed charge density of cF. Everywhere within the tissue under physiological conditions, the electroneutrality condition is maintained by
c+ = c- + cF (3)
where c+ and c- are concentrations of mobile cations and anions, respectively. c+, c-, and cFare measured in moles per unit fluid volume.
The conservation of momentum is satisfied everywhere in the tissue
where υsand υfare the velocities of solid and fluid, k is the hydraulic permeability (open circuit). The continuity condition holds everywhere in the tissue, such that
.(Φfυf+ Φsυs+ Φ- υ- + Φ+ υ+) = 0 (5)
Darcy's law is obtained from the second equation of the conservation of momentum (4).
W i = Φf(υf- υs) = -k p (6)
where W i is the effective fluid flux.
Fluid pressure p is composed of hydrostatic pressure, p 0; and osmotic pressure , Δπ:
p = p 0 - Δπ (7)
The osmotic pressure is associated with the chemical potential and is expressed explicitly
π = ξRT [c + A 0 c2 + A 1 c3 + ··· ];
Δπ = π - π0 (8)
where π0 is the initial state of the system; ξ is the osmotic coefficient describing the character of the ion transport through the semi-permeable membrane (tissue). ξ approaches 1 for an ideal membrane, i.e., no ions can pass through the tissue. R and T are the gas constant and absolute temperature, respectively. c is the effective NaCl concentration within the tissue which is related to the concentration of c- and c+ by :
In a first order approximation, Eq. (8) takes the form of van't Hoff's equation:
π = ξR·T·c;
Δπ = π - π0. (10)
The dependence of the chemical expansion stress, T c , on the ion concentration has not been quantified based on a strict theory. Eisenberg and Grodzinsky  proposed a chemical expansion stress as a function of c
where β0 and c β are material parameters.
Another form of T c as a function of c- and cFwas proposed by Lai et al. 
where a 0 and κ are material parameters; γ and γ* are the activity coefficients of NaCl in the tissue and the external solution, respectively.
Conservation of mass requires
where ρα is the the apparent density for component α; the superscripts α = s, f, +, - represent the solid and fluid phase, and cations and anions, respectively; t is time.
Using the relations  ρα = cα Mα (α = +, -) with M+ and M- being the atomic weight of Na+ and Cl-, respectively, Eq.(13) can be written for electrolytes as
The diffusion of the electrolytes can be described using Fick's law
Dα cα = cα (υα - υs) α = -, + (15)
where Dα is the diffusion coefficient for component α.
Fick's diffusion law (15), combined with mass conservation (Eq.(14)), results in the governing equation of the convective diffusion process
where the effects of the volume strain rate of the solid phase on the diffusion process have been neglected. The diffusion problem (Eq. (16)) can be solved using a convective, thermal analogy, i.e., cα → θ with θ denoting temperature.
The equations of motion for fluid and electrolytes (c+ and c-) in the triphasic model of Lai et al.  [their Eqs. (30–32)] are reduced to Darcy's law (Eq. (6)) and Fick's law (Eq. (15)) when the effects of the electrical potential, the friction between ions and fluid, the friction between ions and anions, and the effects of the volumetric deformation rate of the solid phase on the fluid phase are neglected.
For a practical problem with three components, only c- needs to be solved, because c+ can be obtained from Eq. (3). In order to solve the triphasic problem, the triphasic constitutive relation (1) was converted into a biphasic equation  with thermal expansion,
ε = εe+ εc= εe+ K Δc I (17)
where ε is the total tissue strain, εeis the elastic strain associated with mechanical loading, and εcis the strain associated with the chemical expansion stress ("thermal strain"); K = K(c) is a "thermal expansion coefficient", defined as
where E is the Young's modulus of the tissue.
Eq. (17) is in a form that can be implemented into a general, commercial software package. A practical problem can now be solved using an iterative technique as follows:
The initial c distribution within the tissue is obtained by solving the diffusion equation (16) using υs= 0. The problem is equivalent to the diffusion of c through the tissue in a steady-state. The boundary condition, i.e., the distribution of c on the surface of the tissue, is known.
The triphasic problem [Eqs. (1)] is solved using the distribution of c in the tissue obtained in the previous step. The problem is solved using a thermal analogous procedure, as discussed previously. From this step, the distributions of the stress/strain and the displacement rate of the solid phase (υs) were obtained.
The diffusion equation (16) is solved again using the updated υsobtained in the previous step. This is realized using a convective thermal analogous procedure, as discussed previously.
The iterations [steps (b) and (c)] are continued until the problem has converged. The distributions of the stress/strain and the displacement rate of the solid phase (υs) were extracted from the simulation results via a customized subroutine, and the iterations were run manually. In the present simulations, the effective salt concentration in the bath is assumed to be constant.
In order to validate the proposed FEM approach, our results were compared to those obtained by Simon et al.  and Eisenberg and Grodzinsky . All simulations were performed using a confined compression configuration (Fig. 1). The tissue specimens were compressed in a confining chamber with a rigid, porous platen connected to a load cell. The concentration of NaCl in the bathing solution, (c(t)), and the compression displacement of the sample, (δ(t)), were prescribed, while the compression force (F(t)) was predicted in all simulations. ABAQUS (version 5.8) was used in all simulations . In the analysis, thermal expansion of the fluid was neglected. All numerical simulations were carried out using 8-node, biquadratic, axi-symmetric elements. The distribution of NaCl concentration within the bathing solution was considered to be uniform throughout the tests.
Numerical Test A was designed to reproduce the one-dimensional, finite element results by Simon et al. . The cartilage specimen had a thickness of 5.0 mm and a diameter of 30.0 mm. The specimen was compressed using a prescribed displacement and, at the same time, the ion concentration at the surface of the specimen was changed. The c- concentration was chosen as the reference parameter in the simulations. At t = 0, the specimen was in a strainless equilibrium condition and had an original, uniform c- concentration of 0.001 M. The specimen was compressed via a prescribed step function (Fig. 2a); the c- concentration on the surface of the specimen was increased suddenly and then kept constant (Fig. 2b) for the remainder of the test. All material parameters (D, k, H A , and v) were assumed to be constant (Tables 1, 2, 3).
Numerical Tests B and C were designed to reproduce the confined compression experiments by Eisenberg and Grodzinsky . The articular cartilage specimens had a diameter of 6.40 mm, and had an original thickness of 0.540 mm and 0.565 mm for Tests B and C, respectively. The specimens were compressed initially by an amount of δ0 = 0.0560 mm and 0.0185 mm for Tests B and C, respectively. The initial NaCl concentrations in the bathing solution were c* = 0.10 and 0.005 M for Tests B and C, respectively. The specimen was kept in this state for 2000 seconds until the predicted force reached a steady-state (F 0 = 2.31 N and 2.28 N for Tests B and C, respectively). Once the steady-state force had been reached, the NaCl concentration was changed suddenly from the initial value (c*) to 0.15 M and 0.10 M, for tests B and C, respectively, while the compression of the specimen was kept constant. The predicted and the measured  forces were normalized to the steady-state force F 0, and were plotted as a function of time. The aggregate modulus, H A , was assumed to be linearly dependent on c. Other parameters (D, k, and v) were assumed to be constant (Tables 1, 2, 3) .
Van't Hoff's equation [Eq. (10)] was used in all simulations to evaluate the osmotic pressure.
The numerical results of Test A are illustrated in Figs. 2 and 3. The predicted displacement at the center of the specimen and the c- concentration at different locations across the thickness of the specimen as a function of time are shown in Fig. 2a and 2b, respectively. c- concentration distributions across the thickness of the articular cartilage specimen as a function of time are shown in Fig. 3. In order to compare our three-dimensional (axi-symmetrical) simulation to the one-dimensional FEM simulation by Simon et al. , only the results on the central line of the specimen were used in these plots. The numerical results of Simon et al.  are also shown in these figures.
The numerical results of Tests B are shown in Fig. 4. The time-histories of the c- concentration prescribed on the top of the specimen (z/h = 1), together with those predicted at a depth of z/h = 1/8 and 7/8 on the central line of the specimen, are shown in Fig. 4a. The prescribed displacement at the top of the specimen as a function of time is depicted in Fig. 4b. The predicted specimen reaction force as a function of time, in response to the variations in salt concentration and displacement, are shown in Fig. 4c. Approximately 2000 s after the specimen was exposed to 0.10 M NaCl in the bathing solution and was compressed to δ0 = 0.056 mm, an equilibrium state (F 0 = 2.31 N) was reached (Fig. 4c).
The predicted total reaction force on the specimen as a function of time, normalized to the steady-state force, F 0, were compared with the experimental measurements by Eisenberg and Grodzinsky  in Fig. 5a and 5b, for Tests B and C, respectively. t = 0 (Fig. 5) is the time when the specimens reached equilibrium (Fig. 4c).
The interaction between mechanical and chemical stresses and deformations during cartilage deforma-tion was demonstrated in our numerical simulations. For example, a sudden change in the deformation state of the tissue caused a change in ion concentration across the thickness of the cartilage specimen (Fig. 2a and 2b), while a change in NaCl concentration caused to a change in stress within the articular cartilage (Fig. 4a,4b,4c).
Discussion and Conclusion
Experimental animal models have shown that articular cartilage becomes thicker [17, 18], softer , and more permeable  in the initial stages of osteoarthritis (OA). It is assumed that the adaptation and degeneration of articular cartilage is associated with the stress-strain state of articular cartilage during normal, everyday loading. The stress-strain state of articular cartilage is associated with the swelling and deformation of the cartilage, which is influenced by the changing mechanical and chemical environment. A general FEM tool using the triphasic theory, which bioengineers can use to analyse practical swelling and deformation behaviour for physiologically relevant problems, does not exist. Because of the mathematical identity between thermal and mass diffusion processes , it is possible to transfer the triphasic model into a suitable form for implementation into commercial finite element software to simulate the swelling and deformation behaviour of articular cartilage. The main difference between the proposed approach and the triphasic model  is that the electrical potential and its influence on the transfer of electrolytes (i.e., NaCl concentration) was neglected in the current analysis. The advantage of the proposed thermal analogous technique over previous studies [12, 13] is, that it accounts for the convective diffusion of NaCl concentrations and the Donnan osmotic pressure in the interstitial fluid.
Our simulations (Test A) (Figs. 2b and 3) on the time-dependent distributions of concentrations on the central line of the specimen agree well with the one-dimensional finite element results obtained by Simon et al. . The displacement at the center of the specimen (z/h = 0.5) showed a sudden change (Fig. 2a) when the prescribed displacement at the specimen surface was changed using a step function (Fig. 2a). This behaviour was not predicted by Simon et al. (their Fig. 3) , who gave no explanation why they did not predict this behaviour which should occur when considering the articular cartilage as a continuous material.
In order to reach an initial compression force (F 0), the cartilage specimen must be compressed long enough for the force to reach a steady-state. It was found that the amount of initial compression to reach a given compression force in our numerical tests was smaller than the corresponding compression in the experiments by Eisenberg and Grodzinsky . This difference in the magnitude of compression between theory and experiment is likely associated with the fact that in real tests, there is space between the lateral surface of the cartilage specimen and the chamber wall; while in our numerical simulations, the lateral boundary of the specimen was assumed to be fixed ideally.
The proposed approach differs from the poroelastic finite element formulation by Simon et al.  in that the fluid pressure, which is the sum of the hydrostatic and Donnan osmotic pressure, is governed by Darcy's law; while in the model by Simon et al. , the effect of the Donnan osmotic pressure was included into the fluid pressure through a "generalized Darcy's" law [Eq. (12b) in ].
Compared to the general triphasic model by Lai et al.  and Gu et al. , the effects of the electrical potential, the friction between ions and fluid, the friction between cations and anions, and the effects of the volumetric deformation rate of the solid phase on the fluid phase were neglected in the proposed approach. The effects of friction between cations and anions and the volumetric deformation rate of the solid phase on the fluid phase are believed to be negligible [5, 11].
The major difference between the proposed approach and the rigorous formulation by Sun et al.'s  is that the friction between ions and fluid was neglected in our model formulation. However, the effects of the friction between ions and fluid are included in the diffusion coefficient (Dα) in our model. Our model represents a reasonable approximation of complex biological systems and the associated problems, since the friction between ions and fluid cannot be measured directly, while the diffusion coefficient, which includes the effects of ions/fluid friction, can be determined experimentally.
In the present simulations, a constant osmotic coefficient, hydraulic permeability, and diffusion coefficient was used. However, in reality, these coefficients are not constant; they change as a function of the deformation state of the matrix, the ion concentration, and the electrical potential [6, 7]. Deformation-and ion concentration-dependent coefficients (osmotic coefficient, hydraulic permeability, and diffusion coefficients) can be incorporated into the proposed approach.
In summary, we have developed an approach that converts the triphasic analysis into a problem of convective thermal diffusion coupled with a poroelastic thermal stress analysis. The proposed approach makes it possible for bioengineers to analyse triphasic physiological systems using commercially available FEM software.
Edwards J: Physical characteristics of articular cartilage. Proc Inst Mech Eng 1967, 181: 16–24.
Mow VC, Kuei SC, Lai WM, Armstrong CG: Biphasic creep and stress relaxation of articular cartilage in compression?: theory and experiments. J Biomech Eng 1980, 102: 73–84.
Tombs MP, Peacocke AR: The osmotic pressure of biological macromolecules Clarendon Press, Oxford 1974.
Lai WM, Hou JS, Mow VC: A triphasic theory for the swelling and deformation behaviours of articular cartilage. J Biomech Eng 1991, 113: 245–258.
Gu WY, Lai WM, Mow VC: A mixture theory for charged-hydrated soft tissues containing multi-electrolytes: passive transport and swelling behaviours. J Biomech Eng 1998, 120: 169–180.
Gu WY, Lai WM, Mow VC: Transport of fluid and ions through a porous-permeable charged-hydrated tissue, and streaming potential data on normal bovine articular cartilage. J Biomech 1993, 26: 709–723.
Gu WY: A study of mechano-electrokinetic properties of charged-hydrated-soft tissues and biomechanics of articular cartilage. Ph.D. Thesis, Columbia University 1994.
Bryant MR, McDonnell PJ: A triphasic analysis of corneal swelling and hydration control. ASME J Biomech Eng 1998, 120: 370–381.
Hon YC, Lu MW, Xue WM, Zhou X: Numerical algorithm for triphasic model of charged and hydrated soft tissues. Computational Mechanics 2002, 29: 1–15. 10.1007/s00466-002-0307-1
Simon BR, Liable JP, Pflaster D, Yuan Y, Krag MH: A poroelastic finite element formulation including transport and swelling in soft tissue structures. ASME J Biomech Eng 1996, 118: 1–9.
Sun DN, Gu WY, Guo XE, Lai WM, Mow VC: A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated biological soft tissues. Int J Num Meth Eng 1999, 45: 1375–1402.
Myers ER, Lai WM, Mow VC: A contimuum theory and an experiment for the ion-induced swelling behaviour of articular cartilage. ASME J Biomech Eng 1984, 106: 151–158.
Kaufmann MV, Simon BR, Meyer PJ, Baldwin AL: Abaqus finite element analysis of a porohyperelastic material with swelling by thermal analogy. Advances in Bioengineermg 1995, ASME BED-Vol. 31: 103–104.
Eisenberg SR, Grodzinsky AJ: The kinetics of chemically induced nonequilibrium swelling of articular cartilage and corneal stroma. ASME J Biomech Eng 1987, 109: 79–89.
Katchalsky A, Curran PF: Nonequilibrium Thermodynamics in Biophysics. Harvard University Press, Cambridge 1967.
Wu JZ, Herzog W, Epstein M: Evaluation of the finite element software ABAQUS for biomechanical modelling of biphasic tissues. J Biomech 1998, 31: 165–169. 10.1016/S0021-9290(97)00117-6
Kiviranta I, Jurvelin J, Tammi M, Saamanen AM, Helminen HJ: Weight-bearing controls glycosaminoglycan concentration and articular cartilage thickness in the knee joints of young beagle dogs. Arthritis Rheum 1987, 30: 801–809.
Herzog W, Adams ME, Matyas JR, Brooks JG: A preliminary study of hindlimb loading, morphology, and biochemistry of articular cartilage in the acl-deficient cat knee. Osteoarthritis and Cartilage 1993, 1: 243–251.
Jurvelin J, Kiviranta I, Tammi M, Helminen HJ: Softening of canine articular cartilage after immobilization of the knee joint. Clin Orthop 1986, 207: 246–252.
Sah RL-Y, Grodzinsky AJ, Plaas AHK, Sandy JD: Effects of static and dynamic compression on matrix metabolism in cartilage explants. In Articular Cartilage and Osteoarthritis (Edited by: Kuettner KE, Schleyebach R, Peyron JG, Hascall V). New York, Raven Press 1991, 373–383.
Welty JR, Wicks CE, Wilson RE: Fundamentals of momentum, heat and mass transfer Wiley, New York 1976.
This study was supported by The Alberta Heritage Foundation for Medical Research, NSERC of Canada, the Medical Research Council of Canada, and The Arthritis Society of Canada.