Discrete rearranging disordered patterns, part I: Robust statistical tools in two or three dimensions
Discrete rearranging patterns include cellular patterns, for instance liquid foams, biological tissues, grains in polycrystals; assemblies of particles such as beads, granular materials, colloids, molecules, atoms; and interconnected networks. Such a pattern can be described as a list of links between neighbouring sites. Performing statistics on the links between neighbouring sites yields average quantities (hereafter "tools") as the result of direct measurements on images. These descriptive tools are flexible and suitable for various problems where quantitative measurements are required, whether in two or in three dimensions. Here, we present a coherent set of robust tools, in three steps. First, we revisit the definitions of three existing tools based on the texture matrix. Second, thanks to their more general definition, we embed these three tools in a self-consistent formalism, which includes three additional ones. Third, we show that the six tools together provide a direct correspondence between a small scale, where they quantify the discrete pattern’s local distortion and rearrangements, and a large scale, where they help describe a material as a continuous medium. This enables to formulate elastic, plastic, fluid behaviours in a common, self-consistent modelling using continuous mechanics. Experiments, simulations and models can be expressed in the same language and directly compared. As an example, a companion paper Cartes2D provides an application to foam plasticity.
pacs:62.20.DElasticity and 62.20.FDeformation and plasticity and 83.50.-vDeformation and flow in rheology
Cellular patterns include: liquid foams or emulsions (Fig. 1); crystalline grains in polycrystals; or biological tissues (Fig. 2). Assemblies of particles (Fig. 3) include collections of beads, molecules, or atoms; granular or colloidal materials; sets of tracers dispersed in a material, such as fluorescent probes or passively carried particles. Despite their tremendous diversity of sizes and physical properties, all these patterns have a common point: they are made of a large number of well-identified individual objects. We call them discrete patterns, where the word "discrete" here means the opposite of "continuous". Other discrete patterns include interconnected networks, of e.g. springs, polymers, biological macromolecules, fibers, or telecommunication lines.
We define the pattern as rearranging if the mutual arrangement of the individual objects can change. This is the case if they can move past each other, for instance due to mechanical strain (Figs. 1 or 2a), spontaneous motility (Fig. 2b), or thermal fluctuations (Fig. 3). This is also the case if the number of individual objects can change, for instance due to cell division or death (Fig. 2b), coalescence or nucleation of bubbles, shrinkage during coarsening of polycrystals or foams.
|eq. (11)||eq. (10)|
|eq. (14)||eq. (20)||eqs. (17,18)|
Stimulated by the various imaging techniques, ref. mecke reviews many tools available to describe and quantitatively characterise a pattern (that is, a single image). Another tool is the texture: it appears in various contexts, including the order parameters of nematics, the microstructure of polymers, or the fabric of grains; and has been used to describe mechanical strains by Aubouy et al. (see aub03 and references therein). It describes statistically how the individual objects are arranged with respect to each other. With a few simple measurements performed directly from an image, it extracts quantitative information relevant to the size and anisotropy of the pattern. Since it is based on statistics, it is particularly useful for disordered patterns. Here, we present a coherent set of robust tools (listed in Table 1), with a triple goal.
First, we revisit existing definitions of the texture , as well as the statistical strain aub03 and the rearrangements dollet_local based on it. In fact, cellular patterns are better characterised using cell centers than using their vertices. This remark enabled ref. dollet_local to define a preliminary version of . Here, we show that it also yields a more general definition of the texture, valid for all patterns. In addition, since a cell center is measured as an average over all pixels in a cell (while a vertex is a single pixel), measurements in experiments or simulations are more robust. The companion paper Cartes2D shows that displacements of cell centers (but not of cell vertices) are close to affine displacements. In order to make this paper self-contained, we recall and hopefully clarify the definitions of and . We also present a more general definition of and derive explicitly its prefactor.
Second, thanks to their more general definition, we embed these three tools in a self-consistent formalism, which includes three additional ones: , and . From two successive images in a movie, we extract information regarding the magnitude and direction of strain rate and rearrangements.
Third, we show that the six tools together provide a direct correspondence between a small scale, where , and quantify the discrete pattern’s local distortion and rearrangements; and a large scale, where , and help describe a material as a continuous medium without any details related with the discrete scale. This enables to formulate elastic, plastic, fluid behaviours in a common, self-consistent modelling using continuous mechanics even for a discrete material. Experiments, simulations and theories can be expressed in the same language to be directly compared.
The only requirement is that the image should be of sufficient quality to extract the positions of the centers of each individual object (cell or particle); as well as the list of neighbour pairs (which objects are neighbours). All tools here are either static or kinematic, and rely on the image only; that is, they are independent of dynamics (stresses, masses and forces). They apply to discrete patterns regardless of the size of their individual objects, which can range from nanometers to meters or more. They regard simulations as well as experiments, and should enable quantitative comparison between them. They apply whatever the pattern’s disorder is.
Our equations are valid in any dimensions. For clarity, we write them in 3D, and show that is is straightforward to rewrite them in 2D, see section 2.2.3. We specifically choose to illustrate this paper with 2D images (Figs. 1-3), which are simpler and more common that 3D data.
More precisely, we illustrate each definition on the example of a foam flow (Fig. 1a) rau07 , which is both our original motivation and the most suitable example. Nitrogen is blown into water with commercial dishwashing liquid. Bubbles enter a channel, of length 1 m (only partly visible on the picture), width 10 cm, and thickness 3.5 mm: a monolayer of bubbles forms (area mm), sandwiched between two glass plates (quasi-2D foam, liquid fraction less than a percent). It steadily flows from left to right without vertical component (true 2D flow) until it reaches the free end of the channel. Coalescence and ageing are below detection level. A 3 cm diameter obstacle is inserted into the foam channel. The foam is forced to flow around it, resulting in a spatially heterogeneous velocity field. Different regions simultaneously display different velocity gradients, internal strains, and rearrangement rates, and allow to sample simultaneously many different conditions. Bubbles naturally act as tracers of all relevant quantities; and on the other hand the foam’s overall behaviour appears continuous. The total strain rate is partly used to deform bubbles and partly to make them move past each other; the companion paper Cartes2D studies how it is shared between both contributions.
Sections 2.1 and 2.2 start from the static description of Ref. aub03 and develop it step by step, for pedagogical purpose, while refining it. Section 2.3 describes the changes between two successive images. Section 3 is useful to compare measurements on different patterns; or to compare experiments and simulations. Section 4 is more theoretical and regards specific applications: it discusses how to characterize materials which behave as continuous media, that is, where the quantities vary smoothly with space; and when it is possible to identify our statistical tools with the usual quantities of continuous mechanics. Appendices cover many practical, technical or theoretical details including all standard definitions and notations of matrices used in this paper.
2 Texture and time evolution of links in the discrete pattern
The pattern is a collection of individual objects. Here we are interested in the relative positions of these objects, not in each object’s shape (although both are related in particular cases such as cellular patterns, see Appendix A.3.3). We thus replace each object by a point called "site".
The user should adapt the measurement tools to the pattern under consideration, and the scientific questions to be answered. For that purpose, the user should begin by deciding: (i) what are the relevant links, that is, pairs of sites which are connected (section 2.1.1); and (ii) the averaging procedure (section 2.1.2 and Appendix A.1).
These choices are conventions, and thus rather free. The results of the measurements depends on the chosen definition, but they are much more robust than scalar measurements (see section 2.3.3). Moreover, as long as the same definition is used for all measurements, the equations that relate the measurements of the different quantities (such as eq. 9) are valid independently of the chosen definition. Once conventions are chosen, it is thus important to use them consistently.
2.1.1 Links between neighbouring sites
In a cellular pattern (Figs. 1,2), it is often advisable to choose as sites each cell’s geometrical center (see Fig. 4a). However, alternative choices exist. For instance, a user interested in studies of dynamics might prefer the center of mass, if different from the geometrical center. Similarly, a biologist might be more interested in the cell’s centrosome or nucleus. Note that we do not advise to use a definition based on vertices (see Appendix A.3.2).
When two cells touch each other it defines that their sites are connected. This is unambiguous if cells walls are thin. This is the case for grains in polycristals, cells in an epithelium, or in a foam with low amount of water (Figs. 1a,2).
If cells walls are thick, as is the case in a foam with a higher amount of water (Fig. 1b), different definitions of neighbours are possible. For instance, two cells are defined as neighbours if their distance is smaller than a given cut-off. Or, if they are neighbours on a skeletonized image; that is, after an image analysis software has reduced cell walls to one pixel thick black lines on a white background. If cell walls are too thick, cells are really separated (as in a bubbly liquid, where bubbles are round and far from each other) and can be treated like the particles, which we now discuss.
If each object is a particle (as in Fig. 3) it is natural to choose its center as site (see Fig. 4b). There are various possible choices for the links. Since the tools characterize patterns and not forces, the definition of links is independent of interactions: a link between sites does not mean that sites interact; conversely, sites which interact are not necessarily linked. Whatever the chosen definition, it is important that each particle has only a finite number of neighbours.
In a first case (Fig. 3b), the average distance between particles is comparable to their average radius; for instance, for a dense (also called compact or jammed) colloid or granular material. We then recommend to define that two particles are linked if their distance is less than a chosen cut-off. For hard spheres, this cut-off should be the sphere’s radius plus a small tolerance.
In the opposite case, the average distance between particles is much larger than their average radius (Fig. 3a); for instance, for a decompacted colloid or granular material. We then recommend to recreate a cellular pattern by attributing to each particle its Voronoi domain (the set of points surrounding this particle, closer to it than to any other particle). One then chooses to define that two particles are linked if their Voronoi domains touch; this is called the "Delaunay triangulation" of the particles.
If the pattern is a network, it is natural to choose the nodes as sites. The connexions are physically materialised, and thus unambiguously defined.
The present tools aim at describing the collective properties of links. In what follows, denotes the average over a set of links relevant to the user:
where the sum is taken over the number of such links. Appendix A.1 presents some technical details, especially regarding the boundaries of the averaging region, which can be treated as sharp or smooth.
The scale of study determines the number of links included. Performing the same analysis at different scales (Fig. 6) enables to obtain multi-scale results durand2004 ; dollet_local ; jan05 . For instance, we can measure the dependence of pattern fluctuations with scale for particle assemblies.
Choosing to average over a small number of links yields access to detailed local information. For instance, the local heterogeneity of a sample of ice can be measured by including the links around one single grain, then performing a comparison between different grains durand2004 . Similarly, to study the anisotropy of a biological cell which divides, one can include only the links starting at this cell’s center courtypreprint .
On the other hand, choosing a large number of links yields better statistics. This is the case for instance if the system is homogenous in space. In a homogenously sheared foam (Fig. 1b), it makes sense to consider that all bubbles play a similar role, and average over the whole foam. Averaging over all links contained in the whole image enables to detect the overall anisotropy of an ice sample or an epithelium, and compare it with other samples durand2004 ; courtypreprint .
Even if the system is invariant along only one direction of space, one can average over this direction jan05 . Similarly, in a flow which is invariant in time, one can average over time (Fig. 1a). In what follows, figures are prepared with successive video images, representing millions of bubbles: a time average yields good statistics and details of local variations, even if only a small part (i.e. few links) of each image is included (Fig. 6c).
2.2 Texture : current state of the pattern.
2.2.1 Definition and measurement
A pair of neighbour sites of coordinates and constitutes a link. Reducing the pattern to a set of links sets aside the detailed information regarding the actual positions of each site.
The link vector:
has coordinates . It carries the information on link length and angle. However, and play the same physical role: an average over several s will yield a result which depends on this arbitrary choice of sign (and, in practice, if there are enough links, the average turns out to be close to zero).
On the other hand, the number is invariant under the change and thus has a physically relevant (and non-zero) average: . It reflects the average square link length, but loses the information of angle.
The link matrix combines the advantages of both:
Its trace is . Its average defines the texture aub03 :
It is expressed in m. As required, it stores the same information regarding the current pattern: the square length, readily visible as the sum of diagonal terms; the angle and magnitude of anisotropy, as discussed below.
2.2.2 Diagonalisation and representation
By construction, is a matrix with symmetric off-diagonal terms (, etc…). It can thus be diagonalised (see Appendix B.2 for details):
Its three eigenvalues (, 2 or 3) are positive. Their sum, Tr, is exactly . In practice, they usually have the same order of magnitude, of order of .
In a truly 3D pattern, has strictly positive eigenvalues (except in unphysical examples). Thus its inverse always exists (eq. 55).
can be represented as an ellipsoid, which axes directions are that in which is diagonal, represented as thin lines on Fig. (5). Each ellipsoid’s axes length is proportional to the corresponding . It is expressed in m: this is less intuitive than m, and ellipsoids are more elongated that the actual cell shape; but this is necessary for the consistency with the representation of the other matrices (see Appendix A.3.3). The square link length is reflected in the size of the ellipsoid, more precisely as the square root of the sum of the three axes lengths; it is thus not proportional to the ellipsoid’s volume. The direction in which links are longer is represented by the direction of ellipsoid elongation: the greater the pattern’s anisotropy, the more elongated the ellipsoid. If the texture is measured at several regions of the image, it is represented as several ellipsoids, that is, a map of the texture field (see also section 2.2.3 and Fig. 6).
When the pattern is statistically isotropic, so is its texture. It is thus diagonal in any system of axes, and the three s are exactly equal, :
That is, the texture of an isotropic pattern contains only the information of length: , where is the identity matrix in 3D. It is represented as a sphere. In that case, all axes are equivalent (or "degenerated").
2.2.3 Two-dimensional case
If the pattern under consideration is contained in a plane, as are most experimental images, we turn to a 2D notation. As mentioned, this is straightforward:
There exist two orthogonal axes (eigenvectors) in which would be diagonal:
with strictly positive , ( or 2). If we call the r.m.s. length of links in the direction of elongation (say, 1) and the r.m.s. length of links in the direction of compression (say, 2), then and . Its inverse always exists.
In 2D, is represented by an ellipse (Fig. 5). Measurements can be performed at larger scale to decrease the noise due to fluctuations, or at smaller scale to evidence more details of the spatial variations (Fig. 6).
When the pattern is isotropic, its texture is diagonal with any choice of axes:
That is, , where is the identity matrix in 2D. All axes are equivalent (or "degenerated") and the angle of eigenvectors is not defined. is represented by a circle and the thin lines lose their signification (Fig. 5A).
2.3 Time evolution
Differentiating eq. (3) determines how varies. Appendix C.1, useful for practical calculations, discusses finite size effects due to the time interval between successive images of a movie. Neglecting these effects in eq. (66) yields the simplified time evolution of :
The variation of is negligible in most physical examples; however, it is significant for instance in biological tissues with many divisions courtypreprint or in coarsening systems such as ageing foams wea99 ; weairerivier .
The three terms of the r.h.s. can be measured on a movie, and have the following meaning. In the time interval between two successive images, some links enter or exit the region of averaging; some links change their length or angle; some links are created or destroyed, respectively. They are now discussed one by one (sections 2.3.1, 2.3.2 and 2.3.3; respectively).
In eq. (9), is the flux of advection, that is, the transport of texture. It counts the rate at which links enter or exit throught the sides of the region of averaging. Technically, it is a rank-three tensor (i.e. with 3 indices): for more details see ref. dol05 and Appendix B.3. In a good approximation, , where is the local average velocity, see Appendix C.1. Its divergence counts the net balance between links that enter and exit; it vanishes if is spatially homogeneous, or at least is constant along a flux line; it also vanishes if the local average velocity is zero.
2.3.2 Geometrical texture changes:
describes the changes in the pattern’s overall shape, that is, geometry: at which rate, and in which direction, the pattern deforms. It reflects relative movements: it is insensitive to a global, collective translation.
is based on the links which, on both successive images, exist (i.e. do not undergo topological rerrangement) and belong to the region of averaging (i.e. are not advected). These links may change in length and direction (Fig. 7). We obtain each link’s contribution directly from eq. (2) and average it, like in eq. (3):
Appendix C.1 provides more details; it also defines , which symmetrical part is related with , and will turn useful to define and in section 3.2.1. The right hand side of eq. (10) should not be confused with ; that is, is not equal to the variation of : as eq. (9) shows, the difference between them is due to changes in the links included in the averages.
is symmetric. Its units are in ms. It has a positive eigenvalue in a direction of extension, and a negative one in a direction of compression. In case of dilation, all its eigenvalues are positive. In a region of shear, it is plotted as an ellipse with one thin line drawn on it, similar to a "coffee bean" (Fig. 7). We can plot a map of : it is similar to Fig. (10), data not shown.
2.3.3 Topological texture changes:
reflects the topological changes, namely, changes in the list of links: creation and destruction, that is, source term of the texture. Appendix C.1 provides details. Briefly, each link which has appeared since the preceding image has a contribution given by eq. (2), noted ; and similarly the contribution of a link which disappears before the next image is noted . Averaging over all links which appear or disappear between successive images defines as:
The quantity (resp. ), expressed in s, is not the time derivative of a physical quantity (which we would note ). It is the rate of link appearance (resp. disappearance), per unit time and per existing link. If and are equal, their inverse is the average link’s life expectancy.
is expressed in ms. It characterises the total effect on the pattern of all topological changes occurring between two images. By construction it is symmetric, like : it can thus be diagonalized and represented as an ellipsoid. It is robust to artefacts and errors in determinations of neighbours dollet_local . It is general, and includes information of frequency, size, direction and anisotropy for all contributions of all topological changes: they can be treated indifferently and added together. However, as we now discuss, the user might be interested in studying separately the contributions of the different processes courtypreprint .
The coalescence of two sites (Fig. 8d) corresponds in foams to the breakage of a liquid wall between two bubbles, with a net balance of minus one site ohl98 . The reverse process corresponds in epithelia to a cell division, and results in one more site dub98 . When the number of sites decreases (resp. increases), so does the number of links, and usually has only negative (resp. positive) eigenvalues. The variation in the number of sites and links is thus visible in the trace of .
Neighbour exchanges (Fig. 8a,b), also called "T1" in the case of cellular patterns wea99 ; weairerivier , preserve the sites; in 2D, they also preserve the number of links, . It is for that case that ref. dollet_local introduced a specific definition of . usually is mostly deviatoric (Appendix B.1), with both positive and negative eigenvalues (and a vanishing trace), in directions correlated with the appearing and disappearing links. Note that the eigenvectors are exactly orthogonal, while the appearing and disappearing links need not be: thus eigenvectors are not strictly parallel to links, especially when is measured as an average over several individual topological processes. For our flowing foam example, a map of (data not shown) would be similar to Fig. 18 of ref. dollet_local , or to Fig. 11 below.
If a pattern has a free surface, when sites exchange neighbours the number of links can vary, and thus locally can be far from deviatoric.
The disparition of a site (Fig. 8c) corresponds in foams to a bubble which shrinks, also called "T2" wea99 ; weairerivier ; and in epithelia to a cell which dies, or exits the epithelium plane. The reverse process is a site nucleation. Both processes have an approximately isotropic contribution to .
3 Statistical tools to obtain relative deformations and their time evolution
This section facilitates comparison between different experiments; or between experiments, simulations and theory. This is useful for large scale deformations and flows: e. g. of particle assemblies, of foams and emulsions, of granular materials, or of biological tissues during morphogenesis.
Here we try to link the discrete, local description of Section 2 with the continuous, global description of Section 4. For this continuous description to be self-consistent, it is necessary to get rid of the discrete objects’ length scale, that is, the typical size of links. For each of the three discrete quantities , , defined in section 2, it is possible to construct a continuous counterpart: that is, a tool which is dimensionless (or expressed in s), with no m any longer. This defines , , , respectively (sections 3.1, 3.2.1 and 3.2.2, respectively).
3.1 Statistical internal strain :
The internal strain has been defined by Aubouy et al. aub03 through a comparison between the current pattern and a reference one. We include it here in order to make the present paper self-contained, and to provide additional explanations and examples.
3.1.1 Strain of a single link
Consider first a link of length , and apply to it an infinitesimal variation . Its relative extension, or infinitesimal strain, is , or equivalently tan03 . The "true strain" (also called "Hencky strain" tan03 ) is defined with respect to a state chosen as a reference (often a state without stress) using several equivalent expressions:
We perform these manipulations because the last expression of eq. (12) is the easiest to generalise. It is not a problem to take the log of dimensioned quantities (here, the square of a length) because this cancels out in the final result.
3.1.2 Statistical strain of the pattern
For a whole pattern, replacing by enables to perform statistical averages over links. The logarithm of is unambiguously defined and is easily performed in three standard steps on a computer (Appendix B.2). It suffices to first, switch to the three orthogonal axes (’s eigenvectors) in which is diagonal; second, take the logarithm of its eigenvalues, which are strictly positive (section 2.2.2):
and third, switch back to the original axes. It is necessary to perform first all linear operations such as averaging. This ensures in particular that all s in eq. (13) are non-zero. Taking the logarithm, which is a non-linear operation, has to be performed later.
Eq. (13), like eq. (12), requires to define a reference, expressed in the same units as , so that the difference of their logarithms is well defined and dimensionless. Such a reference texture is discussed in section 3.1.3.
The "statistical internal strain" is defined aub03 as:
Here completely characterises the material’s current strain: relative dilation, amplitude and direction of anisotropy.
3.1.3 Reference texture
Practical details regarding the reference texture are presented in appendix A.2.
Eq. (14) shows that the exact choice of affects the value of but not its variations. It thus does not appear explicitly in the kinematics (eqs. 21,80) nor in the dynamics (for instance in the value of the shear modulus, appendix A.3.3). Moreover, eq. (14) remains unchanged if we multiply both and by a prefactor; this is why the exact unit (e.g. m, mm, m) in which and are expressed is unimportant, as long as it is the same unit for both.
Whatever the choice, the reference is defined by the texture . It suffices to determine 6 numbers (3 numbers, if in 2D); or 1, in the (most common) case where is isotropic. Only the reference texture corresponding to the current state plays a role; past changes of the reference pattern, for instance during an irreversible strain (also called "work hardening" cha87 ), need not be taken into account. It is never necessary to know the details of the corresponding pattern’s structure, such as the positions of each object one by one: it is even not necessary that this pattern exists and is realisable.
This section presents a few examples and particular cases of internal strain.
(i) If the material is uniformly dilated (affine deformation, see section 4.3.1) by a factor in all directions, then . Thus , as is expected for instance for gases; that is, in 3D:
(ii) Conversely, if the material is uniformly dilated by a factor in one direction and compressed by a factor in another direction, then:
(iii) For incompressible materials, ’s diagonal terms are usually both positive and negative, and their sum is usually small: is mostly deviatoric (Appendix B.1). Note that even in incompressible materials the links’ mean square length can vary slightly, so that Tr() is not necessarily strictly zero. For instance, it reaches 0.03 in Fig. (1b) where a foam is sheared while keeping bubble number and total foam area exactly constant (Ataei Talebi and Quilliet, private communication).
(iv) In Fig. (9), most ellipses look circular; deviations from circles occur close to the obstacle. We distinguish regions where extension dominates, and ellipses are stretched like coffee beans, from regions where compression dominates, where the ellipses are flattened like capsules.
where is ’s eigenvalue (e.g. is if we use the definition of eq. 44). Eq. (15) reflects that and have the same eigenvectors: they commute. Eq. (15) also relates the trace of with ’s determinant (product of eigenvalues):
3.2 Kinematics: time evolution
3.2.1 Statistical velocity gradient: and
We want to define the continuous counterpart of the geometrical changes (eq. 10). We use (eq. 55), which is in m, and is always defined. For reasons which appear below (eqs. 27-30), we use (eq. 71) as an intermediate step, and define as:
has the dimension of a strain rate (s): its order of magnitude is the links’ average variation rate. Like , it vanishes when the pattern moves as a whole, with a rigid body translation.
In practice, the most useful quantity is its symmetric part, the "statistical symmetrised velocity gradient":
It is the rate of variation of due to the links’ stretching and relaxation. In cases where and commute, is symmetric and eq. (18) simply writes
Fig. (10) plots an example of . It is large all around the obstacle, but only very close to it; it is almost the same before and after the obstacle. When the material’s density is constant, Tr is small (but not necessarily exactly zero), and the corresponding ellipse is nearly (but not necessarily exactly) circular.
3.2.2 Statistical topological rearrangement rate:
Here we have introduced a factor so that is the term which unloads the statistical internal strain, as will appear in eqs. (21,80). In cases where and commute, such as in the companion paper Cartes2D , is symmetric and eq. (20) simply writes
This "statistical topological rearrangement rate" (eq. 20) has the dimension of s. It measures the frequency and direction of rearrangements: it is of the order of magnitude of the number of changes per unit time and per link. Corresponding ellipses are elongated like coffee beans (resp: flattened like capsules) if the number of links decreases (resp: increases); if the number of links is conserved, ellipses are nearly circular.
As an example, Fig. (11) shows that the rearrangements are more frequent just in front of the obstacle, or in a very narrow region behind it. The rate of rearrangements decreases smoothly with the distance to the obstacle. This is due to the foam’s elasticity. It contrasts with the sharp transition between solid-like and fluid-like regions observed in purely visco-plastic materials ber85 ; bla97 . The companion paper Cartes2D presents an example with a larger spatial distribution of topological events, which enables for a better spatial resolution.
3.2.3 Kinematic equation of evolution
We have thus three independent symmetric matrices: , and . As discussed in Appendix C.2.3, there is a relation between them.
That is, the (statistical) symmetrised velocity gradient is shared between two contributions: one part (which includes advection and rotation) changes the (statistical) internal strain, the other part is the (statistical) topological rearrangement rate. How is shared between both contributions constitutes the main subject of the companion paper Cartes2D .
4 Continuous mechanics
This section applies to materials which dynamics can be described using continuous mechanics (section 4.1), in terms of stresses. We examine whether it is possible to relate the continuous, large-scale, dynamical description on one hand; and on the other hand the statistical measurements based on discrete objects introduced in section 3, which describe the pattern’s connections (topology), shape (geometry) and movements (kinematics).
4.1 Continuous description and RVE
If the material acts as a continuous medium cha87 ; bat00 ; lan86 , it usually has the following properties. First, there exists a range of sizes over which measurements yield the same results dollet_local . In that case, the box is called a representative volume element (RVE). This is usually obeyed if is much larger than the range of interaction between individual objects, and also larger than the correlation length of their disorder (but these conditions are neither necessary nor sufficient). Second, its description can be local in space, that is, its equation of evolution involves partial space derivatives, and the spatial variations of its solutions look smooth. Third, the average quantities have at large scale a role more important than that of fluctuations.
Regarding the choice of the RVE, the discussion of section 2.1.2 applies. Here again, averages on detailed geometrical quantities are performed on a spatial box of volume and over a time selected to suit the problem under consideration. The shape of the box should preferably respect the system’s symmetries.
For the present purpose of a continuous description, there is however the additional requirement that . More precisely, the relative statistical uncertainty should be smaller than the relative precision required by the user. A few tens or hundreds of links are often enough (there is no need for links). This does not set any theoretical lower limit to the size of : it can well be as small as the link size, or even smaller, if there are enough images to average (Fig. 6c).
4.2 Elastic, plastic, fluid behaviours
If the pattern behaves as a continuous material, we can consider a RVE (section 4.1) at position . The velocity field is , that is, an average over the whole RVE. If and are the positions of two RVEs, the velocity gradient is the spatial derivative of the velocity field, and is its transposed (eqs. 60,61), then:
Details on this notation can be found in Appendix B.3. Eq. (22) neglects terms of order of and higher. It describes the velocity field as continuous and affine, that is, a term which varies linearly with position plus a constant term (offset).
One of the key ingredients of continuous mechanics is the velocity gradient’s symmetrical part, that is, the total strain rate:
This is a purely kinematical quantity, but it determines the contribution to the viscous (dissipative) stress bat00 .
For small strain (linear elastic regime), neglecting advection and rotation, the integration of eq. (23) defines a total applied strain, which is a function of the past history of the sample, as:
Here is the gradient of the displacement field , and its symmetrical part.
The total strain rate contributes in part (loading) to change the elastic strain , and in part (unloading) to a plastic strain rate which is defined by their difference:
Alternatively elasticity and plasticity are defined through dynamics. A given region of the pattern is said to be in elastic, plastic or viscous regime, according to the contribution to the stress that dominates locally pinceau . The elastic strain contributes to the reversible part of the stress. Plasticity describes the irreversible contribution to the stress in the low velocity limit (note that rearranging patterns can often deform a lot without breaking). Both are solid behaviours, that is, exist in the limit of very low velocity gradient. The viscous contribution to stress is irreversible: it is due to, and thus increases with, the velocity gradient; that is, relative movements of objects within the material.
4.3 Link between discrete and continuous descriptions
4.3.1 Affine assumption
The affine assumption is analogous to, but much stronger than, eq. (22). It assumes that the velocity of each individual object is affine too:
An affine flow field and a non-affine flow field are plotted on figure 12.
In other words, this affine assumption implies that the continous velocity gradient has a meaning down to the level of individual objects, and that fluctuations around it are small enough to have no effect on the material’s mechanical behaviour.
In cellular patterns, especially in dry ones where there are no gaps nor overlaps, the movement of each individual object is highly correlated with its neighbours’; thus the affine assumption is reasonable Cartes2D . In particle assemblies, it might apply to dense assemblies of repelling particles, which cannot be too close nor too far from each other.
Whenever this assumption is valid, it considerably simplifies the description of the pattern evolution. Consider for instance a link, (eq. 1). Its time derivative is
4.3.2 Velocity gradient and total strain rate
thus appears as a statistical measurement of the symmetrised velocity gradient . When only large scale measurements are possible, only can be measured. However, when the detailed information on links in available to perform statistics, measuring offers several advantages.
(i) The signal to noise ratio is optimal, in the sense that all the local information, and only it, is used. Each link acts as a small probe of the local velocity differences: the spatial derivative is taken naturally at the places where the objects are, not on the larger scale of RVEs.
(ii) is intrinsically based on the material’s structure. It can be defined and measured even if there are only a few objects; or if the standard deviation of their velocities is large. Averaging over all links provides a statistical measurement of the total strain rate. At no point does the definition or measurement of require any affine description.
(iii) Physically, we expect to play a more general role than , because it is based on the individual objects themselves. For instance, we expect to be determinant in yielding, and thus in the description of plasticity (and possibly too) Cartes2D . Similarly, the material’s internal dissipations are probably more closely related to changes in the links than to a large scale velocity gradient: this suggests that the dissipative contribution to the stress arises in general from rather than from .
4.3.3 Strain, in the elastic regime
In this section we consider the particular case where the material is in the elastic regime. There is no plastic strain rate, . Eq. (25) becomes simply:
More precisely, at least in the linear elastic regime, one can identify two quantities: the symmetrised gradient of the displacement field, , which is a function of the past history of the sample; and the elastic strain , which is a function of state. In fact, in elasticity, both quantities are considered as equivalent lan86 .
Thus the elastic strain can be measured using two different methods. When large scale measurements of total strain are possible, can be measured as . When the detailed information on links in available, measuring as offers many advantages, similar to that of (section 4.3.2).
An acceptable definition of strain must coincide with in the linear elastic regime. As a consequence, it also implies that it is a conjugate of stress: the scalar product of stress by an infinitesimal increment of strain equals the increment of energy. This is a dynamical constraint on acceptable definitions; but it is a weak constraint (especially since the conjugate equation is a scalar relation). In itself, it is insufficient to define all components of strain.
There are thus several families of acceptable definitions of internal strain jan05 ; Bagi2006 ; one family contains an infinity of acceptable definitions Farahani2000 . Some definitions are particularly adapted to a discrete pattern’s geometry kru03 or dynamics gol02 .
Here, eq.(14) is a definition of strain which is: (i) one of the definitions acceptable in the whole elastic regime, even at large strain (eq. 35) where it coincides jan05 with a true strain hog87 ; (ii) probably the only definition valid outside of the elastic regime jan05 , when bubbles rearrange and move past each other, that is, when the pattern flows: the main advantage of eq. (14) is that it does not require the detailed knowledge of each object’s past displacement.
4.3.4 Plastic strain rate, in steady flow
In the more general case, there is a plastic strain rate, , and deformations can be strongly non-affine. Eq. (32) does not hold. The current elastic strain and the total strain rate are independent physical quantities; can no longer be measured as (whether it can be measured as is discussed in section 4.3.5).
For instance, if the material flows, the displacement of an object relatively to its neighbours can be arbitrary large; can become much larger than . In the extreme cases of steady flows, independent of time, it is possible (in absence of advection) that , and eq. (25) reduces to:
According to eq. (20), a steady flow with a corotational derivative that vanishes (meaning no advection nor rotation effects, see eq. (81)), implies that all the geometrical strain rate translates into the topological strain rate:
Using the identification of eq. (31), we therefore obtain in that case:
4.3.5 Complete identification
The statistical tools and are always defined and measurable, even out of the elastic regime, or out of the steady regime. If we could identify them with and , respectively, it would make possible to measure the elastic strain in all regimes. This is certainly not possible in general, as shown by both following counterexamples referee .
In granular systems, due to solid friction in the contacts, irreversible plastic strains appear before the list of contacts changes. In solid networks (e.g. solid foams) with no topological change, the bond themselves might behave plastically, or they might perhaps undergo buckling instabilities leading to non-reversible stress-strain curves. Those are examples in which plasticity occurs before the first topological change.
Conversely, consider a set of rigid cables which resist tension, but no compression, and tie them together at knots to form a redundant, hyperstatic network. Under given external forces on the knots, some cables will be taut, others will dangle and transmit no force. Upon changing the forces, the list of taut, tension-carrying cables will change. This can be regarded as a topological change. The response, which implies displacements and strains, is however reversible and might be called elastic. Hence a case for which plasticity begins after the first topological change.
This identification might turn possible in some particular cases where one can express the stress as a function of kinematical quantities. This seems to be the case for foams and emulsions rau07 ; Cartes2D . We hope that in these cases, statistical measurements can constitute a coherent language to unify the description of elastic, plastic and fluid behaviours, as well as facilitate models and tests.
In the present paper, we define tools (Table 1) to extract information from a pattern made of discrete objects, subject to rearrangements, within a wide class of complex materials made of individual constituents such as atoms, molecules, bubbles, droplets, cells or solid particles. They characterise quantitatively the mutual arrangements of these objects, or more precisely the links between neighbouring objects.
Their definition, which can flexibly adapt to the questions to be answered, is operational. That is, given an experimental or simulated pattern, whether in 2D or 3D, there is a well defined method to measure them directly as statistics on individual constituents (links between neighbouring sites). This measurement is easy, and requires only a few basic operations on a computer: multiplication, average, diagonalisation, logarithm. It is robust to experimental noise, even if there is a limited number of links.
, (or ) and characterise the current state of the pattern, its geometrical changes, and its topological rearrangements, respectively. They are explicitly based on the pattern’s discrete structure. They can be measured locally, for instance on a single biological cell, or grain in crystals. But they can also be measured as averages over a larger region in space, or as time averages. In foams, measuring them smoothens out the pattern fluctuations due to the discrete nature of bubbles, and evidences the underlying behaviour of the foam as a continuous medium.
Their statistical counterparts , (or ) and , are independent of the pattern’s discrete length scale. Each of them exists and is valid together in elastic, plastic and fluid regimes : they unify the description of these three mechanical behaviours. They facilitate the comparison between experiments, simulations and theories. In at least the linear, affine, elastic regime, we suggest how to identify them with the quantities which characterise the continous mechanics: elastic strain, total strain rate, and plastic strain rate, respectively. From a practical point of view, this offers the advantage of measuring these continuous quantities with an optimal signal to noise ratio, even with few discrete objects. On a fundamental side, this provides a physical basis to the description of a continuous medium, at any local or global scale, by relating it to the individual constituents. Moreover, it provides a coherent language common to elasticity, plasticity and fluid mechanics.
The companion paper Cartes2D illustrates most of these points on a detailed practical example.
This work was initially stimulated by a seminar delivered by G. Porte. We thank M. Aubouy, S. Courty, J.A. Glazier, V. Grieneisen, M. Hindry, E. Janiaud, Y. Jiang, J. Käfer, S. Marée for discussions. We thank the colleagues who have made constructive comments about the first version of the manuscript.
Appendix A Measurement techniques
This Appendix, aimed at non-specialists, lists practical advices based on our past experience.
a.1 Averaging procedure
The average of any quantity is:
where the sum is taken over all links in the averaging region. Here is the weight of the link: for almost all links in the averaging region, ; at the boundaries of the averaging region, decreases to zero, different choices being possible (section A.1.2). Here we note:
For instance, the texture is:
a.1.2 Choices of weights
There are at least three main possible choices for the averaging procedure (Fig. 13). Once a procedure has been selected, it is important to keep consistently the same for all measurements.
The topology is useful for local information, especially for a single site: particle or cell. This is the case for instance when studying the division of a cell courtypreprint . In that case, each link is either included or excluded ("all or nothing", or 1). One should at least include the links between the site of interest and its neighbours (first shell). Statistics are over a few links only, and are easy to compute, even sometimes by hand. This defines the -th site’s texture as a sum over its neighbours, labelled :
One can also choose to include the second shell (next nearest neighbours), third, or even higher. The total pattern’s texture (eq. 40) appears as the average of the site textures, weighted by the site’s number of links, and with a factor because each link is counted twice (each link belongs to two sites):
The geometry is useful for a continuous description, typically to measure , or in a RVE as a function of space and time. This is the case for the examples of foam flow which illustrate this paper. The average is over all links in a box which respects as much as possible the symmetry of the problem: rectangle, annulus. The dimensions of the box determine the scale of averaging. The measurements are performed on boxes at different positions . The distance between measurements positions cannot be more distant than the box size (it would leave gaps between boxes), but they can be closer (thus boxes overlap). If the box dimension is much larger than a link, one might choose to neglect links which cross the box boundary. But in general, links which cross the boundary require attention (especially near the box corner) if an automatised image analysis is used. An additional choice is required. A first possibility ("all or nothing"), computationnally simpler, is to look where the link’s center lies: if the center lies inside the box, the link is assigned to this box (); else, this link is not counted () asi03 . A variant is to assign half of each link to the two boxes of the two bubble centres it binds dollet_local . A second possibility ("proportional"), which yields a better precision, consists in weighting the link with equal to the fraction of the link which is inside the box (the remaining is outside) jan05 .
The coarse graining is seldom convenient in the practical applications considered here. However, theoreticians use it gol02 , especially for the advection term (see Appendix C.1) dol05 . A link at position is counted in a box at position with a weight . The coarse graining function is a function which is non-increasing, from to , and has an integral equal to 1. Its width at half height (that is, where ) defines the scale of coarse graining. It is a continous and differentiable function, so that advected links smoothly enter and leave the averaging box, without singularity gol02 .
a.2 Choice of
Since the reference texture plays almost no physical role, its choice is not very important. It depends on the problem under consideration, but once its definition is chosen, it should be kept consistently. In practice, the choice depends on the available information. Here are a few possibilities.
(i) The most favorable case is when can be measured directly. In experiment, this is possible when an image can be chosen as reference, for instance a stress free pattern. In simulation (Fig. 3b), this requires to relax the stress under prescribed constraints.
(ii) can be determined theoretically in some cases, such as a set of particles which interaction potential is known. This occurs in Fig. (3a), where the natural reference is the honeycomb pattern with a link size:
Here is the area per link, and ; is the area per particle or cell (section A.3.1).
(iii) In cases such as Fig. (1), no reference state is known in details. If only is known, we suggest to take as isotropic. Although we do not know any fundamental reason for that, it seems to be satisfactory in all practical cases we have encountered. From eqs. (5, 8) it writes, in or 3 dimensions:
(v) The most unfavorable case is when even is unknown. A possibility is to take:
where is the average of the s, ’s eigenvalues. Taking the arithmetic average, , corresponds to the assumption that is conserved: . Taking the geometric average, , corresponds to the assumption that Tr, which is close to assuming that the material is incompressible (see section 3.1.4).
a.3 The case of 2D dry cellular patterns, especially foams
a.3.1 Number of neighbours
In 2D dry cellular patterns, the number of neighbours of each cell is variable; but its average over the whole pattern is always close to 6 neighbours, and thus 6 links, per cell wea99 ; weairerivier . Since each link is shared by two cells, the number of links is 3 times the number of cells. This is also true for a moderately wet cellular pattern, if neighbours are defined on a skeletonized image. It extends to Voronoi/Delaunay definition of neighbours for particles.
Some cells might meet by four ("4-fold vertex"). In that case, we recommend to decide that cells which share only a vertex should not be considered as neighbours. This choice is consistent with the fact that the texture describes the cell shape and arrangement (section A.3.3). Moreover, this avoids many artefacts when measuring the T1s.
a.3.2 Cell centers versus vertices
Aubouy et al. aub03 chose to describe a cellular pattern (such a 2D dry foam) as a network, each site being a vertex (that is, a point where three cells meet). Here, we prefer to use cell centers, for several reasons.
(i) First, and most important: the centers move according to the overall velocity field (while vertices have a highly fluctating displacement), thus the affine assumption (eq. 26) applies.
(ii) It is more robust, because a cell center is measured as an average over several pixels (while a vertex is a single pixel, which position might depend on the image analysis procedure).
(iii) This has the advantage of being more general: it applies to all other discrete patterns; and even within cellular patterns, it generalises to wet foams, and to 3D.
(iv) Finally, the topological rearrangements are well characterised (while on the opposite, if was defined as the vector between two vertices, a T1 occurence would be defined as , so that the contribution of a T1 to eq. (11) would systematically be zero).
a.3.3 Texture and inertia matrices
Marée et al. mar07 propose to measure the shear modulus by considering the variation of cell shapes. Each cell’s shape is characterised by its inertia matrix, : it looks similar to eq. (2), but it is averaged over the position of pixels inside a cell; thus their description is intra-cellular. Ours, averaged over the links between the cells, and thus based on the shape of the overall pattern, is rather inter-cellular.
In dry cellular patterns, where there are no gaps between cells, nor overlaps, the deformation of each cell is highly correlated to the global strain; thus, in this case, both descriptions coincide and yield approximately the same results courtypreprint .
For instance, note that the shear modulus is the variation of elastic stress with respect to infinitesimal variations of . This measurement is robust asi03 ; jan05 . As mentioned in section 3.1.3, it is not affected if we multiply and by a same prefactor; and even if we change , see for instance eq. (14). This is why this particular measurement gives similar results with both inertia and texture.
Here, we prefer to use the texture based on cell centers, which is more general, for several reasons.
(i) It also applies to characterise the strain of wet foams (where bubbles are round, and thus each bubble’s inertia is isotropic).
(ii) It applies to all other discrete patterns, including particle assemblies.
(iii) Centers, rather than shape, are involved in the kinematic description, including eqs. (18,20). It could in principle be possible to define an equivalent of (and even of ) based on inertia matrix, but its physical meaning is unclear; and it is probably not possible to define an equivalent of (and ).
(iv) It extends to more than one cell; while the inertia matrix of several cells can be defined, its physical meaning is not relevant to the pattern description.
Note that in the graphical representation of the inertia matrix, the ellipse axis lengths are the square root of the matrix’ eigenvalues mar07 . The advantage is that the ellipse elongation is the same as that of the actual cell. Here, taking the square root of eigenvalues has no physical signification for any matrix (except for the texture), so that we plot the matrix’ eigenvalues themselves (section 5).
Appendix B Matrices: notations and definitions
This appendix is aimed at readers who are not familiar with the matrices. We list all standard definitions used in the text, from the simplest to the most complicated.
We work here in a space with dimensions. A scalar is a simple number; a vector is a list of numbers; a matrix is an array of numbers. All these objects are tensors, of rank 0, 1 and 2, respectively. In this paper, there also appears , which is a tensor of rank 3 (for which there exists no particular name), that is, an array of numbers.
A matrix is an array with components , where the indices , 2 or 3:
Its trace is the sum of its diagonal terms:
Its transposed has components , and has the same trace. Any matrix can be rewritten as the sum of its symmetric and antisymmetric parts:
A matrix is said to be symmetric if it is equal to its transposed, , that is, (while an antisymmetric matrix is equal to minus its transposed); by definition, the symmetric part of is always symmetric. A symmetric matrix can itself be rewritten as an isotropic term and a traceless (or deviatoric) term:
where is the identity matrix in dimension (appearing in eq. 5).
Itself, the deviatoric part can be decomposed in diagonal components, called normal differences, and off-diagonal ones.
To summarize, a matrix has in general 9 independent components . They can be rewritten as 3 antisymmetric ones, namely , , and ; and 6 symmetric ones, namely 1 trace , 2 normal differences and , 3 off-diagonal terms ,