Cellular automata model for a contact problem

MATHEMATICAL AND COMPUTER MODELLING PERGAMON

Mathematical

and Computer

Modelling

36 (2002)

1099-1114 www.elsevier.com/locate/mcrn

Cellular Automata Model for a Contact Problem M. ABDELLAOUI AND A. EL JAI Systems Theory Laboratory, University of Perpignan 52 Av. Villeneuve, 66860 Perpignan, France @univ-perp.fr M. SHILLOR Department of Mathematics and Statistics Oakland University, Rochester, MI 48309, U.S.A. shillorQoakland.edu (Received

and

accepted

May

2002)

Abstract-The aim of this paper is to apply the cellular automata approach to the description of deformations of an elastic body and the frictionless contact between the body and a rigid foundation The model is based on a microscopic approach and can be easily implemented. A numerical algorlthnl is constructed and simulation results are presented, illustrating the properties of the model and the behavior of its solutions. @ 2002 Elsevier Science Ltd. All rights reserved.

Keywords-cellular

automata, Elastic body, Frictionless

contact

1. INTRODUCTION Models of problems in contact mechanics are usually formulated in terms of partial differential equations (PDE) or inequalities related to macroscopic description of the phenomena. The UIIderlying mathematical structure is quite complex and remains as an unfamiliar jargon for most applied users of such models, such as many design engineers. An alternative approach to describe contact processes consists of a cellular automata model which may be much simpler and morr convenient for numerical implementation. Cellular automata approach (see, e.g., [l]) can be used, and has been used to model various physical phenomena [2-51. Whereas models employing partial differential equations are based on a global or macroscopic approach to the system, cellular automata models are deduced from microscopic description of the behaviour of the particles constituting the system. Both PDE and cellular automata models describe the systems’ evolution in space and time. however, the former assume that all the variables describing the process are continuous (often more regular or smooth), while the latter is based on a discrete approach to all of the components. The development and use of cellular automata models has increased considerably with the availability , of powerful computers.

0895-7177/02/$ - see front matter @ 2002 Elsevier PII: SO895-7177(02)00261-3

Science Ltd. All rights reserved.

Typeset

by d&-T@

1100

M. ABDGLLAOUIet al.

The aim of this paper is to explore the cellular a u t o m a t a approach to model a two-dimensional process of frictionless contact between an elastic body and a rigid foundation. We consider a microscopic description of the phenomena by studying the behaviour of the body's particles, We describe the cellular a u t o m a t a approach in Section 2. A cellular a u t o m a t a model of an elastic body is constructed in Section 3 and numerical simulations are presented in Section 4. A cellular a u t o m a t a model for frictionless contact is presented in Section 5, together with numerical simulations on the evolution of the contact. The simulations illustrate this approach, and are based on an adaptation of Matlab code to the model setting. In the Appendix, we provide a very short description of a P D E model for the contact problem for the benefit of the interested readers.

2. C E L L U L A R

AUTOMATA

APPROACH

The spatio-temporal evolution of dynamical systems is often governed by nonlinear partial differential equations. Due to the nonlinearities, the solutions of these dynamical equations can be very complex and are often very sensitive to initial conditions, leading to chaotic behavior. The cellular a u t o m a t a (CA) approach offers an alternative way to describe and simulate the behavior of these complex systems. As defined by Wolfram [1], CA are mathematical modelling tools for complex dynamical systems where time, space, and states are discrete variables. A CA consists of a regular space, a local neighborhood, a set of admissible states, and a set of evolution rules and is denoted by the quadruplet A = (£, £, 12, y ) where we have the following. • £ is a regular lattice which consists of identical elements called cells of a spatial domain fL Each cell is characterized by a state which takes possible values in a finite set of states updated according to the transition rules. • £ is a finite state set containing the values which may be taken by the cell's state at each time step, when the system evolves from an initial state to a final state according to the transition rules. For the problem considered here, the states represent the displacements of the particles constituting the elastic body. • 12(c) is the neighborhood containing all cells which have an influence on the central ceil (c). In a two-dimensional lattice, the neighborhoods often used are the v o n N e u m a n n , the Moore, and the uniform ones. Other types of neighborhoods can be considered depending on the system nature. The size of the neighborhood is the number n of cells in 12(c) which interact which e. We assume that n is constant (nonvarying in time). • 5r is the transition rule, also called system dynamics, which governs the system evolution and determines how cells change their states, in time, depending on their neighbors. The transition rule may be deterministic or probabilistic and can be given as a mapping defined by g~ ---* g, Y: et(12(C))

--'

gt+l (C),

where the state of the cell (c) at time t + 1, et+l(C), is determined from the state of the neighborhood cells in 12(c) at time t, et(12(c)). S~,,

i

/ /

(a)

(b) Figure 1. Three main types of two-dimensional neighbourhoods.

(c)

Cellular Automata Model

1101

In practice, when simulating CA rules, one cannot consider an infinite domain, thus, the lattice is finite and has a boundary. The state equation must be augmented with an initial state and boundary conditions. The initial condition may be given or randomly generated, and for the boundary conditions, we can refer to [6] where many cases are presented. 3.

CELLULAR

FOR

AUTOMATA

AN

ELASTIC

MODEL

BODY

In a continuous medium, there exists a nominal configuration Co of the body, that is to say, a certain distribution of the particles constituting the medium. When exciting forces are applied, this configuration may change to a new one Ct, considered at time t, t > to when the forces are applied. The system is elastic if removing the forces allows the body to return to its initial configuration. The body is plastic if it remains in its new configuration after the forces are removed. As the body changes from the initial configuration Co to the current one Ct, we ---4 ) consider a displacement vector U = PoPi where P1 is the location at time t of the point P0 of Co. This simple observation allows for the description of the system by a cellular automata model. The CA model describes the possible displacements of the particles of the body under the action of different forces. We simplify the system by considering the medium to be perfect, so there are no gaps between the particles, we assume the problem is two dimensional, and the medium is isothermal, so the heat conduction problem is neglected. In the problem, we consider the displacements are due to the global deformation resulting from the forces applied to the system. The cellular automata dynamics is deduced from the mechanism which generates the movement of the particles. The model is then developed from a microscopic point of view. Vre consider a two-dimensional solid represented by a lattice of particles. An example of such solid structure can be found in various applications such as crystals, ceramics, and organic polymerase. When the particle arrangement is in an equilibrium state, we assume that the distance between the particles is equal to a certain value do (see Figure 3). In the particle distribution above, each particle has four neighbors (east, west, south, and north). \\\\

Figure 2. Displacement of a particle in an elastic body. i

\j

\j

\j

\j

k_J

() () () () () () () () G) d° A-~ r'~

Ah

Ah

f'~

Figure 3. Particles constituting the body.

1102

IVl. ABDELLAOUI et al.

3a. T h e F o r c e s The forces acting on an elastic body generate displacement of the particles, and this yields a new configuration. For small displacements, the force f may be assumed to be proportional to the variation of the distance between the particles and is parallel to the direction of the displacement. Let (x, Y) = (z(i, j, t), y(i,j, t)) be the position coordinates of the particle at the grid point (i, j) at time t, and let Ax(i, j, t) be its displacement from its initial configuration. Then,

f(i,j, t) = CAx(i,j, t).

(3.1)

Here, C is a constant which depends on the links energy. Each particle is submitted to the forces due to its four neighbors. Thus, we can write the following relations (see Figure 4): f l (i,

j, t) =

C A x 1 (i, j, t) ~- C [x(i - 1, j, t) - x ( i , j, t)], - 1,t) - x ( i , j , t ) ] ,

f2(i,j,t)

= CAx2(i,j,t)

= C[x(i,j

fa(i,j,t)

= CAxa(i,j,t)

= C [x(i + 1 , j , t )

f4 (i, j, t) = C A x 4 ( i ,

- x(i,j,t)],

j, t) = C [x(i, j + 1, t) - x(i, j, t)].

4

-

f4 3

1

+

f

2

Figure 4. T h e forces applied by t h e four neighbors.

The fbrce distribution can be written in the following matrix form:

fi(i,j,t)) f~(i,j,t) f3(i,j,t) f4 (i, j, t)

[ Axl(i,j,t) ~ |Ax2(i,j,t)| c i/,x3(i,j,t)

(3.2)

!

\ Ax4(i, j, t) /

3b. Power and P o t e n t i a l E n e r g y Let f be the resultant force acting on a particle. The power generated by f is given by P ~ fv~

where the velocity v is V~---

Ax At"

On the other hand, P may be obtained from the potential energy

P = AEp At Thus,

AEp = PAt = fvAt = f a x , since Ax =

vAt.

Ep

as

Cellular A u t o m a t a Model

1103

The potential energy of a particle depends only on the initial and the, final positions occupied by the particle, and is given by Ep = f A x . (3.3) A system of n particles has the potential energy Eptot = ~

(3.4)

fiAx~.

i=I

3c. P a r t i c l e D i s p l a c e m e n t When a force is applied to the body, it undergoes a deformation and the particles' positions change. To describe these position changes, we let the particle in the central location be represented by a full circle, and its four neighbors by empty circles. We first describe the displacements in the x-direction. A similar argument applies to the motion in the g-direction. The central particle moves under the effect of the x-components of the forces applied by its neighbors, denoted by fl, f'2, f3, and J~4. These lead to the displacements '//~1, U2, U::;, alld LL1 along th(~ OX axis. as depicted in Figure 5. Here, the sign (+) denotes the location of the mass (:enter of the white particles. The movement of the central particle depends on the displacements of its four white neighbors. Let r,,j be the position of the central particle

and the neighboring particles are at the positions ri - 1 ,j ~

7"i+ 1 ,j ~

?'i,j - 1 ~

r i ,j + 1

respectively. We can define the displacements of the neighboring partMes along the x-axis by" U l ( / , . j , ~ ) = X i - - l , j ( t ) --

Xi,j(t)

:

z~Xl,

~2(i,j,t)

= x~,j_l(t)

- x~,~(t) = A x 2 ,

u3(i,j,t)

= x~+l,o(t) - xi,j(t) = Ax3,

ua(i,j,t)

= xi,j+l(t) - xw(t)

= Ax4.

U4

()~

f,

~i )

U3

\,\ l

3

~b

f3

+

f

u~ U 2

Figure 5. Displacements of the central particle under its neighbors' forces.

Similar relations hold along the OY axis. In the same way, we can define the displacements Vl, v2, v3, v4 of the neighboring particles along the OY axis. If we choose the origin of the coordinates system for the four neighboring particles as the position of the central particle, (r~ d (t) = 0), then the mass center of the particle system is now at the point (XcM, YCM), where, 1 X C M = ~ [Ul q- ?*2 1- ~3 q- ~/'4],

1 YCIvI = ~ [Vl n-'~'2 @'t'3 + U41.

1104

M. ABDELLAOUI et al.

Due to the symmetry with respect to the mass center, the new coordinate of the central particle is then

xi,j(t + 1)

= 2XcM

1 = ~ [ul + u2 + u3 + u4].

The motion of black particle results in a new distribution of the displacement u(i,j). Since at time t + 1, the white particles will move and the black will be stationary, the displacement u(i,j, t + 1) now has to be computed with respect to each white particle location. For instance, ul (i + 1, j, t + 1) will be the displacement of the black particle as seen by the particle located at 7*i+ l ,j .

Indeed, the displacement of the central particle is as follows: ul(i

+

1,j,t + 1)

= 2XcM

-- X i + l , j

1

= ~ [ul + u 2 + u a + u 4 ] - us

1 = ~ [ul + u2 - us + u4].

In the same way, we can compute u2(t + 1), u3(t + 1), u4(t + 1), and this leads to the following evolution law for any central particle:

ul(i + 1,j,t + 1) /

(ul(i,j,t)

u2(i,j + 1,l + 11

I! u2(i, , l I

u 3 ( i - 1,j,t + 1)

= W [u3(i'j't) I

u4(i,j-- l,t + l)

\u4(i,j,t) /

(3.5)

where W is the propagation matrix given by

1 w

= [

1

1

1

1

1

1

1

1

1

1

-1

1

-1

-1

1/ -1 '

3d. T h e S t r a i n When a particle is in the neighborhood of a position r0, the particle's links may be considered as spring-like with stiffness K (within linear elasticity). We can evaluate the elasticity modulus E with respect t o / ( and r0 as follows. The restoring force is K ( r - ro), while the stress ~r satisfies cr = E(r - ro)/ro, and thus, K 1"0

It is well known that, for an elastic material, when compressed, the displacement is proportional to the force over a wide range of compressive forces, while under tensile force, the displacement is proportional to the force up to a certain value called the decohesion value. Once this level is reached, experiments show that there are nonlinearities that take place, such as plastic deformation, cracks may appear and microscopically dislocations in the material take place (see Figure 6).

Cellular Automata Model

1105

Stiffness in the neighborhoodof the balance

F maxi

t>°n/

e

r D ~-k

Value of decohesion

Value of balance Figure 6. Equilibrium graphs.

F

f

!

(),.()

IF

~

f

•

•

Figure 7. Lengthening of a polymer submitted to a traction force.

We consider a regular arrangement of the particles, without any body forces, but acted upon by an applied positive traction along the axis OX. Thus, a polymer of initial length g (see Figure 7) will be lengthened by Ag. Its relative lengthening is then Ag g The ratio Ag/g represents the strain of the particles, and we may rewrite ei using the displacement u~, AUl , with Aux = 2XCM -- X,+l,j = u l ( i + l , j , t + l ) , tt 1

thus, ¢1 becomes ¢1

~"

ul(i + 1 , j , t + 1) ul

~

1 2ul [Ul + u2 - ua + u41

"

Analogously, we can compute e2, ~3, ¢4. Finally, the strain law will be described by

¢2(i,j + l , t + l)

ca(i

= W'

i ~(~,j,t)

1,j,t + 1)

/ ua(i,j,t)

¢4(i,j - 1,t + 1)

\u4(~,j,t)

-

,

(a.o)

1106

M. ABDELLAOUI et al.

where IV' is given by

I,'V' = 1 2

1

1

1

1

\

Ul

Ul

Ul

Ul

1

1

1

1

~2

~2

~2

U2

1

1

1

1

1~3

~3

~3

~3

1

1

1

1

~4

~4

U4

~4

3e. T h e S t r e s s Let S be the lattice orthogonal section to the traction axis, then, when the body evolves, the strain is proportional to the ratio f O- = ~ 1

where cr is the longitudinal stress of the lattice (the units are N / m 2 ) . known as Hooke law and may be expressed by

The obtained relation is

~:Ec. E is the longitudinal elasticity modulus. Furthermore, if we consider ~1 as a function of u~, then we obtain

ch = E e l = E u l ( i + l ' j ' t + l )

1

ttl

= 2~lE[u1

4- u 2 - u 3 ~- u 4 ] .

The calculation of or2, ~3, cr4 may be done in the same way that for or1. The evolution taw of the stress is then

~2(i,j + l , t + l) ~a(i-l,j,t+l) (~4(i,j - 1, t + 1) where W " is given by

W"=

_E 1 2

t) = W " | u 2 ( i , j , t) [ u a ( i , j , t) \ u 4 ( i , j , t)

1

1

~1

Ul

1

1

1

u2

1

1

Ul

Ul

(3.7)

'

1

u2

u2

i

i

1

I

u2

U3

U3

U3

U3

1

1

1

1

~4

U4

U4

~4

Of course, there will be no displacement on the boundary F1. 4.

NUMERICAL

APPROACH

Let us define the different components of the cellular automata model. 4a. State

Set

We must define the state space which will be used to introduce the transitions rules for the displacement of the particles. We consider a set of five states defined by g={e~, 0