## Abstract

Under study is the geodesic (i.e. shortest path) character of strain fields occurring in inelastic response of matrix-inclusion composites. The spatially random morphology of composites is created by generating the inclusions centres through a sequential inhibition process based on a planar Poisson point field preventing any disc overlaps. Both phases (inclusions and matrix) are elastic–plastic hardening with the matrix being more compliant and weaker than the inclusions, and perfect bonding holding everywhere, so that the plastic flow occurs between the inclusions. A quantitative comparison of a response pattern obtained by computational mechanics with that found only by mathematical morphology indicates that (i) the regions of plastic flow are close (or even very close) to geodesics and (ii) a purely geometric (and orders of magnitude more rapid than by computational mechanics) assessment of these regions is possible.

## 1. Introduction

It is well known that homogeneous isotropic plastic materials display shear lines at 45° to the principal direction of stress and strain tensors. However, when a material's microstructure is considered, some perturbations away from the perfect 45° angle are observed—these are clearly due to spatial inhomogeneities that act as either obstacles or facilitators of plastic flow. In particular, this is observed in the random matrix-inclusion composites, where shear bands weave their way around the inclusions (e.g. Jiang *et al*. 2001; Doumalin *et al*. 2003; Li & Ostoja-Starzewski 2006). In effect, in any given (deterministic) realization of a random material, the shear bands display geometric patterns that conform to the actual spatial distribution of the material microstructure. Basically, the shear bands appear to follow paths of lowest plastic resistance while avoiding the obstacles of high-yield limit. Conceptually, therefore, they may well behave as *geodesics*, curves of shortest path joining two specific points in space.

Interestingly, over a decade ago, cracks forming and coalescing in random elastic-brittle media were found to follow the geodesics (Jeulin 1993). Following that work, Moulinec & Suquet (1998) in a work on numerical simulations mentioned a possibility for shear bands to follow the shortest paths, but did not check this point. A first application of geodesics to plasticity was proposed by Jeulin & Ostoja-Starzewski (2000; see also Ostoja-Starzewski 2005). More recently, trial fields related to geodesics (although not explicitly stated as such) were used to workout upper bounds in two-dimensional polycrystals (e.g. Goldsztein 2003). A concept more closely related to geodesics has recently been proposed to study strain localization in granular media (Unger 2007). In this paper, an analytical approach, somewhat modified from Jeulin's (1993) and building on that in Jeulin & Ostoja-Starzewski (2000), is applied to a two-dimensional composite material with elastoplastic hardening phases obeying an associated flow rule. First, using finite elements we compute an equivalent plastic strain field in the composite subjected to macroscopic shear. On the other hand, we obtain two-dimensional geodesics on the same composite, without accounting for any mechanics whatsoever—that is, for each and every point in the square-shaped material domain we determine the shortest distance to pass across the specimen (both from top to bottom and from left to right), while avoiding the inclusions. This allows us to assess the geodesic property of the strain field, and, if it turns out that the cross-correlation of both fields is positive, the approach based on geodesics is in orders of magnitude faster than the finite-element analysis.

## 2. Strain field patterns by computational mechanics

The composite body is taken as a random material in the sense that it is an ensemble of realizations *B*(*ω*), i.e. {*B*(*ω*), *ω*∈*Ω*}, where each *B*(*ω*) responds in a deterministic fashion. As commonly done in probabilistic theories and in accordance with the axiomatic definition of Kolmogorov (1933), *ω* denotes an elementary event from a sample space *Ω*, the latter being equipped with a *σ*-algebra and a probability measure *P* (for more details and examples in stochastic mechanics of materials, see Jeulin & Ostoja-Starzewski (2001) and Ostoja-Starzewski (2007)).

In general, is described by a random vector field ** Θ**={

*G*,

*ν*,

*c*}, where

*G*is the shear modulus;

*ν*is the Poisson ratio; and

*c*is the yield limit at a given point in space. In particular, we consider a two-phase elastic–plastic hardening matrix material with inclusions (

*p*=1) perfectly bonded to the matrix (

*p*=2) everywhere. The composites' spatially random morphology is created by generating the disc centres through a

*sequential inhibition process*based on a Poisson point field in plane. The constitutive response of

*B*(

*ω*) of is described by an associated flow rule(2.1)where primes indicate deviatoric tensor components. In (2.1)

*G*

_{p}is the shear modulus;

*ν*

_{p}is the Poisson ratio;

*f*

_{p}is the yield function; and

*c*

_{p}is the yield limit. Clearly,

*G*

_{p},

*ν*

_{p}and

*c*

_{p}form a vector whose each component (described by its indicator function

*Χ*

_{p}) gives rise to a scalar random field, such as(2.2)

We consider two special cases of such composite materials.

Elastic inclusions in the elastic–plastic hardening matrix, with the matrix having elastic response identical to that of the inclusions up to the yield point *σ*_{0} (figure 1*a*). For the matrix, the von Mises yield criterion is used and the rate-independent plasticity with associated flow rule and isotropic hardening is assumed. The matrix' stress–strain curve is characterized by a piece-wise power law (e.g. Dowling 1993)(2.3)with the yield stress *σ*_{0}=*Eϵ*_{0} and *ϵ*_{0} being the yield strain. We take *ϵ*_{0}=1/300, *σ*_{0}=170 MPa, *N*=0.1, and *E*=*σ*_{0}/*ϵ*_{0}=51 GPa and *ν*=0.3. Scale-dependent responses of such a material were studied by Jiang *et al.* (2001), with the figures in that reference clearly showing plastic flow occurring at roughly ±45° to the direction of principal stress and strain between the inclusions under various kinds of loading. We now apply uniform kinematic loading of simple shear type to the boundary of *B*(*ω*)(2.4)where , . The resulting pattern of equivalent (von Mises) plastic strain obtained by the computational mechanics for the particular realization of figure 2*a* is shown in figure 2*b*.

The inclusions and the matrix are elastic–plastic hardening and equation (2.1) applies to both phases. The matrix is more compliant and weaker than the inclusions (figure 1*b*). Values of material properties are given in table 1. Scale-dependent responses for case 2.2 were studied by Li & Ostoja-Starzewski (2006). Patterns of equivalent plastic strain, again obtained by computational mechanics under equation (2.4), are shown for three realizations *B*(*ω*) in the left column of figure 4.

## 3. Geodesics

### (a) Geodesics in the Euclidean space

Many applications in physics are based on some propagation phenomena, possibly with different propagation velocities for heterogeneous media. For example, the propagation of light in geometric optics (or acoustics), according to the principle of Fermat, in the setting of a heterogeneous medium, involves the shortest time paths. When the structure can be modelled by a valued planar graph, it is easy to compute distances to a source on the valued graph, usually called the *geodesic distance* (namely the length of shortest paths). This is done after defining the valuation of the edges of the graph according to the physical problem of interest (e.g. inversely proportional to the light velocity or proportional to the fracture energy when dealing with fracture processes) and after defining a source *W*_{s} and a destination *W*_{d} on images of the structure (for instance, two opposite edges of a two-dimensional field or two opposite faces of a three-dimensional field). The geodesic distance *d*(*V*,*W*_{s}) between *W*_{s} and any vertex *V* is then calculated. From this geodesic distance functions, useful information can be extracted.

The geodesic distance between

*W*_{s}and*W*_{d}(or the length of shortest paths linking the source and the destination); the statistics of this parameter can be investigated (Jeulin 1993).The distribution of the geodesic distances in the image. By application of thresholds to the image of distances, the sequence of propagation of the front is displayed.

The shortest paths connecting

*W*_{s}and*W*_{d}.

According to the geometric optics, the light rays follow paths directed by the Fermat principle, minimizing the travel time between two points. In a heterogeneous, locally isotropic medium with an index of refraction *n*(*x*) depending on the position *x*, the trajectories of rays starting from a source *W*_{s} to a destination *W*_{d} are made of geodesics paths. They can be found as a solution of the eikonal equation(3.1)where ** u** is a unit vector showing the direction of propagation of the light wave. Note for future reference when anisotropic composites will be studied that a generalization of the eikonal equation to locally anisotropic media has been given by Ostoja-Starzewski (2001).

The equation (3.1) can be solved by numerical techniques, using for instance the level sets and fast marching methods (Sethian 1999). In the present case, we employ digitized images of the microstructure and will therefore use the underlying graph structure to extract geodesic paths.

### (b) Graphs, paths and distances

We work on planar-valued graphs {*V*, *E*}, where *V* is the set of vertices (*V*_{i}) and *E* is the set of edges (*E*_{ij}) connecting pairs of nearest neighbours *V*_{i} and *V*_{j}. Each edge carries the value , which is the local metric.

Now consider two vertices which are generally not neighbours: *V*_{s} (a source) and *V*_{d} (a destination). We call a path *P*(*V*_{s},*V*_{d}) between the source and the destination a set of connected edges(3.2)The length *L*(*P*(*V*_{s},*V*_{d})) of that path is obtained through addition of all *d*(*V*_{i},*V*_{j})'s along it(3.3)

Note that for a graph representing a heterogeneous medium, the length associated with each edge depends on its length and of the component in which it is embedded. The distance *d*(*V*_{s},*V*_{d}) between *V*_{s} and *V*_{d} is the minimal value of *L*(*P*(*V*_{s},*V*_{d})) over all the possible paths *P*(*V*_{s},*V*_{d}) on the graph {*V*,*E*}.

Having the distance *d*(*V*_{s},*V*_{d}), it is straightforward to define two distance functions:

a distance from a vertex

*V*_{s}to a set of vertices*W*_{d}(3.4)where the symbol inf indicates the infimum over all the values*d*(*V*_{s},*V*_{d}) anda distance from a set of vertices

*W*_{s}to another set*W*_{d}(3.5)

It is clear from the above definitions that the source and the destination can be interchanged as long as there is no orientation of the edges. Otherwise, the orientation of the edges must be reversed for a backward propagation. In fact, this situation is very common in geometric optics, where use is made of direct and inverse propagations. We use this analogy to extract the geodesic (i.e. shortest) paths between two sets of vertices, *W*_{s} and *W*_{d}: the vertices belonging to the shortest path(s) satisfy the relationship(3.6)The equation (3.6) requires the use of two propagations, exchanging the roles of *W*_{s} and *W*_{d}.

### (c) Algorithms on a hexagonal graph

Different types of graphs can be used on images of a microstructure. A random graph can be generated from random points in the medium and the associated random Voronoï tesselation: we generate edges by connecting each point to the centre of the neighbouring cells in the tesselation. Operating on two-dimensional images, pixels can be arranged on a square or on a hexagonal grid. On the square grid, we can either connect every pixel to its four nearest neighbours or to its eight nearest neighbours, obtaining a 4- or 8-connected graph. On the hexagonal grid, a 6-connected graph is generated by connecting every pixel to its six nearest neighbours. Efficient algorithms for various types of graphs were developed by Jeulin *et al*. (1992) and Jeulin (1993).

In the present case, we operate on images of a two-phase medium (inclusions embedded in a matrix), where we expect that the shear flows weave around rigid inclusions which necessarily act as obstacles. We can therefore build an appropriate geodesic distance from the following construction: the distance *d*(*V*_{i},*V*_{j}) between two neighbouring vertices *V*_{i} and *V*_{j} is the unit distance if they belong to the matrix *D*; otherwise, we have *d*(*V*_{i},*V*_{j})=+∞.

To find the geodesics in this situation, basic geodesic operations of mathematical morphology (Serra 1982) are implemented. We work on a hexagonal grid, which is more isotropic than the square grid. Operating on a grid produces a digitized version of the geodesic distance that is well known to be a good approximation of the Euclidean case for images with a sufficient resolution, as in the present case (Serra 1982; Jeulin *et al*. 1992; Jeulin 1993). Starting from a binary image defining a set *A* (a pixel *x* belonging to *A* has the value 1, else 0), we define the dilated set by a structuring element *K* ( being the reflected set of *K* with respect to the origin: ) as the set of points *x* such as *K*_{x} (obtained by translation of *K* to the point *x*) intersects the set *A*(3.7)

Consider now, in addition to the set *A*, a set *D*. We define the unit geodesic dilation of set *A* in the mask *D* by the unit hexagon *H* as(3.8)By *n* iterations of the operation defined by equation (3.8), we obtain the geodesic dilation of set *A* of size *n*(3.9)By construction, the points belonging to the set are located at a geodesic distance of *A* lower than *n*. By summing in a grey level image each iteration , we obtain a valued image proportional to the negative of the geodesic distance to *A*. The distance increases with *n*, but for images of finite size it is bounded; the algorithm of equation (3.9) is stopped when the upper bound of the distance is reached. To generate the geodesic paths avoiding inclusions, we use two propagations according to equation (3.6). We start with a source image *W*_{s} (here one edge of the image) for the set *A*, the opposite edge *W*_{d} as destination and the image of the matrix as set *D*. The backward propagation is obtained after interchanging *W*_{s} with *W*_{d}.

Note that an opposite situation concerns crack propagation in a homogeneous matrix. In that case, geodesic distances are defined in such a way that cracks (or pores) connected to the current front of propagation with size *n* are added to the front, the ‘distance’ on edges in pores carrying the value 0 (Osmont *et al*. 1987; Vincent & Jeulin 1989). They are also implemented by iterations of binary morphological operations.

## 4. Geodesics vis-à-vis strain fields

Now, in accordance with the ‘shortest path’ hypothesis, the geodesics should join the opposite faces of a specimen *B*(*ω*), such as that in figure 1*a*, to which the simple shear loading is applied. Therefore, as potential shear bands we consider two families of shortest paths obtained by geodesic propagations in two orthogonal directions that avoid (black) inclusions. For a horizontal propagation as source *W*_{s} and destination *W*_{d} we first take the left and right faces and then invert their roles obtaining the set *G*-hor. For a vertical propagation as source *W*_{s} and destination *W*_{d} we first take the top and bottom faces and then invert their roles obtaining the set *G*-ver. In fact, in order to obtain similar graphs for the ‘horizontal’ and ‘vertical’ propagations, we perform horizontal propagations on the image after a *π*/2 rotation, and we rotate back the resulting map of geodesic distances. We consider propagations obtained on hexagonal lattices, described in §3*c*. Final results for the material of case 2.1, following Jeulin & Ostoja-Starzewski (2000), are as follows:

addition of

*G*-hor with*G*-ver to obtain the set*G*-add (figure 3*a*) andsupremum of

*G*-hor with*G*-ver to obtain the set*G*-sup (figure 3*b*).

The greyscale at any given point in figure 2*a*,*b* indicates the (appropriately normalized) distance it takes to connect the opposite edges of the square-shaped domain *B*_{δ}(*ω*): the darker the point, the shorter the distance. Given two fields such as (geodesic) and (equivalent plastic strain), both defined on a square lattice of the same size (256×256), so that *I*=256^{2}, we compare them by using the cross-correlation coefficient(4.1)where *C*_{ge} is the cross-correlation of *g* and *e*. Recall that, in general, . Basically, working with a single realization forces us to assume both fields, *g* and *e*, to be correlation ergodic. Therefore, we estimate *ρ*_{ge} by the following expression for all points of the matrix (outside of inclusions which do not participate in the plastic flow at all):(4.2)It turns out that the cross-correlation of *G*-sup (figure 3*b*) with the true solution, *e*, obtained by finite elements (figure 2*b*) is only approximately 0.2, and a lower number (approx. 0.1) is obtained for the *G*-add geodesic propagation (figure 3*a*).

The situation changes when we move to the material of case 2.2 (see figures 1*b* and 4). Using the same methodology as above, with *g* computed using *G*-sup and including the inclusions which now do participate in the plastic flow, we find *ρ*_{ge} to vary between 0.78 and 0.92. This is a very high cross-correlation!

## 5. Conclusions

A quantitative comparison of strain patterns obtained by computational mechanics with those found simply by mathematical morphology indicates that the regions of plastic flow are very close to geodesics for a composite with both matrix and inclusions being elastic–plastic hardening and only approximately so, when the inclusions are elastic. The goodness of the geodesic approximation is seen to stem from the difference in initial moduli of matrix and inclusions (composite of case 2.2), which allows more deformation, and then plastic flow, to develop in the gaps between the inclusions. In contradistinction to this, the composite of case 2.1 is allowed (as well as forced) to build up stress and strain fields uniformly throughout the entire domain along the common stress–strain curve. By the time the yield point of the soft (i.e. matrix) phase is exceeded, the shear bands in that phase tend to develop along more curved paths than the geodesics, and hence the cross-correlation of both fields is poorer.

Interestingly, the geodesics can be computed much more rapidly than by, say, finite elements—in a matter of seconds on a laptop computer as opposed to about an hour. This may offer a very rapid way of a purely geometric assessment of zones of plastic flow in disordered heterogeneous materials, although much further research needs to be done to assess the appropriate ranges of all the material parameters and to demonstrate mathematically that fields of von Mises plastic strain are close to, but not identical with, the geodesics. Clearly, the latter aspect is a major challenge since the geodesics are scalar fields while the equivalent plastic strain results from a second-rank tensor field.

## Footnotes

- Received November 18, 2007.
- Accepted January 17, 2008.

- © 2008 The Royal Society