Propagation of Elastic Waves in Homogeneous Media: 2D Numerical Simulation for a Concrete Specimen

Abstract

This paper addresses the theoretical foundation of a localization method for crack detection in a concrete sample based on the time of arrival of the elastic wave generated by the crack formation to a group of sensors positioned on the boundary of the sample. The equations of motion for the elastic waves are carefully presented, including a body force term which accounts for the sudden formation of a crack. Then, a localization method based on the detection of acoustic emissions, and specifically on their arrival times, is described. Finally, a discretization scheme for the 2D equations of elasticity is developed, and some numerical experiments are performed to assess the validity of the method.

1. Introduction

In recent years, many university researchers have studied the problem of concrete damage through the use of Acoustic Emission [1,2,3,4,5]. Acoustic emission (AE) is a natural phenomenon that occurs in the presence of a crack, or a dislocation source, in a material. The response of concrete to the load to which it is subjected in both compression and tension is affected by cracking. Researchers first investigated the microscopic behavior of concrete under compressive stress. The findings in [6,7,8] show that the stress–strain response of concrete is closely associated with the formation of microcracks, that is, cracks that form at coarse-aggregate boundaries (bond cracks) and propagate through the surrounding mortar (mortar cracks), as shown in Figure 1.

During the early studies on microcracking, concrete was considered to be composed of two brittle, linear and elastic materials: the paste and cementitious aggregates; moreover, the main causes of the non-linear stress–strain behavior of concrete in compression were considered to be microcracks [6,8]. This approach changed in the 1970s. Cement is a non-linear softening material, as is the mortar that makes up the concrete. The nonlinear compressive behavior of concrete depends largely on the response of these two materials and depends less than originally thought on the bond and microcracks of the mortar [9,10,11]. However a significant portion of the nonlinear deformation of mortar and cement is due to the formation of microcracks several orders of magnitude smaller. The effect of macroscopic cracks on the performance and failure characteristics of concrete is also an important factor. They reduce the characteristic resistance of the material and are formed as a union of microcracks in a given neighborhood. In the analysis of the AE signals generated in concrete, it is assumed that each fracture consists of the sum of many microcracks that join to compose a larger fracture. The largest crack will be characterized by the maximum amplitude of the AE signal. The signal is recorded by special acoustic emission sensors operating at precise resonance frequencies and in the ranges between 50 and 400 kHz [12,13,14]. These frequencies represent the range in which the signal develops within the concrete under stress.

Different methods for studying and analyzing these signals can be found in the literature. One of the most recent and interesting methods is the one that links the b value to AE signals. This method is based on the Gutenberg Richter law (GBR), typically used in seismology to study earthquakes [15,16]. The GBR law defines the relationship between the magnitude and the total number of recognized earthquake events in a region during a predetermined time interval. The b-value parameter, defined within the framework of the GBR law, is used to select AE signals that identify critical damage events [17,18]. Specifically, critical AE signals are selected provided that the b-value has a value in the neighborhood of 1. In fact, in this neighborhood, the AE signals have the highest amplitude and are generated by the most important damage events. This extension of a law of seismology to AE was possible because of the strong similarities between acoustic waves and seismic waves. In fact, both types of waves can be classified as damped waves, characterized by a peak and a progressive damping of the signal. This technique is used to detect major cracks that develop in concrete during different types of tests, such as the compressive test or the three-point test. To date, the interest of the scientific community has focused on the use of this and other methods to study composite materials and concrete under flexural stresses in order to facilitate the identification of the area where damage develops. The study of a brittle material, such as concrete, is more challenging because the entire section of the compression-proof material is subjected to stress, making the area early damage development uncertain. In addition, we briefly mention that another aspect that could influence the propagation of AE in a porous material, such as concrete, is the aggression of the material by some chemical pollutant [19,20].

Another interesting issue is the localization of the damage which is the source of AE. This is one of the most important criteria for damage detection and classification. In [21], the authors investigated a microseismic/acoustic (MS/AE) emission source localization method to predict and control the formation of potentially dangerous fractures in complex structures. In order to avoid localization errors induced by both the irregular shape of concrete structures and the propagation velocity of elastic waves, a velocity-free MS/AE source localization method was developed. The method avoids repetitive manual training by using equidistant grid points for path finding, introducing the A* search algorithm and using grid points to fit complex structures and velocity-free source localization. Other localization methods not aimed at damage localization in fragile materials such as concrete but useful for improving the accuracy of indoor localization are based on the integration of inertial navigation system (INS) with wireless sensor network (WSN) [22].

Several studies have been carried out by researches in relation to AE source localization and are available in the literature [23,24]. These different methods are characterized by uncertainty about the location of cracks and need to be improved for brittle materials subject to compression. They often depend on the wave propagation velocity in the material, and no index exists in the literature for materials such as concrete, but there are uncertain ranges from 5000 to 10,000 m/s. A possible approach to validating and seeking damage in concrete is based on the use of X-ray tomography. This technology is still in the experimental stage, and it could provide 3D reconstructions of the specimens before and after testing, allowing the exact location of cracks to be obtained, with results comparable to those obtained with localization methods. Concrete is an X-ray shielding material, and so far, the experimentation has only been applied to small specimens [25].

This paper addresses the theoretical foundation of a localization method for crack detection in a concrete sample based on the time of arrival of the elastic wave, generated by the crack formation, to a group of sensors positioned on the boundary of the sample. The following Section 2 is devoted to writing the equation of motion for the elastic wave, with a body force term which accounts for the sudden formation of a crack. A localization method is described on Section 3, and, finally, some numerical experiments are performed in Section 4 to assess the validity of the method.

2. Equations of Motion

We consider a sample of concrete, which is subject to pressure stress along the vertical direction. The sample is not homogeneous, but it is possible to extract effective elastic properties, so that it can be modeled as an elastic material with Lamé constants $\lambda$, $\mu$. The compression stress will eventually produce one or multiple cracks inside the sample, which will cause the propagation of elastic waves. The waves can be detected by some sensors located on the boundary of the sample. Here, the relevant information is the time of detection of the wave arrival at each sensor. The problem that we want to solve is the identification of the position of the cracks by using only the above information. Together with the position, we also want to determine the time of the cracks’ formation.

Before discussing the identification procedure, we present the mathematical model of the direct problem, that is, the propagation of an elastic wave in a concrete sample, originated by a single pointsize crack formation inside it at a given time. We give a brief summary of the classical theory, following [26]. Let us consider a displacement field $\mathbf{u}(\mathbf{x},t)\in {\mathbb{R}}^3$, with space coordinates $\mathbf{x}\in \mathsf{\Omega }\subset {\mathbb{R}}^3$ and time $t>0$. Here, $\mathsf{\Omega }$ represent the region occupied by the sample. The displacement field satisfies the dynamic elastic equilibrium equations:

$$\rho \frac{{\partial }^2\mathbf{u}}{\partial t^2}=\nabla ·\mathbf{\sigma }+\mathbf{f},$$

(1)

where $\rho$ is the density, depending on space; $\mathbf{\sigma }$ is the stress tensor; and $\mathbf{f}$ is a body force term, depending on space and time. For an isotropic medium, Hooke’s law provides a linear constitutive relation for the stress tensor $\mathbf{\sigma }$ in terms of the strain tensor $\mathbf{\epsilon }:=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)$, that is,

$$\mathbf{\sigma }=\lambda \mathrm{Tr}\left(\mathbf{\epsilon }\right)\mathbf{I}+2\mu \mathbf{\epsilon },$$

(2)

where $\lambda$ and $\mu$ are Lamé’s first and second parameters, respectively. They can be expressed in terms of the Young modulus E and the bulk modulus K by $\lambda =\frac{3K(3K-E)}{9K-E}$ and $\mu =\frac{3KE}{9K-E}$.

In component form, let $\mathbf{u}=(u_1,u_2,u_3)$, $\mathbf{x}=(x_1,x_2,x_3)$, and the notation $u_{a,t}$ denote the derivative of $u_a$ with respect to t, the notation $u_{a,b}$ denote the derivative of $u_a$ with respect to $x_b$, $a,b=1,2,3$. Using Einstein’s convention of implicit summation for repeated indices, we have:

$${\mathbf{\epsilon }}_{ab}=\frac{1}{2}(u_{a,b}+u_{b,a}),{\sigma }_{ab}=\lambda {\mathbf{\epsilon }}_{cc}{\delta }_{ab}+2\mu {\mathbf{\epsilon }}_{ab}\equiv \lambda u_{c,c}{\delta }_{ab}+\mu (u_{a,b}+u_{b,a}),$$

where ${\delta }_{ab}$ is the Kronecker delta, and Equation (1) becomes:

$$\rho u_{a,tt}=(\lambda +\mu )u_{c,ca}+\mu u_{a,cc}+f_a,a=1,2,3.$$

(3)

The body force components $f_a$ will be chosen to model a sudden separation of two sides of a fault localized in a point, which might be thought of as an infinitesimal planar crack with a prescribed normal direction. Following the classical paper by Burridge and Knopoff [27], let us assume that the displacement and its derivatives present a discontinuity across a surface $\mathsf{\sum }$ embedded in $\mathsf{\Omega }$. Let $\mathbf{n}$ denote the normal direction to $\mathsf{\sum }$, and $\left[u_a\right](\mathbf{x},t)$, $\left[u_{a,b}\right](\mathbf{x},t)$ the discontinuity of $u_a$, $u_{a,b}$ across $\mathsf{\sum }$ in the direction $\mathbf{n}$. Then, the body force equivalent is given by:

$$\begin{array}{cc}\hfill f_a(\mathbf{x},t)=& -{\int }_{\mathsf{\sum }}{n}_d\{\left[u_c\right](\mathbf{\sigma },t)c_{cdab}{\delta }_{,b}(\mathbf{x}-\mathbf{\sigma })\hfill \\ \hfill & +\left[u_{c,b}\right](\mathbf{\sigma },t)c_{adcb}\delta (\mathbf{x}-\mathbf{\sigma })\}\mathrm{d}{\mathsf{\sum }}_{\mathbf{\sigma }},\hfill \end{array}$$

(4)

where $\delta$ is the Dirac delta function and ${\delta }_{,b}$ is its (distributional) derivative with respect to $x_b$. The occurrence we want to describe is explicitly considered in [27], in the special case when $\mathsf{\sum }$ collapses to the origin and $\mathbf{n}={\mathbf{e}}_3\equiv (0,0,1)$. Then, we have $\left[u_1\right]=0$, $\left[u_2\right]=0$, $\left[u_3\right]=H\left(t\right)\delta \left(x_1\right)\delta \left(x_2\right)$, where $H\left(t\right)$ is the Heaviside function, and the equivalent body force has components:

$$\begin{array}{ccc}\hfill f_1(\mathbf{x},t)& =& -\lambda H\left(t\right){\delta }_{,1}\left(x_1\right)\delta \left(x_2\right)\delta \left(x_3\right),\hfill \\ \hfill f_2(\mathbf{x},t)& =& -\lambda H\left(t\right)\delta \left(x_1\right){\delta }_{,2}\left(x_2\right)\delta \left(x_3\right),\hfill \\ \hfill f_3(\mathbf{x},t)& =& -(\lambda +2\mu )H\left(t\right)\delta \left(x_1\right)\delta \left(x_2\right){\delta }_{,3}\left(x_3\right).\hfill \end{array}$$

In compact form, introducing the diagonal matrix $D=\mathrm{diag}(\lambda ,\lambda ,\lambda +2\mu )$, we can write:

$$\left[u_a\right](\mathbf{x},t)={\delta }_{a3}H\left(t\right)\delta \left(x_1\right)\delta \left(x_2\right),f_a(\mathbf{x},t)=-H\left(t\right)D_{ab}{\delta }_{,b}\left(\mathbf{x}\right).$$

(5)

In the general case, we need to transform to new coordinates $x_a^{\prime }$ with respect to a basis ${\mathbf{e}}_a^{\prime }$, such that ${\mathbf{e}}_3^{\prime }=\mathbf{n}$. We then apply the previous result, and finally, we perform the inverse transformation. Introducing the rotation matrix R, with $R^{-1}=R^T$, defined by $x_a=R_{ab}{x}_b^{\prime }$, $n_a=R_{a3}$, we have:

$$\left[u_a\right](\mathbf{x},t)=n_{a}H\left(t\right)\delta \left(R_{b1}{x}_b\right)\delta \left(R_{b2}{x}_b\right),f_a(\mathbf{x},t)=-H\left(t\right)R_{ab}{D}_{bc}{R}_{dc}{\delta }_{,d}\left(\mathbf{x}\right).$$

(6)

Here, we have used the inverse transformation $x_a^{\prime }=R_{ba}{x}_b$ and the property $\delta \left(R{\mathbf{x}}^{\prime }\right)=\delta \left(\mathbf{x}\right)$. The body force term can be translated to an arbitrary point ${\mathbf{x}}_0$ by replacing $\mathbf{x}$ with $\mathbf{x}-{\mathbf{x}}_0$ inside the Dirac delta functions.

Summing up, we model the formation of a crack in a concrete sample by Equation (3), with

$$f_a(\mathbf{x},t;\mathbf{n},{\mathbf{x}}_0)=-H\left(t\right)R_{ab}\left(\mathbf{n}\right)D_{bc}{R}_{dc}\left(\mathbf{n}\right){\delta }_{,d}(\mathbf{x}-{\mathbf{x}}_0),$$

(7)

where the parameter $\mathbf{n}$ is the normal direction to the crack and comprises three angular variables in three dimensions (Euler angles) or an angular variable in two dimensions, and the parameter ${\mathbf{x}}_0$ is the location of the crack, that is, three coordinates in three dimension, or two coordinates in two dimensions, for a total of six parameters in three dimensions, or three parameters in two dimensions. These parameters are treated as random variables. The equations are supplemented with homogeneous Neumann boundary conditions on the boundary of the space domain and with zero initial data.

Once we have a model for the propagation of the elastic wave generated by the formation of a point crack, we need to address another important element for the definition of the localization problem, described in the following section, that is, the velocity of propagation of the elastic wave. To find the characteristic velocities, we look for a plane wave solution for (3), with zero body force, of the form:

$$u_a(\mathbf{x},t)=U_a{e}^{i(\mathbf{k}·\mathbf{x}-\omega t)},\mathbf{k}\in {\mathbb{R}}^3,\omega \in \mathbb{R},$$

(8)

where $U_a$ and $a=1,2,3$ are the amplitudes, $\mathbf{k}$ is the wavenumber vector, and $\omega$ is the angular frequency. Plugging (8) in Equation (3) with $f_a=0$, we find the eigenvalue problem:

$$[(\lambda +\mu )k_a{k}_b+{\mu \left|\mathbf{k}\right|}^2{\delta }_{ab}-\rho {\omega }^2{\delta }_{ab}]U_b=0,a=1,2,3,$$

(9)

for the eigenvalue ${\omega }^2$ and eigenvector $\mathbf{U}=(U_1,U_2,U_3)$. The solutions are:

$$\begin{array}{c}{\omega }^2=\frac{\lambda +2\mu }{\rho }{\left|\mathbf{k}\right|}^2=:c_p^2{\left|\mathbf{k}\right|}^2,\mathbf{U}=\mathbf{k},\hfill \\ {\omega }^2=\frac{\mu }{\rho }{\left|\mathbf{k}\right|}^2=:c_s^2{\left|\mathbf{k}\right|}^2,\mathbf{U}=\left[\begin{array}{c}-k_2\\ k_1\\ 0\end{array}\right],\left[\begin{array}{c}-k_3\\ 0\\ k_1\end{array}\right],\hfill \end{array}$$

corresponding to dilatational and isochore waves, with velocities $c_p$ and $c_s$, respectively. For all the considerations in the next section, we will use the primary velocity $c_p$ to evaluate the ratio between the distance of a sensor from the source of the crack and the corresponding time of travel.

3. Crack Detection

As discussed in the previous section, the location and orientation of a sudden crack are described by the vectors ${\mathbf{x}}_0$ and $\mathbf{n}$. We assume to have N sensors, at the positions ${\mathbf{x}}_i={(x_i^1,x_i^2,x_i^3)}^T$, with $i=1,2,\dots ,N$. The i-th sensor detects the arrival of an elastic wave at time $t_i$, $i=1,\dots ,N$. We assume for simplicity that a crack only develops at the position ${\mathbf{x}}_0={(x_0^1,x_0^2,x_0^3)}^T$ at time $t_0$. As a working hypothesis, we also assume that the orientation of the crack is not relevant for the effectiveness of the localization procedure that we will present. We aim to identify ${\mathbf{x}}_0$ and $t_0$ from the reading of the arrival times $t_i$, $i=1,\dots ,N$.

We know that the speed of the elastic wave in the direction of the front propagation is $c_p$, so we have:

$$c_p(t_i-t_0)=\parallel {\mathbf{x}}_i-{\mathbf{x}}_0\parallel ,i=1,2,\dots ,N.$$

(10)

To solve for $({\mathbf{x}}_0,t_0)$, we need to minimize the residuals

$$r_i=c_p(t_i-t_0)-\parallel {\mathbf{x}}_i-{\mathbf{x}}_0\parallel ,i=1,2,\dots ,N,$$

that is, to minimize the function given by their quadratic norm,

$$S:=\sum _{i=1}^N{r}_i^2=\sum _{i=1}^N(c_p(t_i-t_0)-\parallel {\mathbf{x}}_i-{\mathbf{x}}_0{\parallel )}^2.$$

(11)

This is a nonlinear least squares problem which can be solved by the Gauss–Newton method. Before discussing this method, we present a strategy for choosing a good initial guess by means of the simple intersection algorithm. The starting point is Equation (10), which, after squaring both members, can be written in the form:

$$2{\mathbf{x}}_i^T{\mathbf{x}}_0-2c_p^2{t}_i{t}_0=\parallel {\mathbf{x}}_i{\parallel }^2-c_p^2{t}_i^2+{\parallel {\mathbf{x}}_0\parallel }^2-c_p^2{t}_0^2.$$

(12)

Introducing the vectors ${\mathbf{X}}_i={(x_i^0,x_i^1,x_i^2,x_i^3)}^T$, $x_i^0=c_p{t}_i$, $i=0,1,\dots ,N$, and the matrix $\eta =\mathrm{diag}\{-1,1,1,1\}$, we can write (12) as:

$$\left[\begin{array}{c}2{\mathbf{X}}_1^T\eta \\ ⋮\\ 2{\mathbf{X}}_n^T\eta \end{array}\right]{\mathbf{X}}_0=\left[\begin{array}{c}{\mathbf{X}}_1^T\eta {\mathbf{X}}_1\\ ⋮\\ {\mathbf{X}}_n^T\eta {\mathbf{X}}_n\end{array}\right]+\left[\begin{array}{c}1\\ ⋮\\ 1\end{array}\right]{\mathbf{X}}_0^T\eta {\mathbf{X}}_0,$$

that is:

$$A{\mathbf{X}}_0=\mathbf{b}+\mathbf{1}{\mathbf{X}}_0^T\eta {\mathbf{X}}_0.$$

(13)

Finally, subtracting the first equation from the remaining $N-1$ equations, we can eliminate the nonlinear terms on the right-hand side and obtain an estimated value of ${\mathbf{X}}_0$:

$${\mathbf{X}}_0={\left(DA\right)}^{†}D\mathbf{b},D=\left[\begin{array}{cc}-\mathbf{1}& I_{N-1}\end{array}\right]\in {\mathbb{R}}^{N-1\times N},$$

(14)

where the dagger denotes the Moore–Penrose pseudo-inverse of a matrix, that is:

$$M^{†}={\left(M^{T}M\right)}^{-1}{M}^T.$$

Now, we can go back to the Gauss–Newton method, replacing the original least squares problem for (11) with a sequence of linear least square problems. We write (11) as:

$$S=\sum _{i=1}^N{r}_i^2=\sum _{i=1}^N{(X_i^0-f_i\left({\mathbf{X}}_0\right))}^2,f_i\left({\mathbf{X}}_0\right)=X_0^0+\parallel {\mathbf{x}}_i-{\mathbf{x}}_0\parallel ,$$

(15)

and the minimum occurs when the gradient of S with respect to ${\mathbf{X}}_0$ is zero:

$${\nabla }_{{\mathbf{X}}_0}S=-2\sum _{i=1}^N(X_i^0-f_i\left({\mathbf{X}}_0\right)){\nabla }_{{\mathbf{X}}_0}{f}_i\left({\mathbf{X}}_0\right)=\mathbf{0},$$

(16)

with

$${\nabla }_{{\mathbf{X}}_0}{f}_i\left({\mathbf{X}}_0\right)=\left[\begin{array}{cc}1& -\frac{{\mathbf{x}}_i^T-{\mathbf{x}}_0^T}{\parallel {\mathbf{x}}_i-{\mathbf{x}}_0\parallel }\end{array}\right].$$

(17)

Assuming we know an approximation ${\mathbf{X}}_0^{\left(k\right)}$ of the solution for $k\ge 0$, we set:

$${\mathbf{X}}_0^{(k+1)}={\mathbf{X}}_0^{\left(k\right)}+\mathsf{\Delta }{\mathbf{X}}_0^{\left(k\right)},$$

(18)

and we solve for $\mathsf{\Delta }{\mathbf{X}}_0^{\left(k\right)}$ using the following approximate gradient equations:

$$-2\sum _{i=1}^N\left(X_i^0-f_i\left({\mathbf{X}}_0^{\left(k\right)}\right)-{\nabla }_{{\mathbf{X}}_0}{f}_i\left({\mathbf{X}}_0^{\left(k\right)}\right)\mathsf{\Delta }{\mathbf{X}}_0^{\left(k\right)}\right){\nabla }_{{\mathbf{X}}_0}{f}_i\left({\mathbf{X}}_0^{\left(k\right)}\right)=0.$$

(19)

In compact form, introducing

$$\mathbf{r}\left({\mathbf{X}}_0\right)=\left[\begin{array}{c}{X}_1^0-f_1\left({\mathbf{X}}_0\right)\\ ⋮\\ X_n^0-f_n\left({\mathbf{X}}_0\right)\end{array}\right]\in {\mathbb{R}}^N,J\left({\mathbf{X}}_0\right)=\left[\begin{array}{c}{\nabla }_{{\mathbf{X}}_0}{f}_1\left({\mathbf{X}}_0\right)\\ ⋮\\ {\nabla }_{{\mathbf{X}}_0}{f}_n\left({\mathbf{X}}_0\right)\end{array}\right]\in {\mathbb{R}}^{N\times 4},$$

(20)

we can write:

$$J{\left({\mathbf{X}}_0^{\left(k\right)}\right)}^T\left(\mathbf{r}\left({\mathbf{X}}_0^{\left(k\right)}\right)-J\left({\mathbf{X}}_0^{\left(k\right)}\right)\mathsf{\Delta }{\mathbf{X}}_0^{\left(k\right)}\right)=\mathbf{0},$$

(21)

which yields the gradient equations:

$$\left(J{\left({\mathbf{X}}_0^{\left(k\right)}\right)}^{T}J\left({\mathbf{X}}_0^{\left(k\right)}\right)\right)\mathsf{\Delta }{\mathbf{X}}_0^{\left(k\right)}=J{\left({\mathbf{X}}_0^{\left(k\right)}\right)}^T\mathbf{r}\left({\mathbf{X}}_0^{\left(k\right)}\right).$$

(22)

In conclusion, the proposed localization method consists of an initialization step, using expression (14), and an iterative method defined by Equations (18) and (22).

4. Numerical Simulations

To validate the localization method developed in the previous sections, we perform a two-dimensional numerical simulation of Equation (3), with a body force given by the expression (7). Then, the sensors will be represented by a set of points on the boundary, and the reading of a specific sensor by the the first time at which the displacement vector at that point becomes different from zero.

We consider the two-dimensional version of system (3) in a rectangular sample $\mathsf{\Omega }=[0,l_x]\times [0,l_y]$, with homogeneous Neumann boundary conditions and zero initial value:

$$\left\{\begin{array}{c}\rho u_{tt}={[(\lambda +2\mu )u_x+\lambda v_y]}_x+{[\mu v_x+\mu u_y]}_y+f,\hfill \\ \rho v_{tt}={[\mu v_x+\mu u_y]}_x+{[\lambda u_x+(\lambda +2\mu )v_y]}_y+g,\hfill \\ \mathrm{for}(x,y)\in \mathsf{\Omega },t>0,\hfill \\ u_x=v_x=0,\mathrm{for}x=0,l_x,y\in [0,l_y],t>0,\hfill \\ u_y=v_y=0,\mathrm{for}x\in [0,l_x],y=0,l_y,t>0,\hfill \\ u=v=u_t=v_t=0,\mathrm{for}(x,y)\in \mathsf{\Omega },t=0.\hfill \end{array}\right.$$

(23)

Here, for simplicity, we do not use indices for the dependent and independent variables, so ${(u,v)}^T$ is the displacement vector, depending on space variables $\mathbf{x}={(x,y)}^T$ and time t, and the body force is ${(f,g)}^T$. Let $\theta$ be the angle formed by the infinitesimal crack line $\mathsf{\sum }$ with the x-axis, so that the tangent vector to $\mathsf{\sum }$ is ${(cos\theta ,sin\theta )}^T$, and the normal vector is $\mathbf{n}={(-sin\theta ,cos\theta )}^T$, and let ${\mathbf{x}}_0={(x_0,y_0)}^T$ be the location of the crack. Then, the rotation matrix R and the diagonal matrix D are given by

$$R=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right),D=\left(\begin{array}{cc}\lambda & 0\\ 0& \lambda +2\mu \end{array}\right),$$

and the expression (7) which gives the corresponding body force becomes

$$\left(\begin{array}{c}f\\ g\end{array}\right)=-H\left(t\right)(\lambda I+2\mu \mathbf{n}{\mathbf{n}}^T)\left(\begin{array}{c}{\partial }_x\\ {\partial }_y\end{array}\right)\delta (\mathbf{x}-{\mathbf{x}}_0).$$

(24)

For the numerical simulations, we use a representation of the Dirac delta:

$$\delta (\mathbf{x}-{\mathbf{x}}_0)\approx {\eta }_{\mathbf{\epsilon }}\left(\mathbf{x}\right)\equiv \frac{1}{2\pi \mathbf{\epsilon }}{e}^{-\frac{|\mathbf{x}-{\mathbf{x}}_0{|}^2}{2\mathbf{\epsilon }}},$$

so that (24) will be replaced by:

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}f\\ g\end{array}\right)& \approx & H\left(t\right)(\lambda I+2\mu \mathbf{n}{\mathbf{n}}^T)(\mathbf{x}-{\mathbf{x}}_0)\frac{1}{\mathbf{\epsilon }}{\eta }_{\mathbf{\epsilon }}\left(\mathbf{x}\right).\hfill \end{array}$$

(25)

The maximum of both components of this function is $\frac{c_p}{2\pi e^{1/2}{\mathbf{\epsilon }}^{3/2}}$, so we divide the above expression by this constant to obtain an approximation of the theoretical body force with reasonable magnitude.

To integrate numerically (23) and (24), we use a box-integration method with two staggered grids. We subdivide the interval $[0,l_x]$ in $n_x$ subintervals of equal length $h_x=l_x/n_x$, the interval $[0,l_y]$ in $n_y$ subintervals of equal length $h_y=l_y/n_y$, and introduce the primary grid points:

$$(x_i,y_j)=(i{h}_x,j{h}_y),i=0,1,\dots ,n_x,j=0,1,\dots ,n_y,$$

and the dual grid points:

$$(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})=(x_i+\frac{1}{2}{h}_x,y_j+\frac{1}{2}{h}_y),i=0,1,\dots ,n_x-1,j=0,1,\dots ,n_y-1.$$

The boxes for the primary grid are rectangles $I_{i,j}$ with length $h_x$ and height $h_y$, centered on the internal grid points $(x_i,y_j)$:

$$I_{i,j}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times [y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}],i=1,\dots ,n_x-1,j=1,\dots ,n_y-1,$$

while the boxes for the dual grid are rectangles $J_{i,j}$ with the same size, not touching the boundary, and centered on the internal grid points $(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}})$:

$$J_{i,j}=[x_{i-1},x_i]\times [y_{j-1},y_j],i=2,\dots ,n_x-1,j=2,\dots ,n_y-1.$$

We integrate the first equation in (23) over the primary boxes $I_{i,j}$ and the second equation over the dual boxes $J_{i,j}$, obtaining the equations:

$$\begin{array}{c}{\int \int }_{I_{i,j}}\rho u_{tt}\mathrm{d}x\mathrm{d}y={\int \int }_{I_{i,j}}{((\lambda +2\mu )u_x+\lambda v_y)}_x\mathrm{d}x\mathrm{d}y\hfill \end{array}$$

$$\begin{array}{c}+{\int \int }_{I_{i,j}}{(\mu v_x+\mu u_y)}_y\mathrm{d}x\mathrm{d}y+{\int \int }_{I_{i,j}}f\mathrm{d}x\mathrm{d}y,\hfill \\ {\int \int }_{J_{i,j}}\rho v_{tt}\mathrm{d}x\mathrm{d}y={\int \int }_{J_{i,j}}{(\mu v_x+\mu u_y)}_x\mathrm{d}x\mathrm{d}y\hfill \end{array}$$

(26)

$$\begin{array}{c}+{\int \int }_{J_{i,j}}{(\lambda u_x+(\lambda +2\mu )v_y)}_y\mathrm{d}x\mathrm{d}y+{\int \int }_{J_{i,j}}g\mathrm{d}x\mathrm{d}y.\hfill \end{array}$$

(27)

Introducing the notation

$$\begin{array}{c}{u}_{i,j}\left(t\right)=u(x_i,y_j,t),i=0,\dots ,n_x,j=0,\dots ,n_y,\hfill \\ v_{i,j}\left(t\right)=v(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t),i=1,\dots ,n_x,j=1,\dots ,n_y,\hfill \end{array}$$

we can approximate:

$$\begin{array}{cc}\hfill {\int \int }_{I_{i,j}}& \rho u_{tt}\mathrm{d}x\mathrm{d}y\approx h_x{h}_y\rho \frac{{\mathrm{d}}^2{u}_{i,j}}{\mathrm{d}{t}^2},{\int \int }_{I_{i,j}}f\mathrm{d}x\mathrm{d}y\approx h_x{h}_{y}f(x_i,y_j,t),\hfill \\ \hfill {\int \int }_{I_{i,j}}& {((\lambda +2\mu )u_x+\lambda v_y)}_x\mathrm{d}x\mathrm{d}y={\int }_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}{[(\lambda +2\mu )u_x+\lambda v_y]}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}\mathrm{d}y\hfill \\ \hfill \approx & \frac{h_y}{h_x}(\lambda +2\mu )(u_{i+1,j}-2u_{i,j}+u_{i-1,j})+\lambda (v_{i+1,j+1}-v_{i,j+1}-v_{i+1,j}+v_{i,j}),\hfill \\ \hfill {\int \int }_{I_{i,j}}& {(\mu v_x+\mu u_y)}_y\mathrm{d}x\mathrm{d}y={\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{[\mu v_x+\mu u_y]}_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\mathrm{d}x\hfill \\ \hfill \approx & \mu (v_{i+1,j+1}-v_{i,j+1}-v_{i+1,j}+v_{i,j})+\frac{h_x}{h_y}\mu (u_{i,j+1}-2u_{i,j}+u_{i,j-1}),\hfill \end{array}$$

which leads to

$$\begin{gathered} \begin{array}{c}\hfill \frac{{\mathrm{d}}^2{u}_{i,j}}{\mathrm{d}{t}^2}& =\frac{\lambda +2\mu }{h_x^2\rho }(u_{i-1,j}-2u_{i,j}+u_{i+1,j})+\frac{\mu }{h_y^2\rho }(u_{i,j-1}-2u_{i,j}+u_{i,j+1})\hfill & \\ +\frac{\lambda +\mu }{h_x{h}_y\rho }\left(v_{i+1,j+1}-v_{i,j+1}-v_{i+1,j}+v_{i,j}\right)+f\left(x_i,y_j,t\right).\hfill \end{array} \end{gathered}$$

(28)

In the same fashion, we can approximate:

$$\begin{array}{cc}\hfill {\int \int }_{J_{i,j}}& \rho v_{tt}\mathrm{d}x\mathrm{d}y\approx h_x{h}_y\rho \frac{{\mathrm{d}}^2{v}_{i,j}}{\mathrm{d}{t}^2},{\int \int }_{J_{i,j}}g\mathrm{d}x\mathrm{d}y\approx h_x{h}_{y}g(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t),\hfill \\ \hfill {\int \int }_{J_{i,j}}& {(\mu v_x+\mu u_y)}_x\mathrm{d}x\mathrm{d}y={\int }_{y_{j-1}}^{y_j}{[\mu v_x+\mu u_y]}_{x_{i-1}}^{x_i}\mathrm{d}y\hfill \\ \hfill \approx & \frac{h_y}{h_x}\mu (v_{i+1,j}-2v_{i,j}+v_{i-1,j})+\mu (u_{i,j}-u_{i,j-1}-u_{i-1,j}+u_{i-1,j-1}),\hfill \\ \hfill {\int \int }_{J_{i,j}}& {(\lambda u_x+(\lambda +2\mu )v_y)}_y\mathrm{d}x\mathrm{d}y={\int }_{x_{i-1}}^{x_i}{[\lambda u_x+(\lambda +2\mu )v_y]}_{y_{j-1}}^{y_j}\mathrm{d}x\hfill \\ \hfill \approx & \lambda (u_{i,j}-u_{i-1,j}-u_{i,j-1}+u_{i-1,j-1})+\frac{h_x}{h_y}(\lambda +2\mu )(v_{i,j+1}-2v_{i,j}+v_{i,j-1}),\hfill \end{array}$$

which leads to

$$\begin{array}{cc}\hfill \frac{{\mathrm{d}}^2{v}_{i,j}}{\mathrm{d}{t}^2}=& \frac{\mu }{h_x^2\rho }(v_{i-1,j}-2v_{i,j}+v_{i+1,j})+\frac{\lambda +2\mu }{h_y^2\rho }(v_{i,j-1}-2v_{i,j}+v_{i,j+1})\hfill \\ \hfill +& \frac{\lambda +\mu }{h_x{h}_y\rho }(u_{i-1,j-1}-u_{i-1,j}-u_{i,j-1}+u_{i,j})+g(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t).\hfill \end{array}$$

(29)

Equation (28) has been derived for $i=1,\dots ,n_x-1$, $j=1,\dots ,n_y-1$ and Equation (29) for $i=2,\dots ,n_x-1$, $j=2,\dots ,n_y-1$. We can extend them to the other values of i and j by using ghost nodes and enforcing the boundary conditions. We introduce the ghost nodes $x_{-1}=-h_x$, $x_{n_x+1}=l_x+h_x$, $y_{-1}=-h_y$, $y_{n_y+1}=l_y+h_y$, and the following mirror conditions with respect to the boundary, which enforce homogeneous Neumann conditions for u and v:

$$\begin{array}{c}{u}_{-1,j}=u_{1,j},u_{n_x+1,j}=u_{n_x-1,j},j=0,\dots ,n_y,\hfill \end{array}$$

(30)

$$\begin{array}{c}{u}_{i,-1}=u_{i,1},u_{i,n_y+1}=u_{i,n_y-1},i=0,\dots ,n_x,\hfill \end{array}$$

(31)

$$\begin{array}{c}{v}_{0,j}=u_{1,j},v_{n_x+1,j}=v_{n_x,j},j=0,\dots ,n_y+1,\hfill \end{array}$$

(32)

$$\begin{array}{c}{v}_{i,0}=u_{1,j},v_{i,n_y+1}=v_{i,n_y},i=0,\dots ,n_x+1.\hfill \end{array}$$

(33)

In this way, we obtain $(n_x+1)(n_y+1)$ Equation (28) for the unknowns $u_{i,j}$, $i=0,\dots ,n_x$, $j=0,\dots ,n_y$, coupled with $n_x{n}_y$ Equation (29) for the unknowns $v_{i,j}$, $i=1,\dots ,n_x$, $j=1,\dots ,n_y$.

After a standard finite difference approximation with three points for the second-order time derivative, with time step $\tau$, we obtain the following scheme:

$$\begin{array}{cc}\hfill u_{i,j}^{k+1}& =2u_{i,j}^k-u_{i,j}^{k-1}\hfill \\ \hfill & +\frac{(\lambda +2\mu ){\tau }^2}{h_x^2\rho }(u_{i-1,j}^k-2u_{i,j}^k+u_{i+1,j}^k)+\frac{\mu {\tau }^2}{h_y^2\rho }(u_{i,j-1}^k-2u_{i,j}^k+u_{i,j+1}^k)\hfill \\ \hfill & +\frac{(\lambda +\mu ){\tau }^2}{h_x{h}_y\rho }\left(v_{i+1,j+1}^k-v_{i,j+1}^k-v_{i+1,j}^k+v_{i,j}^k\right)+{\tau }^{2}f\left(x_i,y_j,t_k\right)\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill v_{i,j}^{k+1}& =2v_{i,j}^k-v_{i,j}^{k-1}\hfill \\ \hfill & +\frac{\mu {\tau }^2}{h_x^2\rho }\left(v_{i-1,j}^k-2v_{i,j}^k+v_{i+1,j}^k\right)+\frac{(\lambda +2\mu ){\tau }^2}{h_y^2\rho }\left(v_{i,j-1}^k-2v_{i,j}^k+v_{i,j+1}^k\right)\hfill \\ \hfill & +\frac{(\lambda +\mu ){\tau }^2}{h_x{h}_y\rho }\left(u_{i-1,j-1}^k-u_{i-1,j}^k-u_{i,j-1}^k+u_{i,j}^k\right)+{\tau }^{2}g\left(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}},t_k\right),\hfill \end{array}$$

(35)

where the superscript denotes the index of the time $t_k$. The above scheme provides the update for the unknowns on the grid points at time $t_{k+1}$ in terms of the same unknowns at times $t_{k-1}$, $t_k$. To initialize, we use:

$$\begin{array}{c}{u}_{i,j}^0=u_{i,j}\left(0\right)=0,v_{i,j}^0=v_{i,j}\left(0\right)=0\hfill \\ u_{i,j}^1\approx u_{i,j}\left(0\right)+\tau \frac{\mathrm{d}{u}_{i,j}}{\mathrm{d}t}\left(0\right)+\frac{{\tau }^2}{2}\frac{{\mathrm{d}}^2{u}_{i,j}}{\mathrm{d}{t}^2}\left(0\right)=\frac{{\tau }^2}{2}\frac{{\mathrm{d}}^2{u}_{i,j}}{\mathrm{d}{t}^2}\left(0\right),\hfill \\ v_{i,j}^1\approx v_{i,j}\left(0\right)+\tau \frac{\mathrm{d}{v}_{i,j}}{\mathrm{d}t}\left(0\right)+\frac{{\tau }^2}{2}\frac{{\mathrm{d}}^2{v}_{i,j}}{\mathrm{d}{t}^2}\left(0\right)=\frac{{\tau }^2}{2}\frac{{\mathrm{d}}^2{v}_{i,j}}{\mathrm{d}{t}^2}\left(0\right),\hfill \end{array}$$

where the time second derivatives at $t=0$ can be evaluated from the right-hand sides of (28) and (29) with $u_{i,j}=u_{i,j}^0$, $v_{i,j}=v_{i,j}^0$.

We have solved this system in a square with size 15 cm, with $250\times 250$ grid points, that is, space step $h_x=h_y=0.06$ cm and time step $\tau =0.08$ µs. We have used Lamé’s parameters $\lambda =9.7078$ GPa, $\mu =12.2883$ GPa, and density $\rho =2200$$\mathrm{kg}/{\mathrm{m}}^3$, which are appropriate values for mortar [28]. The numerical scheme has been implemented in Matlab, and the simulation time is approximately ten minutes. A crack formation has been simulated at the position $(x_0,y_0)=(5.3550,8.9550)$ between a couple of adjacent grid points, one on the primary grid and the other one on the dual grid, at a random time $t_0=3.75$ µs, with an angle $\theta =\pi /4$ between the crack line and the x-axis. The initial shape of the crack is shown in Figure 2.

For subsequent times, we show some snapshots of the solution in Figure 3 and Figure 4.

To validate the localization method discussed in Section 3, we have simulated four sensors, located in the positions $(x_1,y_1)=(4.92,0)$, $(x_2,y_2)=(0,9.9)$, $(x_3,y_3)=(9.9,14.94)$, and $(x_4,y_4)=(14.94,4.92)$. The values of u on each sensor, for all simulation times, are shown in Figure 5. The readings of the sensors should be interpreted as the arrival times of the elastic wave generated by the crack. To discriminate between false readings, we have applied the rule to take the first time for which $\left|u\right|$ is greater than 3% of its maximum value on the simulation time. Our preliminary results are quite encouraging. For the simulation performed in this section, the method gives an estimated location of the crack at (5.4438, 8.8606) and an estimated formation time of $11.7156$ µs. These values are remarkably close to the real crack location (5.355, 8.955) and formation time $t_0=12$, with a relative error of about $1\%$. It is worth noting that the linear initialization step (14) provides a poor result, estimating the crack at $(3.6848,10.0165)$ at (negative) time $-6.0384$ µs, while the nonlinear method reaches convergence with good results after three iteration steps.

5. Conclusions

In this paper, we have revisited the equations of elasticity, including a “body force term to be applied in the absence of a dislocation, which produces radiation identical to that of the dislocation”, citing the classical paper [27] where the general expression of this term was first derived. We have paid special attention to a specific kind of dislocation, “which represents a sudden separation of the two sides of the fault as, for example, an explosion in a plane crack” [27], because this mechanism seems adequate to describe what would happen in a concrete sample under pressure tests which produce microcrack formation. In Equation (7), we provide a generalization of the simple example in [27], which depends on the location of the crack and on its orientation. This is particularly relevant since, to our knowledge, no exact source function of acoustic emissions has been identified as of yet. For 2D numerical simulations, an approximated version of (7) has been provided in (25). The numerical results indicate that indeed it is a reasonable source.

The proposed method, through an experimental setup of acoustic emission acquisition, as shown in [17], could ensure the localization of cracks within concrete structures without the use of destructive or nondestructive survey methods. In fact, the latter are either difficult to implement, as in the case of X-ray microtomography, or are inaccurate, as in the case of ultrasonic methods. The data acquired in [17] can be processed by the method proposed in the article and, based on the acquired delay times, provide the location of the crack associated with the acoustic emissions.

This paper also aimed at establishing a theoretical foundation for the above localization method, based on the arrival times of acoustic emissions produced by a crack’s formation to a group of sensors dislocated on the boundary of a concrete sample. The preliminary results that we have obtained are quite encouraging and suggest the validity of the proposed localization method. Still, we believe that it is possible to improve the accuracy of the method. In our opinion, the main point is to improve the simulation of the reading of the sensors. Here, two observations are in order. The first one is that it is not clear how to determine to correct arrival time, due to the complexity of the detected signals, as shown in Figure 5. We have used the 3% rule, as stated at the end of last section, but other rules might need to be assessed. The second observation is that we made a working assumption that the arrival times do not depend on the angle $\theta$ formed by the crack line with the x-axis. Actually, there is numerical evidence that such dependence does exist. In other words, for deriving an improved localization method, it might be necessary to include the parameter $\theta$ as an additional value to be extracted from the sensors’ readings, together with the time of formation and the position of the crack. This extension of the localization method will be the subject of a subsequent paper.