0
Research Papers

# Finite Element Algorithm for Frictionless Contact of Porous Permeable Media Under Finite Deformation and Sliding

[+] Author and Article Information
Gerard A. Ateshian1

Department of Mechanical Engineering, Columbia University, New York, NY 10027

Steve Maas, Jeffrey A. Weiss

Department of Bioengineering, University of Utah, Salt Lake City, UT 84112

The dyadic products of second order tensors $A$ and $B$ are defined such that $(A⊗B):X=(B:X)A$, $(A⊗̱B):X=A⋅X⋅BT$, and $(A⊗¯̱B):X=(A⋅X⋅BT+B⋅XT⋅AT)/2$ for any second order tensor $X$(14).

A semipermeable boundary typically represents a flow restriction specific to the boundary surface, as might occur when a porous-permeable body is enveloped by a semipermeable membrane, whose transport properties are distinct from those of the bulk material.

1

Corresponding author.

J Biomech Eng 132(6), 061006 (Apr 22, 2010) (13 pages) doi:10.1115/1.4001034 History: Received September 16, 2009; Revised October 02, 2009; Posted January 18, 2010; Published April 22, 2010; Online April 22, 2010

## Abstract

This study formulates and implements a finite element contact algorithm for solid-fluid (biphasic) mixtures, accommodating both finite deformation and sliding. The finite element source code is made available to the general public. The algorithm uses a penalty method regularized with an augmented Lagrangian method to enforce the continuity of contact traction and normal component of fluid flux across the contact interface. The formulation addresses the need to automatically enforce free-draining conditions outside of the contact interface. The accuracy of the implementation is verified using contact problems, for which exact solutions are obtained by alternative analyses. Illustrations are also provided that demonstrate large deformations and sliding under configurations relevant to biomechanical applications such as articular contact. This study addresses an important computational need in the biomechanics of porous-permeable soft tissues. Placing the source code in the public domain provides a useful resource to the biomechanics community.

###### FIGURES IN THIS ARTICLE
<>
Copyright © 2010 by American Society of Mechanical Engineers
Your Session has timed out. Please sign back in to continue.

## Figures

Figure 1

Confined compression stress-relaxation analysis in plane strain. The prescribed displacement uy on the top surface has a ramp-and-hold profile, with a ramp rate of −10−4 mm/s and a final displacement of −0.5 mm. (a) The contact model consists of two slabs constrained as shown (slab width=12 mm, height=2 mm). The top slab has 5×12 uniformly spaced elements; the bottom slab has 3×28 elements, with a mesh bias in the vertical direction to accommodate the boundary layer anticipated at the free-draining bottom surface. (b) The no-contact model consists of a single slab (width=12 mm, height=4 mm) with 3×40 elements and a mesh bias in the vertical direction in the bottom half, replicating the mesh distribution of the contact model in (a) along the vertical direction.

Figure 2

Transient response of the confined compression analysis of Fig. 1, at all the nodes of the contact surfaces in the contact model, and the corresponding nodes (same y coordinate) of the no-contact model: (a) nodal fluid pressure p; (b) nodal displacement uy. Symbols are displayed for all the nodes on each contact surface, but since the results are nearly identical on each surface, the graph appears to show only one symbol at each time point.

Figure 3

Unconfined compression stress-relaxation analysis in plane strain. The prescribed displacement uy on the top surface has a ramp-and-hold profile, with a ramp rate of −0.4 mm/s and a final displacement of −0.4 mm. The models are symmetric about the y-z plane. (a) The contact model consists of two slabs constrained as shown (slab width=3 mm, height=0.5 mm). The top slab has 1×41 elements with a mesh bias in the horizontal direction to accommodate the boundary layer anticipated at the free-draining right surface; the bottom slab has 1×40 elements, with a similar mesh bias. The deformed mesh is also shown at the end of the prescribed ramp displacement (t=1 s) and after the response has nearly reached equilibrium (t=105 s). (b) The no-contact model consists of a single slab (width=3 mm, height=2 mm) with 2×41 elements and a similar mesh bias in the horizontal direction, replicating the mesh distribution of the top slab of the contact model in (a) along the horizontal direction.

Figure 4

Spatiotemporal response of the unconfined compression analysis of Fig. 3, across the nodes of the contact surfaces in the contact model, and the corresponding nodes (same y coordinate) of the no-contact model: (a) nodal fluid pressure p; (b) nodal extrapolation of element effective stress component Tyye. The spatial distribution for these variables is shown at three select time points.

Figure 5

Normal contact of saddle-shaped layers. The layers are identical but rotated by 90 deg about the axis passing through the center of, and normal to, the contact surfaces (the y-axis). The contact surfaces were generated from the inner rim of the surface of a torus; their principal radii of curvature at the center point are −3.5 mm and +2.5 mm; their rim projects onto the x-z plane as a circle with a diameter of 3 mm. Each layer has a thickness of 0.4 mm and is supported on a rigid substrate; the bottom layer’s substrate is stationary and the top layer’s substrate has a prescribed ramp displacement uy at a ramp rate of −0.3 mm/s, with a final value of −0.3 mm. The model has a total of 21,420 nodes and 18,784 elements. Each layer’s mesh has a dual bias along the layer thickness (20 elements). The mesh is shown at three select time points.

Figure 6

Fluid pressure distribution for the contact model of Fig. 5, at t=0.5 s. Two sectioned views are presented, using the x-y plane (top) and the y-z plane.

Figure 7

Finite deformation and sliding contact of a semicylindrical slab (radius=2 mm) and a rectangular slab (width=12 mm, height=2 mm) in a plane strain analysis. The cylinder is initially prescribed a displacement profile of ux=2 mm and uy=(−0.375 mm/s)×t for 0≤t≤2 s, followed by sliding with a profile of ux=2+(0.5 mm/s)×t for 2 s≤t≤10 s. The cylinder has 20×100 elements along the radius and circumference, with a biased mesh along the radius to produce a finer mesh at the contact surface. The rectangular slab has 20×100 elements, with a dual bias mesh along the height. The deformed mesh and fluid pressure distribution are shown at two select time points.

Figure 8

A close-up of the model of Fig. 7 at t=2 s, using discrete color contours, demonstrates clearly the continuity of fluid pressure across the contact interface (emphasized with a white line trace)

Figure 9

Steady-state response of the contact traction tn and contact interface fluid pressure p for sliding of a rigid impermeable cylinder over a rectangular slab of biphasic material in a plane strain analysis. The slab (width=40 mm, height=1 mm) is bonded to a rigid impermeable substrate. The model (not shown) has 20×1024 elements, with a dual bias mesh along the height. The applied load intensity is 1 N/mm. Results are shown for three different Peclet numbers (representing the ratio of sliding velocity to characteristic diffusive velocity of interstitial fluid flow): (a) Pe=10−2; (b) Pe=100; (c) Pe=102.

## Errata

Some tools below are only available to our subscribers or users with an online account.

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

• TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
• EMAIL: asmedigitalcollection@asme.org
Sign In