10626

FACTA UNIVERSITATIS

Series: Mechanical Engineering

https://doi.org/10.22190/FUME220327020F

Original scientific paper

NUMERICAL SIMULATION OF DYNAMICS OF BLOCK MEDIA

BY MOVABLE LATTICE AND MOVABLE AUTOMATA

METHODS

Alexander E. Filippov1,2, Valentin L. Popov2

1Donetsk Institute for Physics and Engineering, Donetsk, Ukraine
2Technische Universität Berlin, Berlin, Germany

Abstract. Two versions of modified Burridge-Knopoff model including state dependent

friction, elastic force and thermal conductivity are derived. The friction model

describes a velocity weakening of friction and elasticity between moving blocks and an

increase of both static friction and rigidity during stick periods as well their weakening

during motion. It provides a simplified but qualitatively correct behavior including the

transition from smooth sliding to stick-slip behavior, which is often observed in various

tribological and tectonic systems. Attractor properties of the model dynamics is studied

also. The alternative versions of the model are proposed which apply a simulation of

the motion of interacting elastically connected mesh elements and motion of relatively

large solid blocks, utilizing technique of the movable cellular automata. First version

of the model was already basically studied before. Its advanced version here involves

all components of the real system: state-depending friction and changeable rigidity, as

well as heat production and thermal conductivity. Model based on the movable

automata also involves the components included into traditional lattice model. It has

its own ad-vantages and disadvantages which are also discussed in the paper.

Key words: Burridge-Knopoff model, Friction, Tectonics, Voronoi diagram, Delaunay

triangulation

1. INTRODUCTION

One of the most important applications in the study of the behavioral patterns of

mechanical systems of block structures are the geological faults and block environments.

Of particular interest, of course, is the study of the possible variations of their

deformation, potential prediction or even prevention of the strong seismic shocks. The

Received: March 27, 2022 / Accepted April 23, 2022

Corresponding author: Alexander E. Filippov

Donetsk Institute for Physics and Engineering, R. Luxemburg Str. 72, 83114 Donetsk, Ukraine

E-mail: filippov_ae@yahoo.com

2 A.E. FILIPPOV, V.L. POPOV

fact that the statistics of earthquake force and their time correlations meet the laws of

Gutenberg-Richter [1,2] and Omori [1,3], typical for self-organized critical systems [4,5],

is often the basis for the conclusion that earthquakes are a process formed on all spatial

scales ranging from microscopic to continental plate scale. During last decades [6,10] on

the basis of modeling by movable cellular automatons and full-scale experiments it was

shown some possibility of initiation of the fault wing displacement and, consequently, the

release of a part of the accumulated elastic energy as a result of local effects on the fault

area, including vibration load and watering [11-16].

One of the often applied to study the problems under interest is the well-known

Burridge-Knopoff (BK) model [17,18] initially proposed many years ago to investigate

statistical properties of earthquakes. Numerical studies by Carlson et al. [19,20] have

demonstrated that BK model can reproduce characteristic empirical features of tectonic

processes such as the Gutenberg-Richter law for the magnitude distribution of

earthquakes, or the Omori law for statistics of aftershocks [19-22], both properties

stemming from the so called “self-organized criticality” of this system. It has been

intensively used to simulate different aspects of the problem [21-40] and to discuss

general properties of earthquakes statistics as well as predictability of earthquakes.

Numerical simulations give evidences that the self-organized criticality and the

corresponding fractal attractor of the system is closely related to dynamic structures with

“traveling waves” [23], their ordering and specific “phase transitions” [24] controlled by a

number of parameters (external driving velocity, springs stiffness, number of blocks, their

mutual interaction and so on). It is in particular the dependence of the dynamic properties

of the BK model on spring stiffness which makes it necessary to introduce the

generalization of the friction law proposed in the works [41,17].

The physical reason for the stick-slip instability in the Burridge-Knopoff model is the

assumed decrease of the friction force with the sliding velocity [19,20]. Motion of a single

block with this friction law is always unstable which does not correspond to properties of

real tribological systems. The real law of rock friction is more complicated [29,34]. In

this article we base on the realistic friction laws described in the handbooks [34,35]. The

main qualitative pictures of a realistic law are: (a) approximately logarithmic increase of

the static friction force as a function of contact time – the property, found already by

Coulomb, and (b) a logarithmic dependence of the sliding friction on the sliding velocity.

Both properties can be described in the frame of “state dependent” friction laws by

introducing additional internal variables describing the state of the contact [17,41-43] as

well as shear band propagation effect on the deformation and fracture [44,45], generation

and propagation of slow deformation in geomedia [46-48].

Modified version of the BK model with a state dependent friction force, reproduces in

the simplest way both the velocity weakening friction and the increase of static friction

with time when the block is not moving.

2. METHOD

Let us briefly remind the original BK model. The blocks of mass m are attached to a

moving surface by springs with stiffness k1 and are coupled to each other by springs with

stiffness k2. The moving surface has velocity v, and the blocks are in contact with a rough

Block-media Motion Lattice and Movable Automata Numerical Models 3

substrate. The friction force F between the blocks and the rough surface is assumed to

depend only on velocity. The equations of the BK model can be written in the following

form:

   
2

2 1 1 12
2 ( )

j j j j j

u
m k u u u k vt u F v

t
 


     


, (1)

where v=const is the external driving velocity and vj≡∂uj/∂t is an array of individual block

velocities (j=1, … N). The sliding friction force is supposed to decrease monotonously

from a constant initial value F0. It is further supposed that the static friction F(vj→0) can

possess any necessary negative value to prevent back sliding

0
, 1 ; / 0

1 2 / (1 )( )

( , 0] / 0

F
F u t

vF v

u t


       

    


  . (2)

The parameter α defines a rate of friction decrease when block starts to slide, and σ is the

acceleration of a block at instant when slipping begins [19,20]. The friction equation (2)

is an oversimplification of real properties of static and kinetic friction which is now

explained in detail in the books [34,35]. Dieterich has proposed friction equations with

internal variables [29,30], where the friction force is dependent upon an additional

variable that describes the state of the contact zone. This variable is in a sense a “melting

parameter.” Friction force at non-zero velocity drops down from its initial value due to a

“shear melting” effect which may have different physical nature [17-20]. When the

motion stops, the surfaces start to form new bonds and the static friction increases with

time.

The dynamics of systems with state dependent friction has been investigated in a

number of papers [29-37]. Majority of these studies were devoted to the simple one-

particle version of the model, but it was also done for many-body BK model with a state

dependent friction [17], where the friction force is defined by the additional kinetic

equation:

 1 0 2
[ ( )]

j j

F v t
F F v


  


  . (3)

with β2<0, at vj>0, and F(vj)=-∞ in opposite case vj≤0. The parameters β1 and β2 have

following meaning. When a block starts to slide with vj>0 its friction force monotonously

decreases from an initial value F0. General time scale of this process is determined by the

first parameter β1 and an effectiveness of the melting is given by a relation between

absolute value of negative parameter β2<0 and β1. It should be noted that the kinetic

equations for state dependent friction at sliding are not unique. Different forms of this

equation have been proposed [17,18] leading to qualitatively similar behavior of the

dynamic friction at vj>0. The above form is just simplest and even linear one.

The equations of motion in Eqs. (1-3) are nonlinear wave equations. Last few decades

different variants of such equations have been widely studied, starting from the very early

implementation of the numerical simulations [38]. In particular, it has been shown that N

interacting segments of the nonlinear chain form a collective attractor with energy transfer

4 A.E. FILIPPOV, V.L. POPOV

performed by nonlinear excitations [39,40]. Depending on the relations between total

energy and strength of the interaction between the blocks k2, the chain can form nearly

uniform or strongly non-uniform structures and phase patterns.

Analogous behavior should be expected in the modified Burridge-Knopoff model as

well. Instantaneously moving blocks must be involved with stationary ones in the overall

motion. This causes a “detachment” wave, propagating along the chain. Such waves were

found and studied in the original BK model [18], and the analogous process exists in the

modified one. There are two possible types of traveling waves: (a) after a transition time

all the blocks are provoked to move simultaneously, (b) in steady state some of the blocks

can be found instantaneously motionless while others move. One would expect the

reaction of the chain to depend on the strength of the interaction between the blocks. This

question has been studied for the original BK model in [20], and a specific “phase

transition” between correlated and uncorrelated behavior has been found.

Real contact of two surfaces is two-dimensional. Let us generalize the MBK to the

case. The generalized model is very similar to the one-dimensional model but

incorporates a 2D array of blocks connected by elastic springs in both directions. All

other components of the MBK remain unchanged.

The system of equations of motion takes the form:

   
2

, ,

2 1, 1, , 1 , 1 , 1 , ,2
4 ( );

j n j n

j n j n j n j n j n j n j n

u u
m k u u u u u k vt u F v

t t
   

 
        

 
 (4)

 
, ,

1 0 , 2 , ,

, ,

( ( ))
, 0, 0

( ) 0

j n j n

j n j n j n

j n j n

F v t
F F v v

F v v

  


    


  

. (5)

where j=1, …, Nx and n=1, …, Ny. Nx and Ny are the numbers of elements in the x- and y-

directions. For general applicability to tribological problems Eq.(4) contains also a

viscous term η∂uj/∂t. Our calculations show that attractor properties are weakly influenced

by the viscous term at sufficiently high velocities vj. For generality, we keep small

nonzero value η=0.05. Other parameters are: β1=1, β2=125, F0=1, k1=1. k2=36, v=0.2.

Large scale structure of the wave state in the two-dimensional model with number of

dimensionless grid cells (Nx = 700, Ny = 300) shown by the velocities v=v(x,y;t)≡∂u/∂t is

reproduced in Fig. 1 by the colored surface (MatLab color-map ‘jet’). Dark blue

corresponds to zero. Mutual scattering manifests itself in high sharp peaks of the “events”

reproduced here by the intensive red of the corresponding densities. Instant structure of

the wave state and history of the events are combined in Fig. 2. The subplots (a) and (c)

here represent snapshots of the instantaneous densities of the local displacements

u=u(x,y;t) and velocities v=v(x,y;t)≡∂u/∂t, respectively. Time-space representation of this

process is shown in the subplot (b), and time dependency of total area of events (d). The

information about attractor behavior is accumulated in Fig. 3 and briefly described in the

caption to this figure.

Block-media Motion Lattice and Movable Automata Numerical Models 5

Fig. 1 Large scale structure of the wave state in the two-dimensional model shown by the

velocities v=v(x,y;t)≡∂u/∂t reproduced by the colored surface (MatLab color-map ‘jet’).

Dark blue corresponds to zero. Mutual scattering manifests itself in high sharp peaks of

the “events” reproduced here by the intensive red of the corresponding densities.

Fig. 2 Instant structure of the wave state and history of the events in the model. Subplots

(a) and (c) represent snapshots of the instantaneous densities of the local displacements

u=u(x,y;t) and velocities v=v(x,y;t)≡∂u/∂t, respectively. The colormaps are equidistant and

normalized to the intervals between maximum (red) and minimum (blue) values for the

displacements [-1,1] and velocities [0 1.5]. Time-space representation of this process is

shown in the (b) and (d) by the spatial-temporal map and time dependency of total area of

events.

6 A.E. FILIPPOV, V.L. POPOV

Fig. 3 The attractor of the variable effective friction in projection to the plane {v, Ffriction},

where <Ffriction >=<F[vj(t)]+ ηvj(t)> is shown by color-map ‘jet’ in the subplot (a). The

attractor behavior of the system causes quasi-periodic behavior of the mean area of the

events. Correlation function in the subplot (b) here is calculated for the area of events

shown in Fig. 2a. The same attractor is projected onto the plane {v,u} in subplot (c).

To characterize the difference of the correlated and uncorrelated behavior, an order

parameter has been introduced in ref. [25]. If we define a value hj in block number j by

the condition:

1 u / 0 (the block is moving);

0 u / 0

t
h

  
 

  

, (6)

then the local density of the order parameter H*j can be written as H
*
j=hj(hj+1+hj-1). This

function takes on unit value H*j=hj(hj+1+hj-1)=1 if block j is moving and exactly one of its

nearest neighbors is also moving. Further, H*j=hj(hj+1+hj-1)=2 when both nearest

neighbors of the moving block are in motion. All other cases yield H*j=0. Our

observations with the modified model show that even for blocks at vanishing inter-block

interaction k2→0 (when the motion of neighboring blocks is almost uncorrelated), the

fraction of configurations where sets of neighboring blocks are moving is still relatively

high. All that to say, the combination H*j=hj(hj+1+hj-1) does not vanish in such an

uncorrelated system. It is therefore convenient to construct another simple combination:

1 1

1 all 3 blocks , 1, 1 in contact are moving

0 otherwise
j j j j

j j j
H h h h

 

 
  



. (7)

This makes Hj=hjhj+1hj-1 zero in all cases except when both neighboring blocks of a

central sliding block are also in motion. This combination can be used as an “order

parameter”. To extract integral quantitative information about the steady ordered and

disordered states we calculate also the time evolution of the ensemble-averaged order

parameter: H(t>≡ <Hj(t)>=1/t ∫Hj(t)dt.

As an illustration the dependence of the integrated order parameter on the stiffness k2,

showing a transition from correlated to uncorrelated block motion was calculated for one-

dimensional system. It is shown in Fig.4. Two well defined limiting cases can be seen

easily: the order parameter tends to a constant non-zero asymptote <Hj(t)>→const≠0 at

strong interaction k2>>1 and vanishes <Hj(t)>→0 below the transition point k2≈0.

Block-media Motion Lattice and Movable Automata Numerical Models 7

Fig. 4 Phase transition from correlated to uncorrelated motion of blocks in a one-

dimensional system. The order parameter tends to a constant asymptote at high mutual

interaction k2>>1 and vanishes below transition point k2= k2,critical≈1. Two intermediate

fluctuation regions A and B correspond to the states with short-range and long-range

correlated nonlinear excitations, respectively.

Let us return now to general problem. If the interaction is already varied due to the

motion with nonzero velocity, it sounds natural to modify the elasticity (potential channel

of the interaction) of the system into the state dependent one, which varies at nonzero

velocity and gradually returns to the initial one in static limit:

   
2; ,

2
2; 1, 2; 1, 2; , 1 2; , 1 2; , 1 20 2 ,

4
j nj n j n j n j n j n k j n

k
k k k k k k k v

t
   


       


   . (8)

The same note can be done about the thermalization of the system by a heating produced

by the mechanical work performed by the external force. Local power consumed by the

system in every node of the lattice can be calculated as a work per unit of time:

Q(rj,n)=fj,nvj,n. It can be treated as local source of energy acting in every node and

temperature expansion inside the system should be described by the thermal conductivity

equation:

 , 1, 1, , 1 , 1 , ,4 ( )
j n

j n j n j n j n j n j n

T
T T T T T Q r

t
   


     


 . (9)

It is possible to repeat all the simulations of the previous section, reproduce all the

results presented in Figs.1-3 and demonstrate that these properties are common for all the

modifications of the model. So, let us concentrate below on some additions coming from

generalized version of the model.

Large scale structure of the velocity fluctuations analogous to that plotted in Fig. 1

accompanied by the variations of the elasticity self-consistently calculated according to

Eq. (8) is shown by the colored surfaces (mainly red picks and blue smooth relief,

respectively) in Fig. 5. One can see directly that strong fluctuations of the velocity cause

deep valleys in the surface of elasticity and in turn “prefer” to appear in the same places

again.

8 A.E. FILIPPOV, V.L. POPOV

Fig. 5 Large scale structure of the velocity fluctuations accompanied self-consistently by

the variations of the elasticity are shown by the colored surfaces (mainly red picks and

blue smooth relief, respectively).

Fig. 6 The temperature deviations (bright spots in Matlab color-map ‘hot’ normalized in

the interval [0 2]) caused by the same instant configuration of the system as in Fig. 5.

Green contour lines show median value of the elasticity profile and mark mutual

correlation between its deep valleys and hot regions.

The same correlation exists between the spatially dependent elasticity and the

temperature deviations. In Fig. 6 the temperature deviations (bright spots in Matlab color-

map ‘hot’) caused by the same instant configuration of the system as in Fig. 5 are shown.

Green contour lines show median value of the elasticity profile and mark mutual

correlation between its deep valleys and hot regions.

To extract physically meaningful quantitative information about the earthquakes in the

system different time moments one can calculate a density of power consumed by the

Block-media Motion Lattice and Movable Automata Numerical Models 9

system at each time moment and overlap it other densities of the system like: distributed

rigidity, friction constant and temperature. The local power in the frames of the model

can be calculated as a work produced by the forces acting on each lattice node of the

system pij=fijvij during a unit of time. Calculating spatially distributed power pi=fivi

enables one to evaluate the degree of localization of excitations and hence the fraction of

the regions in which the fragments of the lattice move practically as a whole. As at a

particular site the energy is absorbed and released at different points in time, it is useful to

characterize these sites by the absolute value of the power, |pi|=|fivi|. However, if formally

accumulate the statistics of such “events” using all the nodes as independent realizations,

strongly overestimated value of the small “events” will be obtained. Physically only the

quakes with a magnitude higher than a critical one would be treated as the observable

ones. Moreover, if the amplitude of the power exceeds the critical value for a number of

mutually connected nodes, all of them will be treated by the potential observers as one

common event with an average amplitude in this connected area. So, we organized the

procedure in this manner: all the nodes were found with the amplitudes higher than the

critical one; defined all the connected domains of them and calculated mean amplitudes

for every such a region. To control visually the above procedure, we plotted these

domains by the (bright) ‘jet’ colors over the (blue) map of the rigidity for every time

moment. Typical instant realization of such a picture is presented in Fig. 7. The histogram

calculated for these events, as it is shown in Fig. 8 in double-logarithmic scale.

Fig. 7 Instant connected regions of the “earthquakes” (bright spots) and valleys of the

small rigidity (deep blue) plotted together to illustrate their mutual correlations. The

artificial colors for the magnitude and rigidity are shifted to red and blue parts of the

spectrum to clearly overlap both densities in one picture.

As we already noted above the variations of the friction, elasticity and local velocities

are mutually correlated. In other words, new events are preferably generated in the same

regions (as it is in the reality). As a result, the system partially memorizes such a common

relief for extremely long times. One can calculate, for example, time-averaged relief of

effective rigidity, and we did this. It is quite expected that the variations are found to be

much smaller than local instant variations of the same value, but they exist.

10 A.E. FILIPPOV, V.L. POPOV

Fig. 8 Scaling of the probability to find an earthquake with a particular value of

magnitude calculated for the connected event regions. Physically it corresponds to a

selection of the rare but intensive “events”, which is compatible with the ideology of the

empirical Gutenberg-Richter law.

There are also some correlations between such a memorized structure and preferable

channels of the wave propagation. Again, the correlation is much weaker than it is

observed between instant variations of different densities, but it still exists. These results

are illustrated in Fig. 9, where instant and almost stationary time averaged elastic

coefficient integrated over long time periods ((a) and (b) sub-plots, respectively) shown

by the ‘jet’ color-map. Black contour line depicts mean value of the long-time living

rigidity surface. It marks only partial correlations between the averaged and instant

rigidities, Analogous (only partial) correlation is found between travelling regions of

maximal power production (deep red spots in subplot (b)) and time-averaged relief of the

rigidity.

Fig. 9 Instant and almost stationary time averaged elastic coefficient integrated over long

time periods ((a) and (b) subplots, respectively) shown by the ‘jet’ color-map. Contour

line marks only partial correlations between the averaged and instant rigidities, Analogous

(only partial) correlation is found between travelling regions of maximal power

production (deep red spots in subplot (b)) and time-averaged relief of the rigidity.

Block-media Motion Lattice and Movable Automata Numerical Models 11

3. MOVABLE AUTOMATA MODEL

Generally, the variant of the model which uses movable system is based on a

description of a system of interacting N “particles”. In contrast to described above lattice

model which treats system as deformable elastic medium with varied friction and rigidity,

here each “particle” simulates a complete block (“solid plate”) in some instant position

represented radius-vector ri of its center and the mechanical momentum pi. In simplest

approximation, the plates interact with a potential U(ri-rj) depending on the distance |ri-rj|

between their centers. It formally leads to the following many-body Hamiltonian:

1 , 1

( , ) / 2 (| |) / 2
N N

i i i i i j

i i j

H m U
 

   r p p r r . (10)

In a static configuration the plates are supposed to occupy some equilibrium positions

defined by a competition of the repulsion and attraction forces acting between them. In

the context of this study, it is convenient to represent the interaction potential by a pair of

the Gauss potentials:

2 2(| |) exp{ [( ) / ] } exp{ [( ) / ] }
i j ij i j ij ij i j ij

U C c D d      r r r r r r , (11)

where Cij and Dij define the magnitude, while cij and dij the radii of attraction and

repulsion, respectively. The minimization of the energy corresponding to a desired

equilibrium with some characteristic distance between the centers corresponds to the

inequalities between the parameters: Cij>>Dij, cij<dij. Initially, correct equilibrium

positions of the plates are not known and should be found by a self-assembling transient

process. To obtain it in numerical experiments, the “particles” are initially placed at

random on a rectangle with the side lengths of [0, Lx] and [0, Ly]. During the transient

process the system finds some frozen structure which is a compromise between the

interactions and boundary conditions.

Depending on the particular study, the boundary conditions can be chosen differently.

For example, for the case of shear deformation, it is convenient to employ periodic

boundary conditions on the horizontal axis. It allows simulate potentially infinite run for

relatively small limited numerical arrays. According to the periodic boundary conditions,

a particle leaving the interval [0, Lx] is returned to it at the opposite end and the vectors ri-

rj connecting particles either located within the interval [0, Lx] or the images of “escaped”

particles at the opposite side of the system. On the vertical axis such a system can be

limited by the horizontal plates at y=0 and y=Ly with reflecting boundary conditions:

Uup=C·exp[(y-Ly)/c] and Udown=C·exp[-y/c]. Such a combination of the boundary

conditions simulates together quasi-infinite motion along horizontal axis and confined

motion between two “external” plates in orthogonal direction. The stronger the

inequalities C>>Cij and c<<cij, the more rigid and sharp are the walls in relation to other

forces and lengths relevant to the system.

Keeping in mind typical “tectonic” geometry with the moving continents, in the

present study we are interested in a configuration when along horizontal axis the system is

also limited and moreover, is pressed by a movable “rigid wall”. In analogy with standard

Burridge-Knopoff model, the plate is driven by an external “elastic” force:

12 A.E. FILIPPOV, V.L. POPOV

( )
ext x

F K L V t  , (12)

In this case it is natural, in contrast to the periodic one, to take the system which is limited

along horizontal axis with reflecting boundary conditions: Uright=C·exp[(x-Lx)/c] and

Uleft=C·exp[-x/c], but with the movable right boundary Lx= Lx(t).

Since the parameters Cij, Dij, cij and dij form arrays covering all possible combinations,

it can be said without loss of generality that the Hamiltonian given by Eqs. (10)- (12)

describes systems with any number of components (plates of the different sizes). What is

required is just to specify the parameters Cjk and cjk. In a real system different values of all

the parameters (and what is even more important, the mutual relations between them)

regulate different sizes of the physical plates and their mutual rigidity. Below, to conserve

a generality and get statistically representative results, we use the values Cjk, and cjk

randomly varied with 10 times difference between maximal and minimal ones around unit

Cjk, cjk=1, and reset them from one run to another of every numerical experiment

performed at other fixed or systematically varied parameters.

The simplicity of the equations of motion

/ ( , ) /
i i i i i i

m t H     v x p p f , (13)

is deceptive. In majority of the cases, simulation using such equations is extremely time-

consuming because it involves a summation at every time step over all possible neighbors,

including very distant ones (and even fictitious image particles in the case of periodic

boundary conditions). However, in the case of the solid plates contacting by their

boundaries the number of the mutual interactions can be essentially reduced. The

reduction is based on the following considerations.

It is natural for this particular problem to construct Voronoi diagram of the set of

effective “particles” and associate each cell of this diagram with some of the mutually

contacting plates. According to mathematical definition, a Voronoi diagram is a partition

of a plane into regions close to each of a given set of objects. In the simplest case, these

objects are just set of points in the plane (called seeds). For each seed there is a

corresponding region, called a Voronoi cell, consisting of all points of the plane closer to

that seed than to any other. In our case the central “particles” of every “plate” play a role

of such mathematical “seeds”.

So, the cells of Voronoi diagram give practically the same picture, which we have in

mind when building our dynamic model. Every such a cell is surrounded by the limited

(and in fact quite small) number of the other plates (nearest neighbors). Mathematically

such neighbors are completely defined by the Delaunay triangulation, which is dual to the

Voronoi diagram. According to its definition, Delaunay triangulation for a given set P of

discrete points in a general position is a triangulation DT(P) such that no point in P is

inside the circumcircle of any triangle in DT(P). As result one gets quite realistic

representation of the movable plates and short list of the neighbors interacting with every

particular plate.

Interacting plates exchange momentum, and hence a dissipation channel acting to

equilibrate the relative velocities of particles (that happen to be within the distance cv

close to the equilibrium one) needs to be introduced. This can be done by including an

additive velocity-depending force.

Block-media Motion Lattice and Movable Automata Numerical Models 13

( ) exp [( ) / ]
neigboursN

i i j i j v

f c


     
v

v v r r , (14)

acting on the i-th particle, with some dissipation constant η. As a result, the equations of

motion assume the following form:

i
i i i

m
t


 



r vv
f f . (15)

They can be integrated by using Verlet’s method, which conserves the energy of the

system at each time step, provided no energy is supplied externally through mechanical

work of the movable right boundary which is driven by a balance of the external force

Fext=K(Lx-V0t) and total interaction with the internal plates repulsed by the boundary

potential Uright=C·exp[(x-Lx)/c]

2
1

( ( )) /
N

x
ext right j x j

L
M F U x L t x

t 


    


 . (16)

So, the equations (15) and (16) must be solved self-consistently. Moreover, counting

on the nature of the problem under consideration, this system has to be accompanied by

the equations regulating variation of the effective mutual friction between the plates with

negative derivative in the same spirit as it was done above for the standard BK model on

the elastic lattice:

 1 0 2( , ) ;
j

j j j

F
F F v t v


  


  2 ,<0, 0j nv  and ( ) at 0j jF v v   . (17)

Taking into account an experience which obtained in modification lattice BK model it

seems natural to add here also variable intensity of the interactions between the plates at

nonzero mutual velocities:

 1 1 10 1 2( , ) ;
j

j j j

C
C C v t v


  


  2 ,<0, 0j nv  and 1 10( ) at 0j jC v C v  . (18)

At the same time, here is no need to include directly the temperature conductivity

equation into the model because it is already based on the Brownian dynamic equations

for the plate centres by itself, and in fact automatically involves an exchange by kinetic

energy (temperature) between them. So, plotting the kinetic energy of the “particles” by

the artificial colors one gets a color-map in the representation which is natural for the

model based on the movable automata. The results obtained in the frames of this model

are summarized in the Figs. 10-15 below.

Typical instant configuration of the plates pressed from one of the boundaries by an

external force is shown in Fig. 10. Let us remind that the rigid wall presses the system

from the right side. As a result, it becomes more compressed form this side and areas of

the Voronoi cells here become smaller than in average inside the system far from this

boundary. The mass of every plate conserves in the frame of the model. The only are of it

changes. Gray-scale map of the subplot (a) is used to show the individual areas of the

plates.

14 A.E. FILIPPOV, V.L. POPOV

The compression and motion of the plates causes nonzero mutual velocities of every

plate in relation to their proximities. This fact is well seen in the subplot (b) by the

MatLab ‘jet’ color-map. According to Eq.(17) nonzero mutual velocities result in a

variation of the effective friction between every plate and its surrounding. This

characteristic of the system is depicted in the subplot (c). It is important to note also that

the Delaunay triangulation is recalculated at every time step. As a result, the links

between the plates moving with different velocities can appear and disappear depending

on their physical positions at a particular time moment. However, in sufficiently dense

system (which is used here) such changes of the neighbors happen quite rarely. So, the

model seems quite realistic to describe some features of the tectonic motion.

Fig. 10 Typical instant configuration of the plates pressed from one of the boundaries by

an external force. Gray-scale map of the subplot (a) shows the areas of the plates more

compressed from one side (white color) than inside the system (different intensity of gray

color). The compression and motion of the plates causes different (mean) mutual velocity

of every plate in relation to its proximity which results in the variation of the effective

friction between this plate and surrounding ones. Both these characteristics are plotted by

the ‘jet’ colors.

Block-media Motion Lattice and Movable Automata Numerical Models 15

Fine structure of the model is presented in Fig. 11. Movable Voronoi diagram in

subplot (a) here show the plates in the frame of the model, as above. However, here they

are accompanied by the colored seeds of the Voronoi cells. The colors show different

instant velocities of the plates. Dual to the above diagram Delaunay triangulation shows

connections between the plates and their neighbors counted in the model. According to its

mathematical definition it naturally gives a list of the neighbors for every plate. So, one

can calculate the interactions between them, solve the dynamic equations, so on. Different

sizes of the nodes here reproduce different (randomly distributed) interaction parameters.

Besides, it is convenient to control instant number of them for every cell (symmetry of the

plate polygon) and visualize it by the colors from the same standard ‘jet’ map.

Fig. 11 Movable Voronoi diagram in subplot (a) describes the plates in the frame of the

model. The colors of the seeds show different instant velocities of the plates. Dual to the

above diagram Delaunay triangulation shows connections between the plates and their

neighbors counted in the model. Different sizes of the nodes here reproduce different

interaction coefficients and their colors show different instant number of neighbors.

16 A.E. FILIPPOV, V.L. POPOV

Two projections of the phase portrait on different planes (subplots (a) and (b)) are

presented in the Fig. 13 by the colored circles of different sizes. The colors and sizes of

the circles have the same meaning as in the Fig. 12a. Green line with the arrows shows

typical individual trajectory, which every point passes through the phase space with the

time and reproduces the evolution of the effective friction due to periods of large or small

local velocities. Generally, the subplot (b) here has the same physical meaning as the

phase portrait plotted in Fig. 3a in the lattice variant of BK model.

Fig. 12 The same instant configuration as in Fig. 11 reproducing mutual relations between

the interactions, local rigidities, velocities and friction coefficient directly marked in the

plots. The colors and sizes of the circles in the subplot (a) depict velocities and interaction

coefficients of the individual plates, respectively. Yellow, cyan, magenta and red lines in

the subplot (a) display integrated over y-coordinate averaged x-dependencies of the

velocity, rigidity, effective friction coefficient density and temperature respectively. The

colors of the circles in the subplot (b) correspond to different individual power of the

events in every plate.

Block-media Motion Lattice and Movable Automata Numerical Models 17

To get the local power in the frames of the discrete variant of the model we calculate

the work produced by the forces acting on each plate (seed) of the system pi=fivi.

Calculation of spatially distributed power pi=fivi enables one to evaluate the degree of

localization of excitations and hence the fraction of the regions in which the fragments of

the lattice move practically as a whole. As at a particular site the energy is absorbed and

released at different points in time, it is useful to characterize these sites by the absolute

value of the power, |pi|=|fivi|. It is shown in subplot (b) of the Fig.12.

Summation over all particles yields the total power VFext(t)=P(t)=∑fivi. It gives an

energy which is consumed during every time unit due to motion of rigid right boundary,

which is forced by external flow with constant velocity, V=const. However, the boundary

here is moving with some varied velocity V(t) which is determined by a balance of the

total force, Fext(t) and only its time average coincides with <V(t)>=V0 .

Fig. 13 Two projections of the phase portrait on different planes (subplots (a) and (b)).

The colors and sizes of the circles depict velocities and interaction coefficients of the

individual plates, respectively. Green line with the arrows shows typical individual

trajectory, which every point passes in phase space with the time.

18 A.E. FILIPPOV, V.L. POPOV

Fig. 14 Instant configuration of the density (surface in subplot (a)) and time-space

diagram of the history of events (colormap in subplot (b)). Bold black contours mark

intensive “earthquake” events near the propagating right boundary. Different speed of the

front and some “epochs” of intensive events can be seen from different inclination of the

front and horizontal lines in the recorded history.

In particular, it causes quasi-periodic oscillations of the velocity of boundary which

are typical for the “stick-slip” nature of the processes. It is actually the same phenomenon,

which takes place in ordinary BK model, shown above in Figs. 2-3. In the discrete model

one can additionally observe the waves of the velocities, densities and power of the

events, which are accompanied by quasi-periodic “epochs” of intensive “earthquakes”.

Instant configuration of the density and history of the density integrated along y-axis

corresponding to the magenta curve in Fig. 12a are shown in Fig. 14. The instant surface

is plotted in subplot (a) and corresponding time-space diagram of the history of events is

presented in the subplot (b). Bold black contours mark intensive “earthquake” events near

the propagating right boundary. Different speed of the front and some “epochs” of

intensive events can be seen from different inclination of the front and deeply colored

almost horizontal lines in the recorded history.

Exactly the same history of the events (recorded during the same calculation run) as in

Fig. 14, is depicted by the intensity of power production (shown by the ‘hot’ color map).

The epochs of the intensive earthquakes are much stronger pronounced here than in the

previous figure.

Block-media Motion Lattice and Movable Automata Numerical Models 19

Fig. 15 Exactly the same history of the events, as in Fig.14, depicted by the intensity of

power production (shown by the ‘hot’ colormap). The epochs of the intensive earthquakes

are much stronger pronounced here than in the previous figure.

4. CONCLUSION

It is shown that original Burridge-Knopoff model allows a number of the

modifications which can increase its ability to be adapted to description of the different

practically important features of the tectonic phenomena. Basic modifications of the

model were already made previously. They incorporate naturally introduced state

dependent friction where the friction force is defined by the additional kinetic equation as

well as extension of the model to more realistic two-dimensional case.

In the present paper we continued previous modifications and mainly concentrated on

the following extensions of the model. First of all, elastic force and thermal conductivity

were added into consideration. Besides to the velocity weakening of friction between

moving blocks, we combined an increase of both static friction and rigidity during stick

periods as well their simultaneous weakening during motion. It is expected to give better

and qualitatively correct description of plural transitions from smooth sliding to stick-slip

behavior, which are normally observed in various tribological and tectonic systems. To

20 A.E. FILIPPOV, V.L. POPOV

get easier comparison with the previous versions of the model we initially did this using in

standard approach based on interacting elastically connected mesh elements.

An alternative version of the model is finally proposed. It applies a simulation of the

motion of relatively large solid blocks, utilizing technique of the movable cellular

automata. In this advanced version we also involve all the components of the real system:

state depending friction, changeable rigidity, as well as heat production and thermal

conductivity due to dynamics of the interacting plates, which were already included into

the previous versions based on traditional lattice model.

Proposed modification has its own advantages and disadvantages, which were

discussed in the paper step by step. However, a combination of the different approaches

gives additional degrees of freedom, which allows select more convenient or adequate

approach depending on the particular problem under consideration. We plan to continue

further and combine our studies in both alternative directions.

REFERENCES

1. Gutenberg, B., Richter, C.F., 1994, Frequency of earthquakes in California, Bulletin of the

Seismological Society of America, 34(4), pp. 185-188.

2. Bak, P., Christensen, K., Danon, L., Scanlon, T., 2002, Unified scaling law for earthquakes, Physical

Review Letters, 88, 178501.

3. Žalohar, J., Omori's Law, 2018, Chapter 10 in the book “Developments in Structural Geology and

Tectonics”, Elsevier, pp. 123-134.

4. Bak, P., Tang, C., Wiesenfald, K., 1987, Self-organized criticality: An explanation of the 1/f noise,

Physical Review Letters, 59, 381.

5. Carlson, J.M., Langer, J.S., 1989, Properties of earthquakes generated by fault dynamics, Physical

Review Letters, 62, 2632.

6. Ruzhich, V.V., Smekalin, O.P., Shilko, E.V., Psakhie, S.G., 2002, About nature of “slow waves” and

initia-tion of displacements at fault regions, In: Proceedings of the international conference on new

challenges in mesomechanics, pp.311-318.

7. Ruzhich, V.V., Truskov, V.A., Chernykh, E.N., Smekalin, O.P., 1999, Neotectonic movements in fault

zones of the Baikal region and their origin mechanisms, Russian Geology and Geophysics, 40(3), pp.

356-368.

8. Psakhie, S.G., Shilko, E.V., Astafurov, S.V., 2004, Peculiarities of the mechanical response of

heterogeneous materials with highly deformable interfaces, Technical Physics Letters, 30, pp. 237-239.

9. Psakhie, S.G., Ruzhich, V.V., Shilko, E.V., Popov, V.L., Dimaki, A.V., Astafurov, S.V., Lopatin, V.V.,

2005, Influence of the state of the inter-block interfaces on the character of local displacements in

block-structured media, Technical Physics Letters, 31, pp. 712-715.

10. Ruzhich, V.V., Psakhie, S.G., Bornyakov, S.A., Smekalin, O.P., Shilko, E.V., Chernykh, E.N.,

Chechelnitsky, V.V., Astafurov, S.V., 2002, Investigation of influence of vibroimpulse excitations on

regime of displacements in seismically active fault regions, Phys. Mesomech., 5(5), pp. 85-96.

11. Bhattacharya, P., Chakrabarti, B.K., Samanta, K.D., 2009, Fractal models of earthquake dynamics,

Heinz Georg Schuster (ed), Reviews of Nonlinear Dynamics and Complexity, Wiley-VCH, pp. 107–150.

12. Ford, J., 1992, The Fermi-Pasta-Ulam problem: Paradox turns discovery, Physics Reports, 213(5), pp.

271-310.

13. Lichtenberg, A.J., Lieberman, M.A., 1992, Regular and Chaotic Dynamics, Springer, New York, 692 p.

14. Fillipov, A.E., Hu, B., Li, B., Zeltser, A., 1998, Energy transport between two attractors connected by a

Fermi-Pasta-Ulam chains, Journal of Physics A: Mathematical and General, 31, 7719.

15. Filippov, A.E., Popov, V.L., Psakhie, S.G., 2006, Converting displacement dynamics into creep in block

media, Technical Physics Letters, 32(6), pp. 545-549.

16. Filippov, A.E., Popov, V.L., Psakhie, S.G., 2008, Correlated impacts optimizing the transformation of

block me-dium dynamics into creep regime, Technical Physics Letters, 34(8), pp. 689-692.

Block-media Motion Lattice and Movable Automata Numerical Models 21

17. Filippov, A.E., Popov V.L. 2021, Study of Dynamics of Block-Media in the Framework of Minimalistic

Numerical Models, Chapter 7 in the book «Multiscale Biomechanics and Tribology of Inorganic and

Organic Systems», Springer International Publishing, pp. 143-168.

18. Burridge. R., Knopoff, L., 1967, Model and Theoretical Seismicity, Bulletin of the Seismological

Society of America, 57(3), pp. 341-371.

19. Bruce, S.E., Carlson, J.M., Langer, J.S., 1992, Patterns of seismic activity preceding large earthquakes,

Journal of Geophysical Research, 97, pp. 479-488.

20. Carlson, J.M., Langer, J.S., Shaw, B.E., 1994, Dynamics of earthquake faults, Review of Modern

Physics, 66(2), pp. 657-670.

21. Saito, T., Matsukawa, H., 2007, Size dependence of the Burridge-Knopoff model, Journal of Physics:

Conference Series, 89, 012016.

22. Mandelbrot, B., 1982, The Fractal Geometry of Nature, Freeman and Co., San Francisco, 460 p.

23. Pelletier, J.D., 2000, Spring-block models of seismicity: review and analysis of a structurally

heterogeneous model coupled to the viscous asthenosphere. Geocomplexity and the Physics of

Earthquakes, Geophysical Monographs, American Geophysical Union, pp.27-42.

24. Muratov, C.B., 1999, Trave ling wave solutions in the Burridge-Knopoff model, Physical Review E, 59,

3847.

25. Huisman, B.A.H., Fasolino, A., 2005, Transition to strictly solitary motion in the Burridge-Knopoff

model of multicontact friction, Physical Review E, 72, 016107.

26. Clancy, I., Corcoran, D., 2005, Criticality in the Burridge-Knopoff model, Physical Review E, 71,

046124.

27. Clancy, I., Corcoran, D., 2006, Burridge-Knopoff model: Exploration of dynamic phases, Physical

Review E, 73, 046115.

28. Mori, T., Kawamura, H., 2008, Simulation study of earthquakes based on the two-dimensional

Burridge-Knopoff model with long-range interactions, Physical Review E 77, 051123.

29. Dieterich, J.H., 1974, Earthquake mechanisms and modeling, Annual Review of Earth and Planetary

Sciences, 2, pp. 275-301.

30. Dieterich, J.H., 1978, Time-dependent friction and the mechanics of stick-slip, Pure and applied

geophysics, 116, pp. 790-806.

31. Rice, J.R., Gu, J.C., 1983, Earthquake aftereffects and triggered seismic phenomena, Pure and applied

geophysics, 121, pp. 187-219.

32. Gu, J.C., Rice, J.R., Ruina, A.L., Tse, S.T., 1984, Slip motion and stability of a single degree of freedom

elastic system with rate and state dependent friction, Journal of the Mechanics and Physics of Solids,

32(3), pp. 167-196.

33. Marone, C., 1998, Laboratory-derived friction laws and their application to seismic faulting, Annual

Review of Earth and Planetary Sciences, 26, pp. 643-696.

34. Persson, B.N.J., 2000, Sliding Friction, Physical Properties and Applications, Springer, Berlin.

35. Popov, V.L., 2010, Contact Mechanics and Friction Physics, Springer, Berlin, 362p.

36. Popov, V.L., 2000, A theory of the transition from static to kinetic friction in boundary lubrication

layers, Solid State Communications, 115(7), pp. 369-373.

37. Filippov, A.E., Klafter, J., Urbakh, M., 2004, Friction through dynamical formation and rupture of

molecular bonds, Physical Review Letters, 92, 135503.

38. Fermi, E., Pasta, J., Ulam, S., Tsingou, M., 1993, The Many-body Problem: An Encyclopedia of Exactly

Solved Models in One Dimension, World Scientific, Singapore, 988 p.

39. Braun, O.M., Hu, B., Filippov, A.E., Zeltser, A., 1998, Traffic jams and hysteresis in driven one-

dimensional systems, Physical Review E, 58, 1311.

40. Filippov, A.E., 1994, Mimicry of phase transitions and the large-river effect, JETP Lett., 60(2), pp. 133-

137.

41. Filippov, A.E., Popov, V.L., 2010, Modified Burridge–Knopoff model with state dependent friction,

Tribology International, 43, pp. 1392-1399.

42. Ohmura, A., Kawamura, H., 2007, Rate-and state-dependent friction law and statistical properties of

earthquakes, Europhysics Letters, 77, 69001.

43. Persson, B.N.J., Popov, V.L., 2000, On the origin of the transition from slip to stick, Solid State

Communications, 114(5), pp. 261-266.

44. Balokhonov, R., Romanova , V., Zemlianov, A., 2021, A mesoscopic analysis of a localized shear band

propagation effect on the deformation and fracture of coated materials, Reports in Mechanical

Engineering, 2(1), pp. 6-22.

22 A.E. FILIPPOV, V.L. POPOV

45. Romanova, V., Balokhonov, R., Zinovieva, O., 2021, Mesoscale deformation-induced surface

phenomena in loaded polycrystals, Facta Universitatis-Series Mechanical Engineering, 19(2), pp. 187-

198.

46. Kocharyan, G.G., Kishkina, S.B., 2021, The Physical Mesomechanics of the Earthquake Source, Phys.

Mesomech., 24, pp. 343–356.

47. Makarov, P.V., Khon, Y.A. Autosoliton View of the Seismic Process. 2021, Part 1. Possibility of

Generation and Propagation of Slow Deformation Autosoliton Disturbances in Geomedia, Phys.

Mesomech., 24, pp. 363–374.

48. Makarov, P.V., Smolin, I.Y., Khon, Y.A., et al. 2021, Autosoliton View of the Seismic Process. Part 2.

Possibility of Generation and Propagation of Slow Deformation Autoso liton Disturbances in

Geomedia, Phys. Mesomech., 24, pp. 375–390.