- Email: [email protected]
A nonlinear viscoelastic membrane model applied to the human cervix
A NONLINEAR VISCOELASTIC MEMBRANE MODEL APPLIED TO THE HUMAN CERVIX* D. A. RICEand T. Y. Biomedical
Engineering
YANG
Center. Purdue University. West Lafayette. IN 47907. U.S.A.
Abstract-An exponentially elastic and viscous membrane model is developed and used to represent the human cervix during the first stage of labor in the vertex presentation. An axisymmetric orthotropic membrane finite-element is derived and evaluated. A procedure for solving for large displacement of structures confined to curved loci is derived and extended to include nonlinearly elastic and linear viscous material properties. The cervix is modeled. and the results are compared quantitatively with a previous simpler model of the cervix and qualitatively with published clinical observations. The model suggests a resolution to the controvers! over the time course of cervical dilatation.
1. INTRODUCTION It has become increasingly recognized that the progress of labor, indicated by a plot of cervical diameter as a function of time. is an important tool in the management of labor. Analysis of this plot often permits early separation of those labors which require special attention from those which are normal. The approaches used have been limited to establishing average or ‘normal’ labor progression (e.g. Studd. 1973) or to fitting measured progression to empirical equations (e.g. Friedman and Kroll, 1969). Insufficient work has been done to relate the cause of cervical dilatation, the contractions of the uterus, directly to cervical dilatation. This paper develops a mechanical model which predicts time dependent deformation of the human cervix from intrauterine pressure caused by uterine contraction. The model is in the form of an axisymmetric. orthotropic membrane finite element. This element has applicability to membraneous organs of rotational symmetry. Although this development is directed towards the particular problem of cervical dilatation, it need not be limited to the cervix. It provides the detailed information on the distribution of stresses and displacements in the membrane. The membrine of the cervix is stretched in the process of dilatation so that it conforms to the fetal head as the fetus is forced through it. Therefore. a method is derived to solve for deflection of elastic structures confined by sliding restraints and extended to include curved restraints. An iterative-incremental approach is used for the nonlinearities of geometry and elasticity and for viscoelasticity. Thus this model formuiation provides a basis for a relatively accurate mechanical representation of the cervix. In the first attempt to relate intrauterine pressure to cervical dilatation, Rice ef al. (1975) developed a ring model of the cervix. Although the ring model is an extreme geometric simplification of the birth system, high correlation between model results and
* Received 2 August 1975.
measured data for relatively short labor segments was achieved through a statistical regression technique. The results presented here show that a thin membrane ring is equivalent to the ring model and. further, that the full membrane model of the cervix covering a substantial portion of the fetal head agrees qualitatively with clinical observations of dilatation over large time periods. The following basic assumptions underlie this development : (1) The cervix, fetus and uterus act as an isolated mechanical system which is approximately axisymmetric. (21 Intrauterine fluid pressure provides the force which dilates the passive cervix. (3) Friction between the head and cervix and the inertial forces are negligible. These assumptions were made and justified by Rice er al. in their previous work. In the text of this paper, the cervix refers to all parts of the uterus below the equator of the fetal head.
2. AXISYMMETRIC,
ORTHOTROPIC,
MEMBRANE
FINITE-ELEMENT ?.I
Description of the finite-element
Following the concept of finite-element method, an axisymmetric membrane can be partitioned into a finite number of conical frusta, each of which is an axisymmetric membrane finite-element. This finite-element is described geometrically in Fig. I. All the geometric variables are defined in the figure. The two nodal edges which bound the element are parallel circles. The element is composed of an orthotropic membrane whose principal axes are in the meridional and circumferential directions. The axisymmetric conditions of geometry and loading restrict the circumferential displacement to zero, so the displacement functions are independent of the circumferential coordinate and are described in terms of a set of two-dimensional Cartesian coordinates r and z only. For each nodal circle, two displacement
201
202
D. A. RICE and T. Y.
YANG
degrees-of-freedom u and L‘in the r and : directions, respectively, are assumed. This element thus has four degrees-of-freedom
element and is found by dividing extension by the length L
where L 1 is the symbol for a row matrix.
For an orthotropic membrane in plane stress. the stress-strain relations are
2.2 Stlfiess formulation In this axisymmetric membrane element, the meridional displacement is the only displacement that contributes to the strain energy. Although two displacement degrees-of-freedom are assumed for each no&l circle, they are the r and z components of the meridional displacement. Thus only the meridional displacement function w need be assumed. It is assumed to be of linear variation along the meridional direction x,
x/w, + WJW,,
w(x) = (1 -
(2)
the meridional
where v is the Poisson’s ratio; subscripts nr and c refer to the meridional and circumferential directions, respectively; and D is the membrane rigidity defined as D,
=
E --!!Land 1 - \‘*Vc
DC= A.
E
1 - v,v,,
(7)
Substituting equations (4 and 5) into (6). the stresses in terms of the nodal displacements are obtained
L
I
The strain energy, without shearing stress, is where w, and w2 are the meridional displacements at the nodal circles (1) and (2), respectively. If 4 is used to define the meridional direction of the element, the following relations can be established, u = iw; r = iiw; r = r, + i.x A = cosfj;
z = 2, + /lx;
The circumferential &=
(3)
strain is found as L-x
r
p = sin4.
L(-_) rL
,o,;.o
J id,
(4)
U = ;
sV
(u,E, + u,,,~,,,) d I’,
(9)
where the differential volume dV = 27&(x) dx. Substituting equations (4, 5 and 8) into equation (9), the strain energy expression is obtained in terms of the four nodal degrees-of-freedom. Following Castigliano’s theorem and performing partial differentiation of the strain energy with respect to each degreeof-freedom, the stiffness equations that relate the four nodal forces (f] to the four nodal displacements {q} are obtained.
where { ] is the symbol for a column matrix, and the meridional strain is assumed constant over the
{S} =
Cklid.
(10)
The content of the stiffness matrix [k] is given in Table 1. 2.3 Evaluation
“2
Fig.
I. Description
of
the
axisymmetric
. finite-element.
membrane
The above formulation is evaluated through numerical examples to indicate its efficiency. The first example chosen is an isotropic conical membrane shell. Its apex is closed and its open edge is constrained from meridional movement. It is subjected to uniform internal pressure. The exact solution for this example is available in the text by Kraus (1967). The values assumed for the example are: internal pressure = 1 gf/cm’; modulus of elasticity = loo0 gf/cm* ; thickness = 0.5 cm: meridional length = 10 cm; and Poisson’s ratio = 0.4. The displacement shapes obtained by using 1, 2 and 10 elements are shown in Fig. 2. The lO-element solution virtually coincides with the exact solution.
A nonlinear viscoelastic membrane model Table I. Stiffness matrix for the axisymmetric orthotropic membrane finite-element
---from Kraus -
IO rlrmrnts
Fig. 3. Displacement shapes ofa spherical membrane loaded under its own weight.
R
I
-+_+‘;” i.
i2L t py
I
:-l
I’.
-1
:r
+ S‘
Q = In (rdr,)
i. = cos fj p = sin4
s, = 27rhE,/(l - V,P,)
R = (r, + r2)/2
s, = ZnhE,/(l - \Iml’,)
An indication of the convergent characteristic of this element can be obtained by varying the number of elements in the analysis. The displacements at the open edge are found to be approx. - 35. - 5.5, - 1.85, -0.8, -0.4 and -0.3% off the exact solution for the idealization of 1, 2, 3. 5, 8 and 10 elements, respect-
ively. The similar convergent trend is also noted in the results for meridional membrane stresses. Details are given in the thesis by Rice (1974). The second example is an isotropic spherical membrane loaded under its own weight as shown in Fig. 3. The opening edge is restrained from meridional movement. The values used for the example are: density per unit area = 1 g/cm2: modulus of elasticity = loo0 gf/cm2; thickness = 0.05cm; radius = I cm; opening angle = 120”; and Poisson’s ratio = 0.3. The results are found to have the same convergent trend as that noted in the first example. The displacement shape for a lOelement solution is shown in Fig. 3. Except for the cetitral element, all the nodal displacements are within 1% of the displacements computed from the exact solution in the text by Kraus (1967). The large error in the central element is expected since the element is poor when its nodal radii are small. Because the cervix representation does
also exact soluti
Fig. 2. Displacement shapes of a membrane cone under uniform pressure.
204
D. A. RICEand T. Y.
not require a zero nodal radius, the present element is adequate.
YANG
or symbolically
1Qd = 2m x 1
3. NODAL CONSI-RAINTFORMULATION
This section formulates a method of sliding constraints to force the nodal points of the finite-element to remain on a specified curve, the locus. The method is first developed for linear loci and then used with an incremental-iterative procedure for curved loci. The derivation assumes two degrees of translational freedom for each node, but its generality is not limited by this assumption.
Gel 2tnxm
Similarly, the components can be expressed as
Id)
(15b)
mxl‘
of the reaction forces {s)
cos41
‘XI‘
Yl
sin 4,
X2
cos& sin Cp2
y2
<
>
3.1 Linear loci
=
After the assemblage, the overall force-displacement stiffness equations are in the form that IF)
=
2n x 1
{QI
CKI
2nx
2nx2n
(11)
1,
where n is the number of nodal points and the capital letters are used to distinguish the system equation (11) from the element equations (10). The force vector {F) is separated as the sum of the applied load vector {P} and the reaction force vector {R} . Equation (11) thus becomes
XIII
cosAl
’ XII ’
sin 4,
or symbolically IQ
=
2tn x 1
IPI = 2n x 1
CKcIK,l 2n x 2m
I ii I
2n x 2j
-
Q/
2mxl -
2j x 1
0
’
2j x 1
(13) where the subscripts ‘c’ and ‘f designate the constrained and the free nodes, respectively. A typical constrained node ‘i’, displaced by di along a linear locus, is shown in Fig. 4. The constrained condition requires that Ui = diCOSBj = dili ri = di sin ~9~= dipi Thus the constrained be expressed as
with i = 1, 2,..., m. displacement
(4 -
[PI = CKICel 2n x m
2n x 1 i
components
Cal
2m x m
co1
tn x 1
CKJI 2n x 2,j
2j x m
R,
QC 2mxl
ISI m x 1’
Substituting equations (15b and 16b) into equation (13) yields
113 = {p) + IN = CUQI. (12) Then nodes are numbered sequentially so that the first m nodes are those to be constrained and the other (n-m) = j nodes are those left unconstrained, or free. Equation (12) is partitioned to separate the free from the constrained nodes.
[cl 2tn x m
I(
or symbolically, 10 = CR1 2n x 1 2n x 2n
(4 2n x 1’
{s) mxl
{Q/I 2jx
1I (174 UW
where [K] is the modified stiffness matrix which accounts for the linear locus boundary condition. The constrained nodal displacements {QC} and the reaction forces {RC} can be obtained by substituting the appropriate section of the solution vector {A} into equations (15 and 16), respectively.
(14) can REACTION Si
Initial position
Fig. 4. A nodal point constrained on a linear locus.
105
A nonlinear viscoelastic membrane model LOAD INCREMENT %P
ONE INCREMENTAL
STEP ---I
Fig. 5. The iterative procedure for one linear-incremental step of a nodal point constrained on a curved locus. It must be noted that the above formulation is limited only to the case of small displacements without geometric nonlinearity, such as the examples given in Figs. 2 and 3. It is also noted that such model gives the distributions of displacements and stresses (in both the circumferential and meridional direction) along the meridional direction of the membrane such as those shown in Figs. 2 and 3. The ring model developed earlier by Rice et al. (1975) provides only the dilation and the hoop force at the reference circle. It gives neither the stresses nor the distributions of displacements and stresses in the membrane surface.
3.2 Curved loci The fetal head is of an ellipsoidal shape. The membrane finite-elements are constrained to slide along
the ellipsoidal surface with large displacements. The geometrical nonlinearity caused by large displacement can be treated by a linear-incremental procedure b! revising the geometry at the beginning of every increment. Details of this incremental procedure can be found in the paper by Yang (1973). Because the locus is curved, an iterative procedure is used at ever! incremental step. This iterative procedure is described below with the aid of Fig. 5. The system shown in Fig. 5 is a spring whose lower and upper nodes are free to slide along the horizontal and curved loci, respectively. A load increment iif is applied to the upper node and the displaced position at e is sought. This system is used to simulate the constrained condition for the finite-element systern. In Fig. 5. the point o is the initial position of a nodal point at the beginning of an incremental step. When the desired load-increment is applied. the node is displaced along the tangent at o to the point h’. a trial solution. Point b on the locus is found by taking a path. as a vertical line. from h’ to the locus. A position such as b on the locus found in this way is a target point. A new constrained sliding path is defined by line z. With the node at its initial position, the same load increment is applied. The new displaced position c’ and then the new target point c are found. The process is repeated until an error tolerance is met. The reaction force is specified to be perpendicular to the locus at the target point. When the displaced upper node lies on the locus, the reaction, the applied load. and the spring tension are in equilibrium and the target point is the solution. For a finite-element system, this procedure is applied to the assembled set. 1
AX I AL
FORCE
AXIAL
Fig. 6. LoadAisplacement
\
FORCE
(prrssurc
x 5**,
gf)
relation for a linearly elastic ring sliding along an ellipsoid.
206
D. A. RICE and T. Y. YANG
To evaluate this procedure. an example of a ring constrained to slide along the surface of an ellipsoid is chosen. The ring is loaded uniformly in the tangential direction to the ellipsoid as shown in Fig. 6. This example has been solved by Rice et al. by a simple ring model. Only a single finite-element is used in this evaluation. Since Rice et al. included the stiffness of the ring in the derivation in a form that no cross sectional area need be specified, the cross sectional area in this finite-element should be assumed to be as small as possible. In this finite-element, the meridional length is assumed as 0.07cm. The thickness is implicitly included in the membrane rigidity SCwhich is assumed as 2.12 gf/cm. The average initial radius is 2 cm. The load-displacement solution obtained by both the ring model and the finite-element model are shown in Fig. 6. Total agreement is seen. During the iteration in each increment, the path to the locus from a trial solution is along the radial direction of the ellipsoid. In the unstable region the solution is not obtainable by the method proposed here. The solution in that region is, however, not desired for practical purposes. 4.
NONLINEAR
MATERIAL
In treating the nonlinear stress-strain (a vs e) relationship, the concept of tangent modulus E, is often used. This modulus is defined in a stress-strain curve as the slope of the tangent at the point with the current strain. The use of tangent modulus allows a linear-incremental stiffness formulation while simultaneously satisfying the nonlinear stress-strain relation. This concept has been successfully used by Yang
and Wagner (1973) in treating the nonlinear materials. The tangent modulus is defined in an equation form as E, = E,(E) = da(e)@. (18) In the ring example by Rice et al.. the exponentially elastic material was assumed F = A~“-‘,,’
(19)
where F is the hoop tension; (r - ro) is the extension in radius; and A and B are material constants. Substituting equations (19) into (18), the tangent modulus becomes E, = abe@
for
E > l//I,
where the strain l = (r - r,)/r,; constant fl = Br,; and constant a = A divided by the cross sectional area. Based on the material property described by equation (20), the example described by Fig. 6 is again performed. The constants a and fl are assumed as 0.25 gf/cm’ and 4, respectively, so that they conform to those constants used in the ring example. The two solutions are compared in Fig. 7. Roth virtually coincide. In performing this example, the tangent modulus used in each iterative cycle is based on the average value of the strain at the beginning of the incremental step and the strain at the current target point. 5. CERVIX
MODEL
Since it is beyond the scope of this paper to characterize cervical tissue by formulating its strain energy
5
ring
-
model
O one finite
2
0
2
4
6 AXIAL
Fig. 7. Load-displacement
(20)
6
FORCE
elamrnt
IO
(gf
relation for an exponentially elastic
12
14
I
1 ring sliding
along
an ellipsoid.
A nonlinear
axial
dilatation
--_
ring
-
IO
I
I
20
30
40
AXIAL
FORCE
relation
+ (B,,&‘]’ ‘;
I
for both linearly
function. a modification of the exponential form proposed by Fung (1967) is proposed. The tissue is assumed to satisfy the following conditions: (1) The membrane is orthotropic with zero stress across its thickness. (2) There is no shear strain along the tissue axes. (3) The tissue is incompressible (v = 0.5). (4) The strain energy equals zero for the state of zero strain. The tangent modulus is assumed as (E,); = cc/?,exp[(&,.?
model
force
IO Fig. 8. Load-radial
viscoelastic membrane
i = c. m, (211
where r and /? are material constants; the subscripts c and nr refer to the circumferential and meridional directions, respectively. For the case of isotropic material in uniform biaxial strain, equation (21) reduces to the same form as that of equation (20). The cervix covers the fetal head for 85”” of the axial length below the equator and is partitioned into 10 finite elements. The head below the equator is an ellipsoid with the equatorial radius a = 5.0cm and axial radius h = 7.0cm. Above the equator the locus is assumed to extend indefinitely with a constant radius equaling the equatorial radius. The fetal head and body are thus represented by a circular cylinder with an ellipsoidal end. The cervix is assumed as an elastic membrane with a constant thickness h = 0.5 cm. The lowest node. node 1. has an initial radius of 2.0 cm. The last node. node 11. is located initially 0.5 cm below the equator. Load is applied in increments and the iterative procedure is used to solve for displacement. The material constants used for linear material are: E, = E, = 40 gf/cm* and 1’= 0.5. Those for nonlinear material are: J = 0.25 gf/cm2; 8, = fl,,, = 4; and v = 0.5. The maxi-
model finite-element
I
1
60
50
I
70
00
(gf 1 and nonlinearly
elastic membrane
cervix
mum load for the parameters specified is equivalent to an intrauterine pressure of 1gfjcm’. The results are shown in Fig. 8. A study by varying the number of elements shows that the results converged to the ones presented in Fig. 8 at the 7 element level.
6. COMPARISON
OF RING AND
MEMBRANE
Assuming that the proposed strain energy function of the membrane is representative of cervical tissue it is reasonable to expect that the load-ilisplacement relation of the model is similar to the pressure-dilatation relation of the biological system. Further vcrification of this is shown by Fig. 8 where the pressureradius solution for the ring described by Rice it ul.. is successfully superimposed on the IO element membrane solution. The parameter values of the ring. except for the multiplicative constant, are similar to the values found experimentally by Rice or ill. The! are: A = 0.215gf/cm*; B = 2cm-‘: and I’(,= Zcm.
7. VISCOELASTICITY
The creep present in cervical tissue during dilatation is handled by an extension of the incrementaliterative procedure. To represent this creep as a series linear element, the viscous strain .C is expressed as 1’51 where q is the viscous coefficient
and t
is
lhe time.
D. A. RICE and T. Y. YANG i-
and equatorial radius = Scm. initial node I radius = 2.0cm. In order to compare the qualitative model results with physical observation. a time varying axial load in the form that axial force = 1 + 10 sin2(0.2t)(gf)
0
20
60
60
40 TIME
100
is used to approximate the periodic nature of the intrauterine pressure. This load is applied to the same model as described in Fig. 9. The solution is shown in Fig. 10. The regression technique used with ring model by Rice et al. is applied to this solution. The multiple correlation coefficient is y > 0.99. The use of statistical technique here is justified because the convergence criterion of the membrane model allows some error in node position. As noted by Hildebrandt et al. (1969), the results of biaxial membrane strain are not predictable from uniaxial experiments, so the parameters of the ring model are not equivalent to the membrane parameters.
120
(s*condr)
Fig. 9. Radial dilation at constant loads due to the viscous effect. This viscous strain is obtained by using a trapezoidal
8. RESULTS
approximation, and is used to calculate the viscous displacements. A finite-element linear-incremental procedure has been used by Yang and Lianis (1974) to solve the geometrically nonlinear large displacement problems of viscoelastic beams and frames. The same concept and procedure are used here. Since the membrane is constrained to slide along a curved locus, the aforementioned iterative procedure, associated with the use of equation (17), is used in every incremental step. Two constant loads, 15 and 30% of the maximum loads shown in Fig. 8, are applied until creep is sufficient to allow displacement to proceed to the unstable region. The results are shown in Fig. 9. The following values were used: qC= Q,, = lo4 gf-sec/cn?; a = 0.25 gf/cm’; /3, = /?,,,= 4; ellipsoid axial radius = 7 cm;
DILATION
Friedman and Kroll (1969) suggest that the dilatation rate decelerates to near zero as dilatation approaches its maximum. Hendricks et al. (1970) suggest that the dilatation rate accelerates until full dilatation. Studd (1973) presents average dilatation curves which support neither of these interpretations, but as shown by Figs. 9 and 10, the membrane model supports the interpretation of Hendricks et al. The conflict between these two interpretations is resolved if it is assumed that for the larger dilatations the head has descended sufficiently to begin stretching structures other than the cervix. The model formulation limits the system to the cervix, uterus and fetus so such descent violates the basic assumptions, and the model is no longer applicable. If Hendrick’s patients
RADIUS
100
TIME
(26)
( seconds)
Fig. 10. Ten element simulatjon of tbe exponentially elastic and viscous membrane model.
209
A nonlinear viscoelastic membrane model such that the cervix is higher in the pelvis or the pelvis is larger than that for Friedman’s patients. the model formulation applies to full dilatation and predicts the accelerating curve. The position of the cervix. whether higher or lower in the pelvis. may be the determining factor for the shape of the dilatation time curve. Embrey and Siener (1965) make several observations concerning the nature of cervical dilatation. Only one of these observations. the passive nature of the cervix. is assumed by the model. Model results are compared to each observation below. 1. Small progress is made for each contraction, but dilatation increases during a contraction. The model shows an elastic stretch during a contraction and a small viscous strain during the course of a contraction (Fig. 10). This is due to the viscoelastic nature of cervical tissue. 3. A contraction causes dilatation in proportion to its strength. This statement implies a monotonic relationship between pressure and dilatation, not necessarily a linear relationship. The model shows a monotonic relationship (Fig. 8). 3. A larger total dilatation results in a larger elastic dilatation for a contraction of given strength. Figure 10 definitely confirms this. 4. There is a delay between the onset of a contraction and dilatation which lessens as dilatation increases. The model shows no such delay. 5. The drug oxytosin causes greater contractions and thus faster dilatation. Oxytosin’s action is well known. That greater pressure causes faster dilatation is shown by Fig. 9. The models predict all the observations except (4). The delay could be a_ measurement artifact or the model may be inadequate to describe the system in this respect. Data from Romero-Salinas et al. (1964) and data collected by Rice er al. show no such effect. were selected
9. DISCL’SSION
The model developed herein for membrane elastica problems is restricted only by the requirements of axisymmetry and quasi-static equilibrium. This model should. therefore. have applicability to many biomechanics problems involving a large class of hollow membrane organs. Because the development allows for the finite element nodes to be either free or confined by sliding restraints. these organs could be filled by solid objects. These organs could also be filled by fluids or a combination of solids and fluids if the stiffness formulation is modified to account for the fluid pressure. The geometric shape of the solid contained within the membrane must be axisymmetric although it need not be defined by a mathematically continuous function. The examples presented assumed an ellipsoid, but a tabulated function should do as well. The assumption in the theoretical development that zero friction occurs at the sliding restraints is not as
restrictive as it may seem. By appropriately specifying the angle of the reaction force to the locus (to be other than perpendicular) the effects of static, kinetic. or viscous friction can be simulated. The cervix model uses the membrane development to represent the anatomical structure of the cervix. Given a simulated intrauterine pressure. the model produces a numerical solution for dilatation which agrees well with published qualitative observations which encompass a large portion of labor time. Further, this solution shows agreement with an earlier one-dimensional ring model which is highly significant in explaining quantitative data. Since the parameters that must be assumed in the membrane model are considerably more than those in the ring model. direct quantitative comparison with experimental data is not performed. To find those parameters is beyond the scope of this paper. This cervix model thus completes the bridge between measured data, dilatation and pressure. and the morphology of the cervix during the first stage of labor. The major limitation of the model is that there is presently no way to evaluate its parameters (initial length, elastic coefficient. viscous coefficient. Poisson’s ratio for each direction for each element) noninvasively. Given reasonable assumptions for the parameter values, however, much information about the state of the cervix may be inferred: the distribution of stress and strain; the nodal positions; and the pressure distribution over the fetal head as a function of intrauterine pressure and time, 10. CONCLUDING
REMARliS
A membrane finite element model was developed which has general applicability to axisymmetric hollow organs having solid contents. This model was shown to be suitable for the problem of human cervical dilatation. Using an exponentially elastic and viscous membrane, the model displays many of the qualitative characteristics observed in the cervix. Further, the model correlates well with an earlier simple. but quantitative, ring model of the cervix. These facts lend credence to the soundness of the ring model and to the ability of the membrane model to represent the cervix during dilatation. The model suggests that ccrvital dilatation will accelerate until complete if fetal descent is unimpeded. Acknowledge~~letlts-This work was supported in part by the Department of Obstetrics and Gynecology, Indiana University Medical School, Indianapolis. Indiana. Help from Dr. R. A. Meiss is acknowledged.
REFERENCES
Embrey. M. P. and Siener, H. (1965)Cervical tokodynamo-
metry. J. Obstet. Gynuec. Br. Commonw. 72, 225-228. Friedman, E. A. and Kroll. B. H. (1969)Computer analvsis of labour progression. J. Ohster. G~wnaec.Br. Cormonw. 76, 1075-1079.
210
D. A. RICE and T. Y. YANG
Fung. Y. C. (1967) Elasticity of soft tissue in simple elongation. Am. J. Phpsiol. 213, 1532-1544. Hendricks, C. H.. Brenner. W. E. and Kraus. G. (1970) Norma) cervical dilatation pattern in late pregnancy and labor. Ant. J. Ohsret. Gynec. 106, 1065-1082. Hildebrandt. J.. Fukaya, H. and Martin. C. J. (1969) Simple uniaxial and uniform biaxial deformation of nearly isotropic incompressible tissues. Biophys. J. 9, 781-791. Kraus, H. (1967) I%rr Elastic Shells, pp. 97 and 103. Wiley. N.Y. Rice, D. A. (1974) Mechanism and measurement of cervical dilatation. Ph.D. Thesis. Purdue University. Rice, D. A., Yang. T. Y. and Stanley, P. E. (1975) A simple model of the human cervix during the first stage of labor. J. Biomechanics 9. 153-163.
Romero-Salinas, G.. Pantle. J., Aramburti, G., Figueroa. J. G.. Bienarz. J. and Caldeyro-Barcia, R. (1964) Registro continua de dilatdcion cervical en el parto human0 con
metodo electronico directo. 4rh L~ruguaran Congress of Obstetrics ant/ Grnecologr. Vol. 2. p. 718. Cited by Caldeyro-Barcia, R. and Poseiro. J. J. (1965) The powers and the mechanism of labor. In Ohsterrics. 13th Edn. (Edited by Greenhill. J. P.). p. 278. Saunders. Philadelphia. Studd. J. (1973) Partograms and nomograms of cervical dilatation in management of primigravid labor. Br. med. J. 4, 451455.
Yang, T. Y. (1973) Matrix displacement solution to elastica problems of beams and frames. Int. J. Solids Srruct. 9, 829-842. Yang, T. Y. and Lianis, G. (1974) Large displacement analysis of viscoelastic beams and frames by the finiteelement method. J. appi. me& 41. 63M40. Yang. T. Y. and Wagner. R. J. (1973) Snap buckling of nonlinearly-elastic finite element bars. J. Compur. and Srrucr. 3, 1473-1481.
Engineering
YANG
Center. Purdue University. West Lafayette. IN 47907. U.S.A.
Abstract-An exponentially elastic and viscous membrane model is developed and used to represent the human cervix during the first stage of labor in the vertex presentation. An axisymmetric orthotropic membrane finite-element is derived and evaluated. A procedure for solving for large displacement of structures confined to curved loci is derived and extended to include nonlinearly elastic and linear viscous material properties. The cervix is modeled. and the results are compared quantitatively with a previous simpler model of the cervix and qualitatively with published clinical observations. The model suggests a resolution to the controvers! over the time course of cervical dilatation.
1. INTRODUCTION It has become increasingly recognized that the progress of labor, indicated by a plot of cervical diameter as a function of time. is an important tool in the management of labor. Analysis of this plot often permits early separation of those labors which require special attention from those which are normal. The approaches used have been limited to establishing average or ‘normal’ labor progression (e.g. Studd. 1973) or to fitting measured progression to empirical equations (e.g. Friedman and Kroll, 1969). Insufficient work has been done to relate the cause of cervical dilatation, the contractions of the uterus, directly to cervical dilatation. This paper develops a mechanical model which predicts time dependent deformation of the human cervix from intrauterine pressure caused by uterine contraction. The model is in the form of an axisymmetric. orthotropic membrane finite element. This element has applicability to membraneous organs of rotational symmetry. Although this development is directed towards the particular problem of cervical dilatation, it need not be limited to the cervix. It provides the detailed information on the distribution of stresses and displacements in the membrane. The membrine of the cervix is stretched in the process of dilatation so that it conforms to the fetal head as the fetus is forced through it. Therefore. a method is derived to solve for deflection of elastic structures confined by sliding restraints and extended to include curved restraints. An iterative-incremental approach is used for the nonlinearities of geometry and elasticity and for viscoelasticity. Thus this model formuiation provides a basis for a relatively accurate mechanical representation of the cervix. In the first attempt to relate intrauterine pressure to cervical dilatation, Rice ef al. (1975) developed a ring model of the cervix. Although the ring model is an extreme geometric simplification of the birth system, high correlation between model results and
* Received 2 August 1975.
measured data for relatively short labor segments was achieved through a statistical regression technique. The results presented here show that a thin membrane ring is equivalent to the ring model and. further, that the full membrane model of the cervix covering a substantial portion of the fetal head agrees qualitatively with clinical observations of dilatation over large time periods. The following basic assumptions underlie this development : (1) The cervix, fetus and uterus act as an isolated mechanical system which is approximately axisymmetric. (21 Intrauterine fluid pressure provides the force which dilates the passive cervix. (3) Friction between the head and cervix and the inertial forces are negligible. These assumptions were made and justified by Rice er al. in their previous work. In the text of this paper, the cervix refers to all parts of the uterus below the equator of the fetal head.
2. AXISYMMETRIC,
ORTHOTROPIC,
MEMBRANE
FINITE-ELEMENT ?.I
Description of the finite-element
Following the concept of finite-element method, an axisymmetric membrane can be partitioned into a finite number of conical frusta, each of which is an axisymmetric membrane finite-element. This finite-element is described geometrically in Fig. I. All the geometric variables are defined in the figure. The two nodal edges which bound the element are parallel circles. The element is composed of an orthotropic membrane whose principal axes are in the meridional and circumferential directions. The axisymmetric conditions of geometry and loading restrict the circumferential displacement to zero, so the displacement functions are independent of the circumferential coordinate and are described in terms of a set of two-dimensional Cartesian coordinates r and z only. For each nodal circle, two displacement
201
202
D. A. RICE and T. Y.
YANG
degrees-of-freedom u and L‘in the r and : directions, respectively, are assumed. This element thus has four degrees-of-freedom
element and is found by dividing extension by the length L
where L 1 is the symbol for a row matrix.
For an orthotropic membrane in plane stress. the stress-strain relations are
2.2 Stlfiess formulation In this axisymmetric membrane element, the meridional displacement is the only displacement that contributes to the strain energy. Although two displacement degrees-of-freedom are assumed for each no&l circle, they are the r and z components of the meridional displacement. Thus only the meridional displacement function w need be assumed. It is assumed to be of linear variation along the meridional direction x,
x/w, + WJW,,
w(x) = (1 -
(2)
the meridional
where v is the Poisson’s ratio; subscripts nr and c refer to the meridional and circumferential directions, respectively; and D is the membrane rigidity defined as D,
=
E --!!Land 1 - \‘*Vc
DC= A.
E
1 - v,v,,
(7)
Substituting equations (4 and 5) into (6). the stresses in terms of the nodal displacements are obtained
L
I
The strain energy, without shearing stress, is where w, and w2 are the meridional displacements at the nodal circles (1) and (2), respectively. If 4 is used to define the meridional direction of the element, the following relations can be established, u = iw; r = iiw; r = r, + i.x A = cosfj;
z = 2, + /lx;
The circumferential &=
(3)
strain is found as L-x
r
p = sin4.
L(-_) rL
,o,;.o
J id,
(4)
U = ;
sV
(u,E, + u,,,~,,,) d I’,
(9)
where the differential volume dV = 27&(x) dx. Substituting equations (4, 5 and 8) into equation (9), the strain energy expression is obtained in terms of the four nodal degrees-of-freedom. Following Castigliano’s theorem and performing partial differentiation of the strain energy with respect to each degreeof-freedom, the stiffness equations that relate the four nodal forces (f] to the four nodal displacements {q} are obtained.
where { ] is the symbol for a column matrix, and the meridional strain is assumed constant over the
{S} =
Cklid.
(10)
The content of the stiffness matrix [k] is given in Table 1. 2.3 Evaluation
“2
Fig.
I. Description
of
the
axisymmetric
. finite-element.
membrane
The above formulation is evaluated through numerical examples to indicate its efficiency. The first example chosen is an isotropic conical membrane shell. Its apex is closed and its open edge is constrained from meridional movement. It is subjected to uniform internal pressure. The exact solution for this example is available in the text by Kraus (1967). The values assumed for the example are: internal pressure = 1 gf/cm’; modulus of elasticity = loo0 gf/cm* ; thickness = 0.5 cm: meridional length = 10 cm; and Poisson’s ratio = 0.4. The displacement shapes obtained by using 1, 2 and 10 elements are shown in Fig. 2. The lO-element solution virtually coincides with the exact solution.
A nonlinear viscoelastic membrane model Table I. Stiffness matrix for the axisymmetric orthotropic membrane finite-element
---from Kraus -
IO rlrmrnts
Fig. 3. Displacement shapes ofa spherical membrane loaded under its own weight.
R
I
-+_+‘;” i.
i2L t py
I
:-l
I’.
-1
:r
+ S‘
Q = In (rdr,)
i. = cos fj p = sin4
s, = 27rhE,/(l - V,P,)
R = (r, + r2)/2
s, = ZnhE,/(l - \Iml’,)
An indication of the convergent characteristic of this element can be obtained by varying the number of elements in the analysis. The displacements at the open edge are found to be approx. - 35. - 5.5, - 1.85, -0.8, -0.4 and -0.3% off the exact solution for the idealization of 1, 2, 3. 5, 8 and 10 elements, respect-
ively. The similar convergent trend is also noted in the results for meridional membrane stresses. Details are given in the thesis by Rice (1974). The second example is an isotropic spherical membrane loaded under its own weight as shown in Fig. 3. The opening edge is restrained from meridional movement. The values used for the example are: density per unit area = 1 g/cm2: modulus of elasticity = loo0 gf/cm2; thickness = 0.05cm; radius = I cm; opening angle = 120”; and Poisson’s ratio = 0.3. The results are found to have the same convergent trend as that noted in the first example. The displacement shape for a lOelement solution is shown in Fig. 3. Except for the cetitral element, all the nodal displacements are within 1% of the displacements computed from the exact solution in the text by Kraus (1967). The large error in the central element is expected since the element is poor when its nodal radii are small. Because the cervix representation does
also exact soluti
Fig. 2. Displacement shapes of a membrane cone under uniform pressure.
204
D. A. RICEand T. Y.
not require a zero nodal radius, the present element is adequate.
YANG
or symbolically
1Qd = 2m x 1
3. NODAL CONSI-RAINTFORMULATION
This section formulates a method of sliding constraints to force the nodal points of the finite-element to remain on a specified curve, the locus. The method is first developed for linear loci and then used with an incremental-iterative procedure for curved loci. The derivation assumes two degrees of translational freedom for each node, but its generality is not limited by this assumption.
Gel 2tnxm
Similarly, the components can be expressed as
Id)
(15b)
mxl‘
of the reaction forces {s)
cos41
‘XI‘
Yl
sin 4,
X2
cos& sin Cp2
y2
<
>
3.1 Linear loci
=
After the assemblage, the overall force-displacement stiffness equations are in the form that IF)
=
2n x 1
{QI
CKI
2nx
2nx2n
(11)
1,
where n is the number of nodal points and the capital letters are used to distinguish the system equation (11) from the element equations (10). The force vector {F) is separated as the sum of the applied load vector {P} and the reaction force vector {R} . Equation (11) thus becomes
XIII
cosAl
’ XII ’
sin 4,
or symbolically IQ
=
2tn x 1
IPI = 2n x 1
CKcIK,l 2n x 2m
I ii I
2n x 2j
-
Q/
2mxl -
2j x 1
0
’
2j x 1
(13) where the subscripts ‘c’ and ‘f designate the constrained and the free nodes, respectively. A typical constrained node ‘i’, displaced by di along a linear locus, is shown in Fig. 4. The constrained condition requires that Ui = diCOSBj = dili ri = di sin ~9~= dipi Thus the constrained be expressed as
with i = 1, 2,..., m. displacement
(4 -
[PI = CKICel 2n x m
2n x 1 i
components
Cal
2m x m
co1
tn x 1
CKJI 2n x 2,j
2j x m
R,
QC 2mxl
ISI m x 1’
Substituting equations (15b and 16b) into equation (13) yields
113 = {p) + IN = CUQI. (12) Then nodes are numbered sequentially so that the first m nodes are those to be constrained and the other (n-m) = j nodes are those left unconstrained, or free. Equation (12) is partitioned to separate the free from the constrained nodes.
[cl 2tn x m
I(
or symbolically, 10 = CR1 2n x 1 2n x 2n
(4 2n x 1’
{s) mxl
{Q/I 2jx
1I (174 UW
where [K] is the modified stiffness matrix which accounts for the linear locus boundary condition. The constrained nodal displacements {QC} and the reaction forces {RC} can be obtained by substituting the appropriate section of the solution vector {A} into equations (15 and 16), respectively.
(14) can REACTION Si
Initial position
Fig. 4. A nodal point constrained on a linear locus.
105
A nonlinear viscoelastic membrane model LOAD INCREMENT %P
ONE INCREMENTAL
STEP ---I
Fig. 5. The iterative procedure for one linear-incremental step of a nodal point constrained on a curved locus. It must be noted that the above formulation is limited only to the case of small displacements without geometric nonlinearity, such as the examples given in Figs. 2 and 3. It is also noted that such model gives the distributions of displacements and stresses (in both the circumferential and meridional direction) along the meridional direction of the membrane such as those shown in Figs. 2 and 3. The ring model developed earlier by Rice et al. (1975) provides only the dilation and the hoop force at the reference circle. It gives neither the stresses nor the distributions of displacements and stresses in the membrane surface.
3.2 Curved loci The fetal head is of an ellipsoidal shape. The membrane finite-elements are constrained to slide along
the ellipsoidal surface with large displacements. The geometrical nonlinearity caused by large displacement can be treated by a linear-incremental procedure b! revising the geometry at the beginning of every increment. Details of this incremental procedure can be found in the paper by Yang (1973). Because the locus is curved, an iterative procedure is used at ever! incremental step. This iterative procedure is described below with the aid of Fig. 5. The system shown in Fig. 5 is a spring whose lower and upper nodes are free to slide along the horizontal and curved loci, respectively. A load increment iif is applied to the upper node and the displaced position at e is sought. This system is used to simulate the constrained condition for the finite-element systern. In Fig. 5. the point o is the initial position of a nodal point at the beginning of an incremental step. When the desired load-increment is applied. the node is displaced along the tangent at o to the point h’. a trial solution. Point b on the locus is found by taking a path. as a vertical line. from h’ to the locus. A position such as b on the locus found in this way is a target point. A new constrained sliding path is defined by line z. With the node at its initial position, the same load increment is applied. The new displaced position c’ and then the new target point c are found. The process is repeated until an error tolerance is met. The reaction force is specified to be perpendicular to the locus at the target point. When the displaced upper node lies on the locus, the reaction, the applied load. and the spring tension are in equilibrium and the target point is the solution. For a finite-element system, this procedure is applied to the assembled set. 1
AX I AL
FORCE
AXIAL
Fig. 6. LoadAisplacement
\
FORCE
(prrssurc
x 5**,
gf)
relation for a linearly elastic ring sliding along an ellipsoid.
206
D. A. RICE and T. Y. YANG
To evaluate this procedure. an example of a ring constrained to slide along the surface of an ellipsoid is chosen. The ring is loaded uniformly in the tangential direction to the ellipsoid as shown in Fig. 6. This example has been solved by Rice et al. by a simple ring model. Only a single finite-element is used in this evaluation. Since Rice et al. included the stiffness of the ring in the derivation in a form that no cross sectional area need be specified, the cross sectional area in this finite-element should be assumed to be as small as possible. In this finite-element, the meridional length is assumed as 0.07cm. The thickness is implicitly included in the membrane rigidity SCwhich is assumed as 2.12 gf/cm. The average initial radius is 2 cm. The load-displacement solution obtained by both the ring model and the finite-element model are shown in Fig. 6. Total agreement is seen. During the iteration in each increment, the path to the locus from a trial solution is along the radial direction of the ellipsoid. In the unstable region the solution is not obtainable by the method proposed here. The solution in that region is, however, not desired for practical purposes. 4.
NONLINEAR
MATERIAL
In treating the nonlinear stress-strain (a vs e) relationship, the concept of tangent modulus E, is often used. This modulus is defined in a stress-strain curve as the slope of the tangent at the point with the current strain. The use of tangent modulus allows a linear-incremental stiffness formulation while simultaneously satisfying the nonlinear stress-strain relation. This concept has been successfully used by Yang
and Wagner (1973) in treating the nonlinear materials. The tangent modulus is defined in an equation form as E, = E,(E) = da(e)@. (18) In the ring example by Rice et al.. the exponentially elastic material was assumed F = A~“-‘,,’
(19)
where F is the hoop tension; (r - ro) is the extension in radius; and A and B are material constants. Substituting equations (19) into (18), the tangent modulus becomes E, = abe@
for
E > l//I,
where the strain l = (r - r,)/r,; constant fl = Br,; and constant a = A divided by the cross sectional area. Based on the material property described by equation (20), the example described by Fig. 6 is again performed. The constants a and fl are assumed as 0.25 gf/cm’ and 4, respectively, so that they conform to those constants used in the ring example. The two solutions are compared in Fig. 7. Roth virtually coincide. In performing this example, the tangent modulus used in each iterative cycle is based on the average value of the strain at the beginning of the incremental step and the strain at the current target point. 5. CERVIX
MODEL
Since it is beyond the scope of this paper to characterize cervical tissue by formulating its strain energy
5
ring
-
model
O one finite
2
0
2
4
6 AXIAL
Fig. 7. Load-displacement
(20)
6
FORCE
elamrnt
IO
(gf
relation for an exponentially elastic
12
14
I
1 ring sliding
along
an ellipsoid.
A nonlinear
axial
dilatation
--_
ring
-
IO
I
I
20
30
40
AXIAL
FORCE
relation
+ (B,,&‘]’ ‘;
I
for both linearly
function. a modification of the exponential form proposed by Fung (1967) is proposed. The tissue is assumed to satisfy the following conditions: (1) The membrane is orthotropic with zero stress across its thickness. (2) There is no shear strain along the tissue axes. (3) The tissue is incompressible (v = 0.5). (4) The strain energy equals zero for the state of zero strain. The tangent modulus is assumed as (E,); = cc/?,exp[(&,.?
model
force
IO Fig. 8. Load-radial
viscoelastic membrane
i = c. m, (211
where r and /? are material constants; the subscripts c and nr refer to the circumferential and meridional directions, respectively. For the case of isotropic material in uniform biaxial strain, equation (21) reduces to the same form as that of equation (20). The cervix covers the fetal head for 85”” of the axial length below the equator and is partitioned into 10 finite elements. The head below the equator is an ellipsoid with the equatorial radius a = 5.0cm and axial radius h = 7.0cm. Above the equator the locus is assumed to extend indefinitely with a constant radius equaling the equatorial radius. The fetal head and body are thus represented by a circular cylinder with an ellipsoidal end. The cervix is assumed as an elastic membrane with a constant thickness h = 0.5 cm. The lowest node. node 1. has an initial radius of 2.0 cm. The last node. node 11. is located initially 0.5 cm below the equator. Load is applied in increments and the iterative procedure is used to solve for displacement. The material constants used for linear material are: E, = E, = 40 gf/cm* and 1’= 0.5. Those for nonlinear material are: J = 0.25 gf/cm2; 8, = fl,,, = 4; and v = 0.5. The maxi-
model finite-element
I
1
60
50
I
70
00
(gf 1 and nonlinearly
elastic membrane
cervix
mum load for the parameters specified is equivalent to an intrauterine pressure of 1gfjcm’. The results are shown in Fig. 8. A study by varying the number of elements shows that the results converged to the ones presented in Fig. 8 at the 7 element level.
6. COMPARISON
OF RING AND
MEMBRANE
Assuming that the proposed strain energy function of the membrane is representative of cervical tissue it is reasonable to expect that the load-ilisplacement relation of the model is similar to the pressure-dilatation relation of the biological system. Further vcrification of this is shown by Fig. 8 where the pressureradius solution for the ring described by Rice it ul.. is successfully superimposed on the IO element membrane solution. The parameter values of the ring. except for the multiplicative constant, are similar to the values found experimentally by Rice or ill. The! are: A = 0.215gf/cm*; B = 2cm-‘: and I’(,= Zcm.
7. VISCOELASTICITY
The creep present in cervical tissue during dilatation is handled by an extension of the incrementaliterative procedure. To represent this creep as a series linear element, the viscous strain .C is expressed as 1’51 where q is the viscous coefficient
and t
is
lhe time.
D. A. RICE and T. Y. YANG i-
and equatorial radius = Scm. initial node I radius = 2.0cm. In order to compare the qualitative model results with physical observation. a time varying axial load in the form that axial force = 1 + 10 sin2(0.2t)(gf)
0
20
60
60
40 TIME
100
is used to approximate the periodic nature of the intrauterine pressure. This load is applied to the same model as described in Fig. 9. The solution is shown in Fig. 10. The regression technique used with ring model by Rice et al. is applied to this solution. The multiple correlation coefficient is y > 0.99. The use of statistical technique here is justified because the convergence criterion of the membrane model allows some error in node position. As noted by Hildebrandt et al. (1969), the results of biaxial membrane strain are not predictable from uniaxial experiments, so the parameters of the ring model are not equivalent to the membrane parameters.
120
(s*condr)
Fig. 9. Radial dilation at constant loads due to the viscous effect. This viscous strain is obtained by using a trapezoidal
8. RESULTS
approximation, and is used to calculate the viscous displacements. A finite-element linear-incremental procedure has been used by Yang and Lianis (1974) to solve the geometrically nonlinear large displacement problems of viscoelastic beams and frames. The same concept and procedure are used here. Since the membrane is constrained to slide along a curved locus, the aforementioned iterative procedure, associated with the use of equation (17), is used in every incremental step. Two constant loads, 15 and 30% of the maximum loads shown in Fig. 8, are applied until creep is sufficient to allow displacement to proceed to the unstable region. The results are shown in Fig. 9. The following values were used: qC= Q,, = lo4 gf-sec/cn?; a = 0.25 gf/cm’; /3, = /?,,,= 4; ellipsoid axial radius = 7 cm;
DILATION
Friedman and Kroll (1969) suggest that the dilatation rate decelerates to near zero as dilatation approaches its maximum. Hendricks et al. (1970) suggest that the dilatation rate accelerates until full dilatation. Studd (1973) presents average dilatation curves which support neither of these interpretations, but as shown by Figs. 9 and 10, the membrane model supports the interpretation of Hendricks et al. The conflict between these two interpretations is resolved if it is assumed that for the larger dilatations the head has descended sufficiently to begin stretching structures other than the cervix. The model formulation limits the system to the cervix, uterus and fetus so such descent violates the basic assumptions, and the model is no longer applicable. If Hendrick’s patients
RADIUS
100
TIME
(26)
( seconds)
Fig. 10. Ten element simulatjon of tbe exponentially elastic and viscous membrane model.
209
A nonlinear viscoelastic membrane model such that the cervix is higher in the pelvis or the pelvis is larger than that for Friedman’s patients. the model formulation applies to full dilatation and predicts the accelerating curve. The position of the cervix. whether higher or lower in the pelvis. may be the determining factor for the shape of the dilatation time curve. Embrey and Siener (1965) make several observations concerning the nature of cervical dilatation. Only one of these observations. the passive nature of the cervix. is assumed by the model. Model results are compared to each observation below. 1. Small progress is made for each contraction, but dilatation increases during a contraction. The model shows an elastic stretch during a contraction and a small viscous strain during the course of a contraction (Fig. 10). This is due to the viscoelastic nature of cervical tissue. 3. A contraction causes dilatation in proportion to its strength. This statement implies a monotonic relationship between pressure and dilatation, not necessarily a linear relationship. The model shows a monotonic relationship (Fig. 8). 3. A larger total dilatation results in a larger elastic dilatation for a contraction of given strength. Figure 10 definitely confirms this. 4. There is a delay between the onset of a contraction and dilatation which lessens as dilatation increases. The model shows no such delay. 5. The drug oxytosin causes greater contractions and thus faster dilatation. Oxytosin’s action is well known. That greater pressure causes faster dilatation is shown by Fig. 9. The models predict all the observations except (4). The delay could be a_ measurement artifact or the model may be inadequate to describe the system in this respect. Data from Romero-Salinas et al. (1964) and data collected by Rice er al. show no such effect. were selected
9. DISCL’SSION
The model developed herein for membrane elastica problems is restricted only by the requirements of axisymmetry and quasi-static equilibrium. This model should. therefore. have applicability to many biomechanics problems involving a large class of hollow membrane organs. Because the development allows for the finite element nodes to be either free or confined by sliding restraints. these organs could be filled by solid objects. These organs could also be filled by fluids or a combination of solids and fluids if the stiffness formulation is modified to account for the fluid pressure. The geometric shape of the solid contained within the membrane must be axisymmetric although it need not be defined by a mathematically continuous function. The examples presented assumed an ellipsoid, but a tabulated function should do as well. The assumption in the theoretical development that zero friction occurs at the sliding restraints is not as
restrictive as it may seem. By appropriately specifying the angle of the reaction force to the locus (to be other than perpendicular) the effects of static, kinetic. or viscous friction can be simulated. The cervix model uses the membrane development to represent the anatomical structure of the cervix. Given a simulated intrauterine pressure. the model produces a numerical solution for dilatation which agrees well with published qualitative observations which encompass a large portion of labor time. Further, this solution shows agreement with an earlier one-dimensional ring model which is highly significant in explaining quantitative data. Since the parameters that must be assumed in the membrane model are considerably more than those in the ring model. direct quantitative comparison with experimental data is not performed. To find those parameters is beyond the scope of this paper. This cervix model thus completes the bridge between measured data, dilatation and pressure. and the morphology of the cervix during the first stage of labor. The major limitation of the model is that there is presently no way to evaluate its parameters (initial length, elastic coefficient. viscous coefficient. Poisson’s ratio for each direction for each element) noninvasively. Given reasonable assumptions for the parameter values, however, much information about the state of the cervix may be inferred: the distribution of stress and strain; the nodal positions; and the pressure distribution over the fetal head as a function of intrauterine pressure and time, 10. CONCLUDING
REMARliS
A membrane finite element model was developed which has general applicability to axisymmetric hollow organs having solid contents. This model was shown to be suitable for the problem of human cervical dilatation. Using an exponentially elastic and viscous membrane, the model displays many of the qualitative characteristics observed in the cervix. Further, the model correlates well with an earlier simple. but quantitative, ring model of the cervix. These facts lend credence to the soundness of the ring model and to the ability of the membrane model to represent the cervix during dilatation. The model suggests that ccrvital dilatation will accelerate until complete if fetal descent is unimpeded. Acknowledge~~letlts-This work was supported in part by the Department of Obstetrics and Gynecology, Indiana University Medical School, Indianapolis. Indiana. Help from Dr. R. A. Meiss is acknowledged.
REFERENCES
Embrey. M. P. and Siener, H. (1965)Cervical tokodynamo-
metry. J. Obstet. Gynuec. Br. Commonw. 72, 225-228. Friedman, E. A. and Kroll. B. H. (1969)Computer analvsis of labour progression. J. Ohster. G~wnaec.Br. Cormonw. 76, 1075-1079.
210
D. A. RICE and T. Y. YANG
Fung. Y. C. (1967) Elasticity of soft tissue in simple elongation. Am. J. Phpsiol. 213, 1532-1544. Hendricks, C. H.. Brenner. W. E. and Kraus. G. (1970) Norma) cervical dilatation pattern in late pregnancy and labor. Ant. J. Ohsret. Gynec. 106, 1065-1082. Hildebrandt. J.. Fukaya, H. and Martin. C. J. (1969) Simple uniaxial and uniform biaxial deformation of nearly isotropic incompressible tissues. Biophys. J. 9, 781-791. Kraus, H. (1967) I%rr Elastic Shells, pp. 97 and 103. Wiley. N.Y. Rice, D. A. (1974) Mechanism and measurement of cervical dilatation. Ph.D. Thesis. Purdue University. Rice, D. A., Yang. T. Y. and Stanley, P. E. (1975) A simple model of the human cervix during the first stage of labor. J. Biomechanics 9. 153-163.
Romero-Salinas, G.. Pantle. J., Aramburti, G., Figueroa. J. G.. Bienarz. J. and Caldeyro-Barcia, R. (1964) Registro continua de dilatdcion cervical en el parto human0 con
metodo electronico directo. 4rh L~ruguaran Congress of Obstetrics ant/ Grnecologr. Vol. 2. p. 718. Cited by Caldeyro-Barcia, R. and Poseiro. J. J. (1965) The powers and the mechanism of labor. In Ohsterrics. 13th Edn. (Edited by Greenhill. J. P.). p. 278. Saunders. Philadelphia. Studd. J. (1973) Partograms and nomograms of cervical dilatation in management of primigravid labor. Br. med. J. 4, 451455.
Yang, T. Y. (1973) Matrix displacement solution to elastica problems of beams and frames. Int. J. Solids Srruct. 9, 829-842. Yang, T. Y. and Lianis, G. (1974) Large displacement analysis of viscoelastic beams and frames by the finiteelement method. J. appi. me& 41. 63M40. Yang. T. Y. and Wagner. R. J. (1973) Snap buckling of nonlinearly-elastic finite element bars. J. Compur. and Srrucr. 3, 1473-1481.