advancing contact of structural elements simplified

computers & Svucrurrs Vol. 33, No. 2. pp. 50%512, Printed in Great Britain.

1989 0

0045-7949189 53.00 + 0.00 1989 Pergamon Press plc

RECEDING/ADVANCING CONTACT OF STRUCTURAL ELEMENTS SIMPLIFIED N. J.

SALAMON

Department of Engineering Science and Mechanics, The Pennsylvania State University, University Park,

PA 16802, U.S.A. (Receioed 30 Augusf 1988) Abstract-An extremely simple algorithm is developed for simulation of structural elements in receding/advancing, unilateral contact with independent constraints. The theory is systematically formulated to prove that the nonlinear solution depends upon a single free variable which can be kinematically determined. A FORTRAN subroutine is given which illustrates the elementary nature of the new development. Results from its implementation tally with those given previously.

INTRODUCTION

The numerical solution of unilateral problems specific to structural elements which feature receding and advancing contact has been reported by Mahmoud et al. [l]. They provide solutions for beams and plates based upon linear function theory and a rather elaborate algorithm for activation and deactivation of contact constraints. This algorithm requires computation of an ‘interpenetration’ factor for detection of contact, and a data structure which records the history of contact events and their type (receding or advancing). The scheme has evolved from the theory given by Dundurs [2]. Upon further reflection and in teaching nonlinear mechanics, it became evident that in situations where contact can be modeled through discrete, independent (compressive) supports, a considerably simpler scheme is possible. The result is a single parameter solution for contacts and separations with two switch variables; one denotes position of the constraint surface with respect to the structure, whether above or below it, and the other denotes the constraint state, either contact or separation. Others have applied variations of this idea, but none in so simple a manner. This paper provides the theory, an algorithm and a succinct FORTRAN subroutine which, with minor modifications to a standard linear code, makes the solution of formidable contact problems elementary. Indeed, the concepts are appropriate for an undergraduate mechanics class. Results from sample problems tally precisely with those given previously [I].

THEORY

Consider a linear elastic structure under a load applied along R and constrained by permanent boundary conditions along P and unilateral contact conditions along C. The latter are unknown a priori in the sense that along this boundary, initial prescripcA s

m*-hi

tions may change, contacts may separate and vice versa, dependent on the load scenario and structural response to it. We call such separations or contacts contact events. Within this framework, over a typical increment, consider an applied load {R*} which yields a structural response via the stiffness equation [K] {tv *} = {R *} where {w *} contains the deflections d* and rotations and [K] is the usual stiffness matrix, but updated to reflect previous contact events. Clearly the next event (or non-event) will be incipient for some scale factor s such that [K](w) =s{R*)

for

0
(1)

because the problem is linear over each increment. We may refer to {R *} as a test load and {w *} as the test response, but after each increment it is preferable to refer to subsequent values of {R*} as the residual applied load. The response for an incipient event is alternatively obtained from {w} = s{w*}

(2)

and the amount of load consumed by the process is {R } = s { R *}. Hence the new residual load becomes {R*} -{R} = {R*}(l -s). For a non-event, the entire residual load may be applied within that increment and s = 1. Axiom. If a contact event is to occur under an applied load s{R*}, then 0 es < 1. It remains to establish s as the single parameter to determine contact events and, importantly, to do this solely kinematically. A kinematic relation for s in terms of geometry and response is obtained as follows. At any one candidate point on C, let c be the absolute position of the constraint, D be the current position of the structure, d* be the deflection response to residual load {R*} measured from the unloaded reference position of the structure and n = f 1, + for a constraint above the structure and - for one below it. (A typical

509

510

N.

J.

SALAMON

Y Undeformed

ref.

constraints (n=-1).

Fig. I. Typical response d*(x) to a test load showing interpenetrations of spring constraints; solutions D@“(s) at least satisfy the constraints at position k and D(*)(x) is the one corresponding to s,,, where D”‘x = P(x) + D(x) and &‘(x) = aid*. configuration is illustrated in Fig. 1.) Then, upon application of load, we can define the quantities

C4n(D+d*-C)9

-s

(3)

where C and S are measures of interpentration between bodies and extension across the contact interface respectively. Alternatively, C > 0 or S > 0 is a precursor for contact or separation. (In [l] C > 0 serves to detect contact.) At the incipient occurrence of either, d* + s .d* and C = 5’ = 0 so that (3) yields s = (c - D)/d*.

(4)

The upper and lower surface variable n is irrelevant to calculation of s. In kinematic terms, over what range of s is an event possible? If we define the gap (non-negative) existing at some point between the current position of the structure and an inactive constraint as gap 4 n(c - D) = nsd* > 0

(5)

where the equality assumes an increment of load has been applied and s determined, then for a contact event we have C > 0 so that 0~Cxgap=n2(d*)2s(l-s)~O~s~1.

(6)

Following a similar argument, the same range for s is found for a separation event. Any other value for s implies no event is possible. Lemma. For a contact or separation event to occur, O 1, there is no change in contact status. Since this range in s detects either a subsequent contact or separation, to distinguish between them, a current contact status flag must be set for each constraint. Then s in [0, l] signals the opposite event. It is noted that s = 0 is a trivial case since no further evolution of the problem is then possible.

For multiple candidate constraints, several possible events may be detected. Over an increment in load, the correct solution is the one compatible with (does not violate) the current set of constraints; it generates neither interpenetration nor interface tension (interface extension). If we first construct a matrix of possible valid incremental solutions d at any one constraint, i.e. solutions at upper or lower constraints with respect to their current status-active or inactive--then two inequality conditions arise: (1) d 3 c - D E sd* for upper, active or lower, inactive constraints, (2) d 6 c - D = sd* for lower, active or upper, inactive constraints. But for condition 1, c - D < 0 and the range of interest in d and d* for occurrence of an incipient event is d < 0 and d* < 0; hence multiplication through by -1 reduces condition 1 to condition 2 and d/d* 6 s for a valid solution where s is in [0, 11. In general, for NC total constraints along C we have d@‘(x) & skd*(x),

0 6 sk 6 1, k=1,2,...,NC

(7)

where solutions d(“) are valid at least at constraint position k and x is an undeformed reference coordinate along the structure (see Fig. 1). In light of (7), the correct solution d(x) must satisfy constraint conditions at all points x,, i.e.

4x,,) =.d(x)/d*(x)

= Min{s,, s2. .

, sNc}.

(9)

R~eding/advancing contact of structural elements simplified Lemma. For all possible events such that 0 6 s 6 1, the correct solution and choice of subsequent constraint action is the one for which s is the minimum. Hence the a~umulated solution sebsequent to the current increment is

O(X)+

D(X)

+

d(X)

s

D(X)

+

where smiffis given by the ~ght-hand each s is computed from (4). ALGORITHM

S,ind*(X)

(10)

side of (9) and

DEVELOPMENT

Apparently the algorithm follows directly from the theory: solve for w*, extract d* and compute s for each constraint, determine s,,,~,,and its location, set the event, update the response and continue. This is true except for the case of receding contact [2] which we prefer to call stationary contact because we view the contact event process as follows: advancing * the activation (connection) of additional constraints, receding * the deactivation (separation) of additional constraints, and stationary *no change in contact status under change in load. For the case of stationary contact, D = c initially; most often D = c = 0. Then from eqn (4), s = 0 for all constraints making it impossible to choose the succeeding event. (Alternatively, receding contact can be viewed as a case where the gaps, given by (S), are zero.) A bright idea is to set the gaps not to zero, but to a very small number and treat the receding problem the same as an advancing contact problem. But then in practice one encounters precison problems due to a loss of numerical significance. For robustness it would be nice to augment the theory with an algorithm which can accept gaps either zero or not. To do this, Dundurs’ proof[2] that stationary contact must precede advancing contact is employed. First we define a stationary constraint as one with a gap < 6 where c is a small prescribed value. Second we evaluate a scale factor check-variable .rchkat each stationary constraint such that . * ‘%hk =

[sep G”,f i

‘ii*

(11)

where iEep= + 1, + for separation, - for contact, each taken currently, irurr= n as defined under the theory section and d* provides a value of relative merit, i.e. the larger d*, the greater the interpenetration or interface extension, and the smaller schk. But (11) Only holds for &,k > 0; hence we finally compute the scale factor for a stationary constraint as if

S&k> 0,

then

S =

S&k8

fl2)

511

where 6 is a very small value to insure (i) that in competition with scale factors for advancing or receding contact constraints, stationary values will be minimum and win the contest and (ii) any contribution to accumulated D will be nil (beyond computing machine significance). The latter is necessary because the numerical solution to a stationary problem is iterative until an equilibrium configuration is found, then D is set (see [I]). Notably, despite this nuance, scale factors for stationary, advancing and receding constraints compete directly and none require further special treatment such as a data structure to distinguish between them. The algorithm then becomes: 1. 2.

Set For 3. 4. 5.

smin= 1 and Location = 0 each contact candidate position K, do Compute the gap: eqn (5) ffd*=O, thens-I If gap < 6, the event is stationary, then 6. s = 1 and 7. Compute s,,,: eqn (11) 8. If schk> 0, compute s: eqn (12) 9. Else 10. Compute s: eqn (4), the event is not stationary 11. End If 12. If 0 < s < Smintthen s,;, = s at Location = K 13. End do 14. ;ftLo;;~e~,O and isfz;mmtz}, then set it

I in contact I Thus the algorithm supplies the new scale factor and event sites to the main code for further processing. RESULTS AND CONCLUSION

The result is a strikingly simple computer program. For instance, a FORTRAN subroutine CONTCT is given in the Appendix. When it was appended to a conventional finite element code modified to accept the update information from SUBROUTINE CONTCT, example problems from [I] were solved and the results are precisely the same. (We note a typographical error in [1], p. 635, Plate over settled support: K = 1390 should be K = 11,120.) The Appendix requires some explanation. The variables L>,c and s are doubly repeated as DD, CC and SS respectively. Other notation is defined in the code or above. One departure from the above algorithm is that the code permits multiple events per load increment; if the position status of a candidate pair is less than a prescribed small value, it is accepted as incipient and action is taken even if its corresponding value of s is not the minimum. We conclude that this class of problem is now elementary. The theory is clearly founded, the algorithm follows directly and the code is concise and easily understood. Excess baggage associated with

512

N. J.

data structures for the history of contact and calculation of interpenetration or interface extension is discarded. The entire procedure is a function of the scale factor and on/off and above/below ‘soft-toggle’ switches.

SALAMON REFERENCES

1. F. F. Mahmoud, N. J. Salamon and T. P. Pawlak, Simulation of structural elements in receding/advancing contact. Comnut. Struct. 22. 629-635 (1986). 2. J. Dundurs, Properties of elastic bodies incontact. In

Acknowledgement-The authorisgrateful to T. P. Pawlak who helped debug an early version ofthesubroutineintothe

The Mechanics of Contact Between Deformable Bodies

(Editedby J. J. Kalkerand A. D. dePater). Delft University Press(1975).

smallhoursof one morning. APPENDIX

SUBROUTINE CONTCT C C C C C

C C C

C C C

C C C c C

(DD,

DTEST, CC, SS,

SMIN,

ID, ISEP, ISURF, NC)

****** DETECTION AND SETTING CONTACT EVENTS NJS 28 AUG 87 ***t** Notation: DD, DTEST, CC, SS are D, d*, c, s in the text respectively ID(I) identifies the node on the structure associated with contact event site I NDOFPN = Number of degrees-of-freedom (dof) per node INDOF sets the degree-of-freedom of interest, e.g., in the present application, the number associated with deflection DELTA = very small value (see Eq. 12 in the text) EPSGAP. EPSMUL set small cut-off values for 'zero' qap and admittance.of multiple events per load increment, respectively Dimension DD and DTEST at least to total system dof Arrays: Dimension CC, SS, ID, ISEP and ISURF at least to the total number of constraints NC DD, DTEST, CC, ID, ISEP, ISURF and NC Supply: Returns: SS and SMIN

10

20

30 40 50 60

COMMON /DOF/ INDOF,NDOFPN DIMENSION DD(l), DTEST(l), CC(l), SS(l), ID(l), ISEP(l), ISLJRF(1) DATA DELTA/l.E-20/, EPSGAP/l.E-lO/, EPSMUL/.390625E-02/ SMIN = 1.0 DO 10 I = 1, NC NODE = ID(I) NDOF = fNODE-1) * NDOFPN + INDOF GAP = CC(I) - DD(NDOF) IF (DTEST(NDOF).EQ.O.) THEN SS(1) = 1.0 ELSE IF (ABS(GAP).LT.EPSGAP) THEN SS(1) = 1.0 SSCHK = ISEP(1) * ISURF( I) / DTEST(NDOF) IF (SSCHK.GT.0.), SS(1) . . = SSCHK * DELTA ELSE . SS(I) = GAP / DTEST(NDOF ENDIF IF fSS(I).GT.O. .AND. SS(1) .LT.SMIN) SMIN=SS(I) CONTINUEIF (SMIN.EQ.l.) THEN WRITE(26,20) SMIN SCALE FACTOR = I, E10.3) FORMAT(ZX, 'NO EVENTS. RETURNING. RETURN END IF DO 30 I = 1, NC DIFF = ABS(SS(1) - SMIN) / SMIN IF (DIFF.LT.EPSMUL) THEN NODE = ID(I) IF (ISEP(I).EQ.l) THEN ISEP(1) = -1 WRITE(26,40) ELSE ISEP(I) = 1 WRITE126.50) END IF . ' WRITE(26,60) I, NODE, ISURF( SMIN END IF CONTINUE FORMAT(2X, a****** MAKING CONTACT ********I) FORMAT(2X, I****** BREAKING CONTACT *******I) FORMAT(3Xr 'CONTACT PAIR: ',13,' NODE: ',13, ' SURFACE: ',13, ' SCALE FACTOR = ',E11.4) & RETURN END