Journal of Applied Mathematics

Volume 2012, Article ID 985349, 18 pages

http://dx.doi.org/10.1155/2012/985349

## Strength Analysis Modelling of Flexible Umbilical Members for Marine Structures

^{1}Division of Structural Engineering, MARINTEK, P.O. Box 4125 Valentinlyst, 7450 Trondheim, Norway^{2}Department of Marine Technology, NTNU, 7491 Trondheim, Norway

Received 10 February 2012; Accepted 30 April 2012

Academic Editor: Carl M. Larsen

Copyright © 2012 S. Sævik and J. K. Ø. Gjøsteen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This paper presents a 3-dimensional finite element formulation for predicting the behaviour of complex umbilical cross-sections exposed to loading from tension, torque, internal and external pressure including bending. Helically wound armours and tubes are treated as thin and slender beams formulated within the framework of small strains but large displacements, applying the principle of virtual displacements to obtain finite element equations. Interaction between structural elements is handled by 2- and 3-noded contact elements based on a penalty parameter formulation. The model takes into account a number of features, such as material nonlinearity, gap and friction between individual bodies, and contact with external structures and with a full 3-dimensional description. Numerical studies are presented to validate the model against another model as well as test data.

#### 1. Introduction

Flexible risers such as flexible pipes and umbilical cables represent a crucial element of a floating production system. The concept is an attractive alternative to a rigid riser since it does not require heave compensation and tensioning devices at the top or riser manifold at the seabed. At the same time, it offers ease of installation, retrieval and usage elsewhere.

Flexible pipes and umbilical cables both work as composite pipes that are compliant and highly deformable in bending but strong and stiff in response to internal pressure, external pressure, tension, and torque.

For the flexible pipe, concentric polymeric layers are used to provide sealing. These layers are supported by interlocked metallic layers to resist pressure loading and tensile armour layers to resist tension, torque, and the pressure end cap effect. Umbilicals have homogenous layers, cables, hydraulic hoses, and metal tubes at their central core. Both helical armours and helical multifunction tubes are used to resist tension and torque; see Figure 1.

Predicting the structural response in such composite structures is a problem that has been dealt with by many authors. Knapp [1] addressed the axisymmetric response of aluminium core with steel reinforcement cables. Raoof and Hobbs [2] proposed to use an orthotropic material formulation to treat the helix armours as shells. Costello [3] developed a model for wire ropes that is suitable for cables with few large tendons in lateral contact. Witz and Tan [4] presented expressions for cables and flexible pipes under moderate loads. Computer models for the same range of application on flexible pipes were also presented by McNamara and Harte [5]. Feld [6] investigated the response of power cables, material nonlinearity included, but assuming a constant factor between tension and core radial contraction. Custódio and Vaz [7] studied the axisymmetric response of umbilical cables, introduced improvements to previously published models, and compared estimations with experiments. Sævik and Bruaseth [8] presented a fairly general 2-dimensional model, however, also assuming axisymmetric behaviour.

The response of flexible pipes subjected to bending has been less focused on. Lutchansky [9] and Tan et al. [10] allow axial movement of the helical reinforcing only. Ferét and Bournazel [11] suggest that the unbonded helically wound tendon will follow the geodesic of the curved cylinder. Leclair and Costello [12] used Love's equations and assumed wire geometry to calculate the local and global response of wires in a bent rope. Sævik [13] developed a curved beam element to study the stress and deformation in one single armour tendon exposed to a given curvature distribution. Most models presented in the literature, however, assume that the curvature is constant and given from global analyses assuming elastic bending stiffness, excluding the effect of the friction work done due to relative sliding between layers. Sævik et al. [14] presented a model considering the full cross-section in bending for flexible pipes; however, the assumptions applied in this model utilised the benefit from the layered structure of such pipes, thus not very well suited for the umbilical cross-section, often consisting of a few individual structural elements with a significant bending stiffness.

The main purpose of the present work has therefore been to formulate a model that can describe stresses and slip of umbilical structural elements for both axisymmetric and bending loads and to validate the model against available test data.

#### 2. Model

##### 2.1. General

Based on a comprehensive literature review of models, Custódio and Vaz [7] listed the most common assumptions. For clarity and comparison purposes, these assumptions are repeated here using their definition.(1)Regularity of initial geometry: (a) the homogenous layers are long and uniform cylinders; (b) the wires are wound on a perfectly cylindrical helix; (c) the wires are equally spaced; (d) the wires of an armour are numerous, hence the forces they exert on the adjacent layers may be replaced by uniform pressure; (e) the structure is straight.(2)Reduction to simple plane analysis: (f) there are no field loads such as self-weight; (g) end effects may be neglected; (h) the material points from any layer have the same longitudinal displacement and twist; (i) all wires of an armour present the same stress state; (j) the wires maintain a helicoidal configuration when strained; (k) the angle between the wire cross-section principal inertia axis and a radial vector linking the centre of structure's cross-section and the centre of wire's cross-section is constant; (l) there is no overpenetration or gap spanning.(3)The effects of shear and internal friction are neglected: (m) the wires are so slender that the movements of the material points are governed only by their tangent strain and not by the change in curvature.(4)Linearity of the response: (n) the materials have linear elastic behaviour; (o) the changes in armour radii and pitch angles are linearly small; (p) the wires in one armour never touch laterally or are always in contact; (q) there are no voids between layers nor among cables in the functional core; (r) the homogeneous layers are thin and made of soft material so they simply transfer pressure; (s) the umbilical's core responds linearly to axisymmetric loading; (t) both loading and response are not timedependent.

The model presented by Custódio and Vaz avoids the assumptions (n)–(s). The present model, however, avoids the assumptions (c), (d), (f), (g), (h), (i), (j), (k), (m), (n), (o) and (p) enabling contact to be handled on component level as well as allowing end effects and friction to be included. However, in order to limit the number of unknown parameters in a 3D model, the governing equations are formulated to only include 1st order helices. It is also assumed that local deformations in the cross-section can be handled by a surface stiffness penalty parameter, so that the structural response can be described by beam theory. The finite element approach has further been selected as it allows structural elements and contact interaction effects to be handled on an individual basis.

##### 2.2. Finite Element Formulation

In order to establish the equilibrium equations, the principle of virtual displacement is applied. For an arbitrary equilibrium state where is the stress tensor, is the initial stress tensor, is the strain tensor, is the surface traction, and is the displacement vector. may be obtained from the initial strain tensor by applying the material law. In the case of nonconservative loading such as pressure, the resulting load will change as a function of the area change. Hence, the change in the surface area of the volume has to be formulated as a function of the strain. The above equation is used for equilibrium control. However, the equation also needs to be formulated on incremental form to allow equilibrium iterations to be carried out. The incremental form is obtained as where is the Green-St. Venant strain tensor; see Belytschko et al. [15]. Equation (2.2) gives the incremental equilibrium equation to be used as basis for the stiffness matrix. The first term gives the material stiffness matrix, whereas the second term gives the geometric stiffness matrix.

With respect to choice of reference configuration from which deformations are measured, the so called *corotational ghost reference* formulation has been chosen. The basic idea is to separate the rigid body motion from the local or relative deformation of the element. This is done by attaching a local coordinate system to the element and letting it continuously translate and rotate with the element during deformation; see Horrigmoe and Bergan [16]. The basis for this procedure is by attaching a local coordinate system to each node of the structure. Along the helix, the following transformation describes the position of the local coordinate system, , relative to the right-handed Cartesian coordinate system, , located at the cross-section centre; see Figure 2:
where is the lay angle and is the angular position.

During deformation the helix nodes will translate and rotate and in order to update the position of the local node base vector system from one configuration, , to the next, +1 the assumption of small incremental rotations is applied. This gives the following relation between the updated local node base vector system relative to the global base vectors: where the incremental rotation matrix is defined by the rotation increments about the global axes for each node and where the rotation is carried out in subsequent order 1-2-3, which by multiplication of the three associated rotation matrices gives where denotes and denotes .

The element types needed can be divided into two categories: (1)beam elements to model copper conductors, tubes, fill materials or user defined structural elements,(2)contact elements to model contact between different element combinations such as contact between core and 1st helix layer, between two helix layers and between outer helix, and other interacting structural elements are examples.

With reference to (2.1) and (2.2), the following is needed in order to develop equations for each element type that can be implemented into a computer code:(1)Kinematics description, that is, a relation between the displacements and rotations and the strains at a material point. (2)A material law connecting the strain with resulting stresses.(3)Displacement interpolation, describing the displacement and rotation fields by a number of unknowns on matrix format that can be programmed.

In the following the above will be separately addressed.

##### 2.3. Beam Kinematics

The beam contribution to the equilibrium equation is established assuming that the Bernoulli-Euler and Navier hypotheses apply. The Green-St. Venant strain tensor is used as strain measure when formulating the incremental equilibrium equations. The second order longitudinal strain term in the Green-St. Venant strain tensor is neglected to avoid shear locking. However, all terms related to coupling between longitudinal strain and torsion are included.

The displacements of an arbitrary point , defined by local coordinates in the cross-section as seen in Figure 2, may be expressed as Introducing the previously mentioned assumptions, the longitudinal Green-St. Venant strain is found to be where is the displacements along the respective axis and is the rotation about the same axis . The subscript “0” quantities represent the displacement field in the respective directions. The shear deformation, , is assumed according to St. Venant’s principle assuming the beam to be long and slender with no end section warping as The beam deformations resulting from the motion of the nodal system are referred to the corotational ghost reference -system. For an arbitrary equilibrium state the ghost reference system is defined by the following procedure, where the superscript index refers to element end node: The local beam deformations at end 1 and referred to as the local system, is found by transformation of the node system into the element system: The elemental rotational deformations in node are determined by where represents the resultant rotation at node .

##### 2.4. Contact Kinematics

Reference is given to Figure 3.

A contact element with two nodes is developed by considering two bodies *A* and *B*1 for this case body *B*2 is ignored. Each of the bodies occupies a region and has a boundary where or . For a virtual time interval [], the displacement fields are denoted by , where . When bodies and 1 are brought into contact, let be the unknown contact surface, which satisfies the relationship and , where denotes part of the surface with prescribed surface tractions and denotes part of the surface with prescribed surface displacements. Also let be the outward surface normal vector of body *A* at and be the corresponding tangent vector. At the beginning of a time increment, an initial gap at in the direction of is defined as
where , or , represents the updated coordinates of a point at time .

In order to describe contact between two helical layers, contact is considered between 3 bodies. One body is from the inner layer, and is denoted . The remaining two bodies are from the layer outside, and are denoted *B*1 and *B*2. With reference to Figure 3, the initial gap can be expressed by:
where the non-dimensional parameter is defined as:
Note that in this case, the direction of the surface normal vector **n** is from the center of the umbilical to the center of body *A*. It is directed outwards, to the next layer containing bodies *B*1 and *B*2.

The current gap at time in the direction of can be described by: for a two-body contact and by for a three-body contact.

Two contact conditions may occur:(1)gap opening, if ,(2)contact, if .

Further, if contact has been established, relative slippage including friction work will occur for a two-body contact when
where is pointing along the centre line of body *A* and is obtained by taking the cross-product between and . The three-body contact is treated in the same manner:

###### 2.4.1. Patch Test

The meshing in the longitudinal direction is restrained to be identical for all the helical components. Hence, the 2-noded contact elements operate between coincident nodes, and the contact element will pass the patch test [17]. For the 3-noded contact elements, the force balance will be satisfied within each separate contact element. A contact force that is acting on body , Figure 3, will have a reaction force that will be linearly distributed with between body *B*1 and *B*2. The patch test is therefore also passed for this contact element.

##### 2.5. Material Models

###### 2.5.1. Beams and Tubes

Denoting the longitudinal direction by index 1 the circumferential direction by index 2 and assuming that the normal direction is governed by a prescribed value, , Hooke’s law for linear elastic materials reads: For the beam and tube elements (2.19) is further simplified by prescribing by application of thin shell theory for the tube element.

When the stresses exceed the elastic limit of the stress/strain relation, an elastoplastic formulation is required. In this case, it is of great importance to use an elastoplastic formulation that takes into account the two-dimensional stress state, that is, both the stresses in the axial and hoop directions for the tubes. The applied elastoplastic material model is based on expressing(i)the yield criterion as a yield surface specified by the scalar function ,(ii)the 2nd Piola-Kirchhoff stress tensor as a measure for stress together with the energy conjugate Green strain tensor as a measure of strain.

-flow theory of plasticity is applied, assuming that the yielding is independent of the first and third deviatoric stress invariant, and any combination of kinematic and isotropic hardening is allowed for. The constitutive relation reads: where where is the elastic material law, is the plastic work done, and is the equivalent stress. For further details, references are made to Levold [18], McMeeking and Rice [19] and Belytschko et al. [15].

###### 2.5.2. Contact Element

In nonlinear finite element analysis, there are three commonly used principles when dealing with contact problems see Shyu et al. [20, 21]. These are(i)the lagrange multiplier method (LM),(ii)the penalty method (PM),(iii)the mixed method (MM).

In LM, the constraint conditions for a contact problem are satisfied by introducing Lagrange parameters in the variation statement, where both the displacements and Lagrange multipliers are treated as unknowns which again leads to an increased number of finite element equations.

In PM, the contact pressure is assumed proportional to the amount of penetration by introducing a pointwise penalty parameter. The final stiffness matrix does not contain additional terms. However, since the resulting contact forces are of the same order as the assumed displacement field, this may lead to violation of the local Babuska-Brezzi stability condition, see Chang et al. [20], and the success of using MM is highly dependent on the selected order of the contact pressure. A guideline on this item is given in [20] where an MM contact element with excellent numerical properties was presented.

In this case, however, the focus has been on formulating a contact element that allows easy and flexible modelling of contact effects in helical structures. Therefore a node-to-node contact formulation based on PM has been adopted.

The constitutive relation used to model friction in the contact elements consists of two major ingredients:(i)a friction surface,(ii)a slip rule.

The constitutive model used is established based on the work by Shyu et al. [21] but is adjusted to include material hardening to improve numerical stability. The friction surface is assumed to be a function of the normal traction , the tangential tractions and , and the friction coefficient , that is: Assuming a linear surface, we have The slip increment is divided into two parts: where is the 2-dimensional increment of elastic slip and is the corresponding increment of plastic slip. The elastic slip increment is determined as where is a proportional constant. It is seen that as grows towards infinity, a classical Coulomb law is obtained. Assuming symmetry, that is, an associative slip rule, the following slip rule is postulated: where is a still unknown proportional constant. By assuming a constant normal traction, the proportional constant is determined from the consistency condition: where the last term is associated with material hardening. The expressions are therefore developed including the hardening parameter , that is:

By combining the previous equations the following constitutive relation is obtained:

##### 2.6. Solution Procedure

The developed finite element equations have been implemented into a computer code using standard procedures for matrix operations. For the beam elements, linear interpolation is used in the axial direction, whereas cubic interpolation has been applied in the transverse direction. For the contact elements, linear interpolation is used for the 3-noded element used between helical layers to allow the helix to distribute forces if the helix is positioned between two bodies in the next layer. Between umbilical helix and 1st helix layer and between last helix layer and external structures, two-noded contact elements are used, and hence no interpolation is applied for these elements. With respect to numerical integration, closed form expressions are developed for the beam elastic case; that is, no numerical integration is needed. In the elastoplastic case, for pipe elements numerical integration is performed by dividing the beam into a number of subvolumes each characterized by its strain/stress state.

The numerical procedure used to solve the finite element equations was based on a Newton-Rapson iterative scheme controlled by applying displacement, force, and energy norms to ensure convergence, see Belytshcko et al. [15].

#### 3. Case Study

The purpose of the case study was to investigate the performance of the developed model compared to others. The case selected is the umbilical test case presented by Lutchansky [9], and the basic properties are summarized in Table 1.

The contact between the beam and the armour was modelled by setting the surface stiffness penalty parameter to 7 MPa/mm, except for the armour to armour contact, where 100 MPa/mm was applied. Table 2 shows the performance of the present model as compared to test data, the Custódio and Vaz model, and the 2-dimensional model presented by Tan et al. [10]. It is seen that the predicted values are similar to the other model values with respect to stiffness parameters. The present 3-dimensional model suggests stronger torsion coupling. However, few data points on coupling parameters were included in the measurements for this case study.

For the test cases applied in this study a linear convergence was observed for the complete model.

#### 4. Experimental Studies and Analyses

The experimental studies were performed for validation of a 2D umbilical model and have been reported in this context by Tan et al. [10].

Different umbilical cross-sections was tested with respect to the axial stiffness under tension, coupling between axial strain and torsion and torsion stiffness in both directions. The different umbilical cross-sections cover the range from torsion balanced dynamic umbilicals to torsion unbalanced static umbilicals and with variable complexity. The cross-section details required to define input to numerical models are attached in the appendix.

Measurements, predictions, and comparison between experimental data and model results are presented in Table 3. In all analyses a surface stiffness parameter of 100 MPa/mm has been applied to all metal-to-metal interfaces. For the other surfaces, a stiffness parameter of 7 MPa/mm has been applied. An exception was done for cross Section 1, which contains layers of fillers creating a soft core. For this cross section the 7 MPa/mm surface stiffness value was reduced to 3.5 MPa/mm. When obtaining the torsional stiffness and corresponding coupling factors, 100 kN tension was applied in the analyses. The numerical model was 1.5 m long, and the element lengths were 1.5 cm. A typical element meshing is shown in Figure 4.

An example of measured and predicted data is shown in Figure 5, where the axial force is plotted against axial strain for umbilical 1 and umbilical 4. Figure 5 shows three graphs, measured, computed, and analytical where the analytical is based on simply adding together the axial stiffness of all structural elements, taking the lay angle into account but neglecting the effect of radial interaction. It is seen that good correlation is obtained between the measured and the computed results in terms of stiffness. For umbilical 1, the analytical stiffness is considerable higher than the simulations and the tests. This is likely to be due to radial interaction effects, which are included in the model, but cannot be taken into account in the analytical solution. Umbilical 1 has a soft core, and radial contraction is therefore more significant than for a stiffer cross-section, as umbilical 4.

It is seen that the measured stiffness is low for the first load cycle of umbilical 1, which most probably is due to gap closure.

In addition to test and analysis results, Table 3 gives a summary of the relation tested value/predicted value for all umbilical cross-sections for all test parameters, including the statistics in terms of average correlation and standard deviation.

It is seen that in terms of axial and torsion stiffness good correlation is found, with a standard deviation within a 10% range. This must be considered reasonable for such complex structures. It is noted that the axial stiffness in average is predicted to be 7-8% lower than the measurements. With regard to tension/torsion coupling the number of valid measurements was limited, and therefore the statistical basis is too weak to conclude. This was due to the fact that some of the specimens were torsion balanced; hence very small inherent test rig friction forces would influence the measured rotation. The reported values include only torsion unbalanced cross-sections where the amount of rotation is larger. For these cross-sections reasonable correlation was found with a standard deviation within a 20% range.

It should be noted that a better match of analysis to the experiment could have been obtained if the surface stiffness parameter had been adjusted for each cross-section. This would also be a reasonable approach, as the surface stiffness represents the components ability to get a local depression at the contact point. This is clearly different for different materials and geometries. However, the results indicate that using the suggested stiffness is sufficient for most purposes.

#### 5. Conclusions

In this paper, a 3-dimensional model for simulation of complex umbilical cross-sections were presented and tested with respect to published data for axisymmetric response. The model was based on the finite element approach, enabling contact interface effects to be handled on individual component level and in 3 dimensions. Material nonlinearities, gap formation, friction, lateral contact between wires, contact with external structures, and wires' curvature change are taken into account. The model can estimate the stresses and the displacements of individual structural element as well as the overall structural response. The influence of gap formation effects between structural elements was handled by a surface stiffness penalty parameter applied to the contact elements, and generally good correlation was found with test data. This demonstrates model robustness at least for the range of umbilical cross-sections used in the experiments. The model enables generality with respect to external loading, and work is ongoing with respect to investigating the model performance with respect to bending and friction stresses as well.

#### Appendix

#### Cross-Section Model Input Details

For more details see Tables 4, 5, 6, 7, 8, 9, 10, 11, and 12.

#### Acknowledgments

The authors wish to acknowledge Statoil, Shell Norge, BP, Petrobras, and Nexans Norway for allowing the paper to be published.

#### References

- R. H. Knapp, “Derivation of a new stiffness matrix for helically armoured cables considering tension and torsion,”
*International Journal for Numerical Methods in Engineering*, vol. 14, no. 4, pp. 515–529, 1979. View at Google Scholar · View at Scopus - M. Raoof and R. E. Hobbs, “Analysis of multilayered structural strands,”
*Journal of Engineering Mechanics*, vol. 114, no. 7, pp. 1166–1182, 1988. View at Google Scholar · View at Scopus - G. A. Costello,
*Theory of Wire Rope*, Mechanical Engineering Series, Springer, New York, NY, USA, 1990. - J. A. Witz and Z. Tan, “On the axial-torsional structural behaviour of flexible pipes, umbilicals and marine cables,”
*Marine Structures*, vol. 5, no. 2-3, pp. 205–227, 1992. View at Google Scholar · View at Scopus - J. F. McNamara and A. M. Harte, “Three dimensional analytical simulation of flexible pipe wall structure,” in
*Proceedings of the 8th International Conference on Offshore Mechanics and Arctic Engineering (OMAE '89)*, pp. 477–482, March 1989. View at Scopus - G. Feld,
*Static and cyclic mechanical behaviour of helically-wound subsea power cables [Ph.D. thesis]*, Faculty of Engineering, Herriot-Watt University, November 1992. - A. B. Custódio and M. A. Vaz, “A nonlinear formulation for the axisymmetric response of umbilical cables and flexible pipes,”
*Applied Ocean Research*, vol. 24, no. 1, pp. 21–29, 2002. View at Publisher · View at Google Scholar · View at Scopus - S. Sævik and S. Bruaseth, “Theoretical and experimental studies of the axisymmetric behaviour of complex umbilical cross-sections,”
*Applied Ocean Research*, vol. 27, no. 2, pp. 97–106, 2005. View at Publisher · View at Google Scholar · View at Scopus - M. Lutchansky, “Axial stresses in armor wires of bent submarine cables,”
*Journal of Engineering for Industry, Transactions of the ASME*, vol. 91, pp. 687–693, 1969. View at Google Scholar - Z. Tan, J. A. Witz, G. J. Lyons, J. Fang, and M. H. Patel, “On the influence of internal slip between component layers on the dynamic response of unbonded flexible pipe,” in
*Proceedings of 10th International Conference of Offshore Mechanics and Arctic Engineering Conference (OMAE '91)*, 1991. View at Scopus - J. J. Ferét and C. L. Bournazel, “Calculation of stresses and slip in structural layers of unbounded flexible pipes,”
*Journal of Offshore Mechanics and Arctic Engineering*, vol. 109, no. 3, pp. 263–269, 1987. View at Google Scholar · View at Scopus - R. A. Leclair and G. A. Costello, “Axial, Bending and Torsional Loading of a Strand with Friction,” in
*Proceedings of the 6th International Offshore Mechanics and Arctic Engineering Symposium (OMAE ’87)*, Houston, Tex, USA, March 1987. - S. Sævik,
*On stresses and fatigue in flexible pipes [Ph.D. thesis]*, NTH, 1992. - S. Saevik, E. Giertsen, and G. P. Olsen, “New method for calculating stresses in flexible pipe tensile armours,” in
*Proceedings of the 17th International Conference on Offshore Mechanics and Arctic Engineering (OMAE '98)*, Lisbon, Portugal, July 1998. View at Scopus - T. Belytschko, W. K. Liu, and B. Moran,
*Nonlinear Finite Elements for Continua and Structures*, John Wiley & Sons, Chichester, UK, 2000. - G. Horrigmoe and P. G. Bergan, “Incremental variational principles and finite element models for nonlinear problems,”
*Computer Methods in Applied Mechanics and Engineering*, vol. 7, no. 2, pp. 201–217, 1976. View at Google Scholar - N. El-Abbasi and K. J. Bathe, “Stability and patch test performance of contact discretizations and a new solution algorithm,”
*Computers and Structures*, vol. 79, no. 16, pp. 1473–1486, 2001. View at Publisher · View at Google Scholar · View at Scopus - E. Levold,
*Solid mechanics and material models including large deformations [Ph.D. thasesis]*, Division of Structural Mechanics, The Norwegian Institute of Technology, NTH, Trondheim, Norway, 1990. - R. M. McMeeking and J. R. Rice, “Finite-element formulations for problems of large elastic-plastic deformation,”
*International Journal of Solids and Structures*, vol. 11, no. 5, pp. 601–616, 1975. View at Google Scholar · View at Scopus - T. Y. Chang, A. F. Saleeb, and S. C. Shyu, “Finite element solutions of two-dimensional contact problems based on a consistent mixed formulation,”
*Computers and Structures*, vol. 27, no. 4, pp. 455–466, 1987. View at Google Scholar · View at Scopus - S. C. Shyu, T. Y. Chang, and A. F. Saleeb, “Friction-contact analysis using a mixed finite element method,”
*Computers and Structures*, vol. 32, no. 1, pp. 223–242, 1989. View at Google Scholar · View at Scopus