Crack identification in hydro-mechanical systems with applications to gravity water dams

Abstract

This article investigates an inverse problem for hydro-mechanically coupled systems. The task is to identify damaged zones in water dams out of combined hydrologic and mechanical measurements. Even though water dams are generally well monitored, there is little data available to precisely describe the fluid and stress fields in the structure due to inhomogeneous distributions of the material. Over the lifetime of the dams, which can exceed 100 years, the condition of the structures may change remarkably due to aging, chemical weathering and cracking. Highly fissured regions may arise. Our goal is to identify the damaged zones by comparing the numerically computed fluid and stress fields with measurements along the boundaries. The identification is shown to be an ill-posed problem. For the inversion, regularizing iterative methods need to be applied, in which modified implementations of the Landweber method are used. The computation of the adjoint of the linearization of the parameter-to-solution mapping is done analytically. Due to the consideration of a multi-field problem, numerical results show that the identification results are significantly improved compared to the results of single-field problems.

Keywords:

• hydro-mechanical interaction
• flow in porous medium
• water dams
• inverse problem
• parameter identification
• regularization
• adjoint state method

AMS Subject Classifications::

• 86A22
• 35R30
• 65M32
• 35Q86
• 76S05

1. Introduction

Construction processes in civil engineering are increasingly supported by simulation models. The solution of many of these models consists of several physical quantities. For example, numerical models for water dams, embankments and subsoils include the mechanical, hydraulic and thermal fields. The precision of the output of these models depends on both the correct choice of the model and the correct adaptation of the model parameters. A careful model and material parameter adaptation to available measurements needs to be conducted. This article is devoted to an inverse problem for a hydro-mechanical problem, namely the identification of material inhomogeneities within the structures by evaluating combined hydraulic and mechanical measurements.

The system identification considered in this article can also be regarded as a method to monitor the embankment's condition. Undamaged systems have a homogeneous distribution of the material properties according to their construction process. However, damaged systems consist of one or several subdomains where the material properties differ remarkably within the structure due to partial cracking, weathering, elution or chemical reactions. Early recognition of such damaged zones is particularly essential in order to avoid seepage, piping and potential failure of the structure.

In the following, the identification of damage corresponds to the simultaneous determination of two material coefficients in a coupled parabolic system of partial differential equations. Related works treating coefficient determination problems for single-field parabolic problems can be found in 1–9 and references therein. For coupled-field parameter identification problems, we refer in particular to 10, where coupled groundwater flow and mass transport problems and general coupled problems are discussed.

The identification of damage, inhomogeneous stiffness properties and calibration of Finite Element models in dams is reported (e.g. 11–17). All these works, however, are restricted to measurements and observations of static or dynamic mechanical quantities.

To the author's knowledge, the consideration of both hydraulic and mechanical data to assess the dam's properties and the treatment of damage by heterogeneous parameter distributions is a new approach and has not yet been proposed in the literature.

The Landweber iteration, which is used in this work to solve the inverse problem of crack detection, was proven to converge for linear compact operators by Landweber 18. The convergence theory for the nonlinear problem has been developed by Binder et al. 19 and Hanke et al. 20. The application of Morozov's discrepancy principle 21 as a stopping rule is treated in 22,23. The nonlinear Landweber iteration with flexible steplength as used in this work is investigated in 24,25. For a survey of the Landweber iteration and its modifications, see 26. Landweber iteration is applied in many inverse problems occurring, for example, in image restoration 27 and inverse conductive scattering problems 28. The method has also been successfully applied to the recovery of parameter nonlinearities, e.g. heat transfer coefficient in 29 or nonlinear permittivity and piezoelectric coupling in 30 or nucleation and growth rates of polymer crystallization processes 31.

The article is organized as follows. Section 2 describes a simple case of fluid-solid interaction for small strain situations. A set of partial differential equations is derived. In Section 3, a numerical scheme for the inverse problem is presented, in which we apply an iterative regularizing scheme. The gradient is computed by solving a coupled adjoint state problem during each iteration step. In Section 4, numerical results are given.

2. Coupled flow and solid deformation

In this section, balance and constitutive equations are introduced which describe a simple form of fluid-flow through porous media. Out of these equations, a set of partial differential equations is then derived. In order to model cracks and damaged zones in the computational domain, a smeared crack model is considered, i.e. the most highly influencing material parameters are assumed to be continuously spatially varying quantities. For works that discuss discrete modelling of cracks in water dams, see, e.g. 32–34.

In the sequel, we consider a one-phase flow (water) through a deforming, fully saturated medium. Although we are aware of the fact that this scenario is not perfectly suitable for the simulation of water dams, we wish to keep the representation as simple as possible. Future work will be devoted to incorporating the empirical functions of van Genuchten 35, which describe the nonlinear relation between permeability and saturation, as well as saturation and suction. Since the task of this work is to develop a scheme which is capable of identifying material inhomogeneities, this simple approach serves the purpose sufficiently. Possible thermal effects are also neglected and the material is assumed to be isotropic.

The hydromechanical system is characterized by the deformation of the solid together with the flow of the fluid through the pores and the effects of pore water pressure distribution on the stress state of the structure. This kind of coupling is usually referred to as primary consolidation in soil mechanics 36.

Setting up the mathematical model of fluid flow-structure interactions involves the momentum conservation law, i.e. (1) where σ denotes the macroscopic total stress tensor (N m⁻²) and g is the acceleration due to gravity (m s⁻²). The averaged density ρ_a ≔ ((1 − n)ρ_s + nρ_w) (kg m⁻³) of the porous medium is determined from the densities of the solid ρ_s and the fluid ρ_w; n denotes the porosity. The operator ℬ is the strain-displacement differential operator Additionally, the mass balance of the liquid needs to be fulfilled, i.e. (2) Here , where K_s and K_w are the bulk moduli (N m⁻²) of the solid and water, respectively, and α denotes Biot's constant (dimensionless). Furthermore, u is the mechanical displacement vector (m), p_w liquid phase pore pressure (here water) (N m⁻²) and (see, e.g. 36). The term in (2) describes the effect of the volumetric strain rate on the pore pressure distribution.

The flow of the fluid is mainly characterized by the phenomenological law of Darcy. It states that the relative mass flux of the liquid is proportional to the negative gradient of the liquid pore pressure p_w, 36,37 (3) Here, q_rw (kg m⁻²s⁻¹) denotes the relative mass flux of the fluid. The isotropic intrinsic permeability tensor is given by (4) where k_rw is the relative liquid permeability (m²), μ_w is the liquid's viscosity (N m⁻² s⁻¹). On the mechanical side, the coupling between the hydraulic and mechanical quantities is described by the effective stress law which considers the additional contribution from the pore water pressure to the normal components of the mechanical stresses (5) In (5) σ′′ denotes the effective stress tensor.

The elastic behaviour of the structure is mainly described by Hooke's law where is the modulus of elasticity and δ_ij is Kronecker's delta.

The Lamé parameters (6) are themselves determined by Young's modulus E (N m⁻²) and Poisson's ratio ν.

The damaged zones in the dam are described by a smeared crack model, i.e. the main parameters, Young's modulus E in (6) and hydraulic conductivity k in (4) are spatially varying functions with H¹ regularity (one-time weakly differentiable) (7) Expressing the theory above solely in dependence of the unknowns, mechanical displacement u and liquid pore pressure p_w, we obtain the following set of partial differential equations: (8) with boundary conditions (9) and initial conditions (10) In (9) and N = (n_x, n_y, n_z) is the outer unit normal vector.

Solving the system in (8) with (9) and (10) will be referred to hereafter as the forward problem.

In this work a situation of plain strain is considered. Figure 1 shows the schematic shape of the computational domain Ω, which represents a profile of a gravity water dam. The upstream side is on the left-hand side. The boundary is divided into five subboundaries, i.e. ∂Ω = Γ₁ ∪ Γ₂ ∪ Γ₃ ∪ Γ₄ ∪ Γ₅ is indicated, too. The additionally sketched boundary Γ₆ consists of just one point and is regarded in Section 3 where an additional measurement is taken. The inhomogeneous Dirichlet boundary condition h in (9) simulates different filling levels of the reservoir.

Figure 1. Profile of gravity water dam. The location of the damaged zone is indicated with the boundaries labelled accordingly. The height of the modelled dam is about 29 m, the crest has a length of 166 m and the radius of curvature is about 120 m.

The numerical solution of the forward problem is done with the Finite Element Method where standard nodal finite elements are applied for the spatial discretization.

In order to derive the Finite Element formulation, the first equation in (8) is tested with vector valued functions having compact support, where and The second equation in (8) is tested with scalar functions. In order to obtain an algebraic equivalent of (8), Green's first integration theorem is applied to the resulting equations. The following FE approximation for the unknowns is then introduced (11) where N_u and N_p are elements of bases of (S_h)³ and T_h which are finite dimensional subspaces of and . The number of nodes in the FE-mesh, on which no homogeneous Dirichlet boundary conditions are defined, is denoted by n_eq.

The nodal displacements and nodal pore pressures are now obtained by solving or after symmetrization according to 36 The element matrices are computed as follows

	Qe = ∫Ωe (ℬNu)T αmNp dΩ (hydro-mechanical coupling matrix)
	Ke = ∫Ωe(ℬNu)T c(x)ℬNu dΩ (mechanical stiffness matrix)
	(permeability matrix)
	(capacity matrix).

For the discretization in time, the generalized trapezoidal method – an implicit time-stepping scheme – is applied 36,38. Figure 2 shows the pore pressure distribution and the horizontal and vertical component of the mechanical displacement at the final time t = T, where the water level lies between Γ₄ and Γ₅.

Figure 2. Solution of the forward problem: The first picture shows the distribution of the pore pressure for a fixed time t. The second and third pictures plot the mechanical displacement in the horizontal and vertical direction, respectively.

3. The inverse problem

This section deals with the inverse problem, i.e. the detection of the damaged zone within the dam structure out of the combined pore pressure and mechanical displacement measurements.

3.1. Measurements

It is assumed that over time the evolution of the mechanical displacement at the dam crest (two points) and on one point along the downstream side are recorded as follows (12) Additionally, it is assumed that there are five available hydraulic pore pressure measurements at the dam's foot (13) In the sequel, the measurements are stored in . Here and also in Equations (12) and (13), the superscribed δ emphasizes the fact that measurements are generally distorted by a certain amount of noise. The δ is introduced as L² measure of the noise level, i.e. Whenever identification is done with noiseless data, the superscribed δ is dropped: y = y^δ for δ = 0.

The actual discrete measurements given along the boundaries Γ₁ and Γ₃ are linearly interpolated. This is why measurements are considered not only at discrete points, but along the boundaries, as well. At Γ₆ a single measurement point is considered. The hydraulic quantities are recorded by pore water pressure transducers which are located at the dam's foundation. The values obtained from these transducers are generally used to determine the reduction of the hydraulic potential in the water dam (see, e.g. 39,40). The mechanical displacements at the crest of the dam and the downstream side are dictated either by plumb line or hanging pendulum (horizontal displacement), geometrical precision, differential levelling or geodetic positioning.

3.2. Mathematical formulation

The measurements ū^δ and are assumed to be functions in L²([0, T], L²(Γ_i)), i ∈ {u, p}. The most important instrument in inverse analysis is the so-called parameter-to-solution mapping (14) which relates the sought-for quantities to the measurements. The inverse problem with given measured data y^δ corresponds to solving (15) where p ≔ (E, k)^T.

During the inversion, i.e. the solution of (15) one computes values in H¹ out of noisy measurements in L². Due to the lack of regularity in the measured data and the fact that the measurements may not be attainable by F, the inversion needs to be regarded as an ill-posed problem. The lack of stability is a severe problem, as small error components in the measured data may lead to large identification errors.

3.3. Inversion

For the stable solution of the inverse problem, iterative regularizing methods are well proven, e.g. the nonlinear Landweber method 20 and its modifications 24–26, for n = 1, … , n_*(δ, y^δ) (16) Even though the measured data y^δ may not be in the range of F and the iterates p^n,δ do not converge to the exact solution p^†, the method in (16) allows for a stable approximation of p^†.

In order to implement the Landweber method efficiently, the following three questions must be answered:

1.	How do we choose the relaxation parameter ω in (16) to guarantee convergence which is sufficiently fast?
2.	How do we determine the stopping index n*(δ, yδ) which influences the trade-off between accuracy and stability?
3.	How do we compute the adjoint F′(p)*(r) of the linearized operator F′(p)[dp] applied on a residual r?
4.	Which rates of convergence can we expected?

3.3.1. Relaxation parameter

In order to guarantee convergence of the method in (16), the relaxation parameter ω needs to be chosen such that (17) for all p ∈ ℬ_2ρ(p⁰), where ℬ_2ρ(p⁰) denotes a closed ball of radius 2ρ around the initial guess p⁰. By guaranteeing that (17) is fulfilled, one implements the classical Landweber iteration 20,26.

Alternatively, one can apply the minimal error method which is a variant of Landweber's iteration with a flexible choice of the relaxation parameter ω ≔ ω^n,δ 24,25 with (18) The numerator and denominator in (18) are quantities which need to be computed anyway during the iteration (16). This can effectively improve the convergence of the Landweber method without increasing the computational efforts.

3.3.2. Regularization by a posteriori stopping

In order to avoid an amplification of the data error, the iteration in (16) is stopped by a generalized discrepancy principle, i.e. when for the first time (19) holds for τ > 2 (see 20,26). The last index n, for which (19) holds, is n_*(δ, y^δ). It is shown that an early stopping according to (19) renders the Landweber method a regularizing method 20,26,41, i.e. it guarantees that the solution is stably dependent on the measured data.

3.3.3. Computation of the adjoint problem

The adjoint F′(p)*(r) is given implicitly via the solution of a coupled hydro-mechanical set of differential equations.

Proposition 3.1

If we assume that the operator F defined in (14) is Fréchet differentiable, then the adjoint of the linearization of F is given by (20) where (ϕ, ψ) are solutions of the following adjoint system of partial differential equations (21) The operator Φ : 𝒟(Φ) = {ρ ∈ H²(Ω) | ∇ρ = 0 on ∂ Ω} → L²(Ω) in (20) is the adjoint of the embedding operator from H¹(Ω) into L²(Ω) where Φ is defined as (22)

Proof

The Fréchet derivative of F is given by the limit (23) where (du, dp_w)^T is the solution of the following system: (24) with the following boundary conditions: (25) and initial conditions (26) In (24) and (25) dc and dk are the stiffness and conductivity tensors computed with an incremented Young's modulus E + dE and hydraulic permeability k + dk, respectively. The adjoint of F′ is now given implicitly via the following identity: (27) In order to obtain an explicit representation of the adjoint, we multiply the two equations in (24) with a sufficiently smooth vector-valued function ϕ ∈ Ω × [0, T] and a scalar function ψ ∈ Ω × [0, T] respectively, integrate over space and time and add the resulting equations (28) Integrating with respect to space and time and rearranging the resulting equations in such a way that the terms du and dp_w can be factored out, gives (29) (30) (31) (32) (33) (34) (35) The terms in the rounded brackets in Equations (29) and (30) define the partial differential equations of the adjoint system given in (21). The entries in Equations (32) and (34) vanish due to the definition of the Neumann boundary conditions in (25). The contribution of the integral in (33) equals zero due to the prescribed end condition ψ(T) = 0 in (21) and the prescribed initial conditions dp_w(0) = 0 and du(0) = 0 in (26).

Writing down the remaining terms results in the following integral equality: (36) Let us rewrite the residual in (16) as follows: F(p^n,δ) − y^δ = (r_u, r_p)^T, where and , set N^T(cℬϕ) = r_u and along the boundaries Γ_u and Γ_p, respectively, and N^T(cℬϕ) = 0 and at the boundaries ∂Ω ∖ Γ_u and ∂Ω ∖ Γ_p.

Then (37) A comparison of (37) with (27) and the application of Φ⁻¹ implies (21).

The definition of Φ as in (22) is due to the fact that for any b ∈ H¹, the following holds: (38) ▪

Note, in particular, that the adjoint system has an end condition. It is, therefore, a well-posed initial-boundary value problem which is solved backwards in time. Introducing the following time transformation τ = T − t, the adjoint problem can be transferred into an initial-time problem with time variable τ. The solutions are then retrieved by ϕ(x, t) = ϕ(x, T − τ) and ψ(x, t) = ψ(x, T − τ).

If and in the case of noiseless data, where u and p_w solve (8) with E(x) and k(x), then F′(E, k)*(r_u, r_p) = 0 and E(x) and k(x) are minimizers of ‖F(E, k) − y‖².

For a detailed discussion on the computation of an adjoint problem for the identification of the diffusion coefficient in a linear (non-coupled) parabolic equation, see 8 and for coupled problems see 10. The Finite Element Method is also applied to solve the derived adjoint system which results in the algebraic system with right-hand sides according to the natural boundary conditions in (21).

We want to emphasize that, during each iteration step, all nodes of the parameterized distributions of E and k are updated by just one forward and adjoint problem call. This is especially superior to any gradient computation with finite difference approximations where the forward problem needs to be solved M + 1 times, when M is the number of nodes of the parameterization of E and k.

3.3.4. Convergence

In order to study the convergence of the Landweber method, one needs to assume that F is Frèchet-differentiable with ‖F′(·)‖ ≤ 1 and that (39) holds in 𝒟(F). Then, for any solution p ∈ 𝒟(F) of F(p) = y a sufficient condition for p^n+1,δ to be a better approximation of the exact solution denoted by p^† is that which also gives an estimate for the factor τ in (19). However, the condition in (39) most likely cannot be shown rigorously for the coupled identification problem here. Nonetheless, it is proven to hold for the parameter estimation problem in an uncoupled elliptic problem 26 and a harmonic problem in 30.

Convergence rates of Landweber iteration can only be shown if a source condition of the following type (40) is additionally fulfilled 20,26. The source condition implies a certain smoothness for the difference p^† − p⁰ for many examples if μ > 0. Under a further restriction of the nonlinearity of F, namely F′(p) = R_p F′(p^†), where R_p is a family of bounded linear operators with ‖R_p − I‖ ≤ C‖p − p^†‖, the rate of convergence can be shown to be (41) The stopping index n_*(δ, y^δ) is obtained by (19) and can be estimated in the following manner (42) For the special case in (40), the terms in (41) and (42) reduce to (43) which is the best possible rate of convergence to be expected. The condition in (40) means that one has to know all rough parts of the sought-for quantities up to a certain regularity in order to obtain the rate . However, the smoothness assumed on the difference p^† − p⁰ will hardly be fulfilled in the application example given here. Therefore, we expect a slower convergence during the numerical computations.

For the minimal error method convergence rates have only been proven in the case of exact data, i.e. for δ = 0. In this case, for the modified Landweber method with flexible relaxation parameter according to (18), one has (see, e.g. 25). (44)

4. Numerical results

This section discusses some numerical results obtained by the application of the inversion scheme. All results are obtained by evaluating synthetically generated data. These data are obtained by solving the forward problem with known parameter distributions. The solutions of the forward problem are then distorted by a one percent normally distributed noise term in order to model measurement errors. The generation of these synthetic measurements is performed on a fine mesh with 576 finite elements, whereas the reconstruction of the damaged zones is performed on a coarser mesh with 144 elements in order to avoid an inverse crime 42. An inverse crime is committed when almost the same computational models are used to both synthesize and invert data in an inverse problem. This must be avoided as the results may appear too optimistic.

The initial guess is a homogeneous distribution of the material parameters E and k. The exact solution, i.e. the location of the crack, is depicted in Figure 3(b). In the following, we applied two strategies. In the first, both parameters E and k have been updated independently, i.e. we solved the adjoint problem and obtained parameter increments for both, the permeability and Young's modulus in each iteration step in (16). Figure 4 shows the uncorrelated reconstruction results. The damaged location, i.e. the distribution of Young's modulus in Figure 4(a) and the distribution of the permeability in Figure 4(b), is identified acceptably well. However, as we see, the shape and location of the damage according to the distribution of E and k are not exactly the same. Also the reconstruction of Young's modulus could lead us to guess that there is another crack close to the dam crest, which is not the case according to the results of the permeability. In the second strategy, we consider the correlation of the sought-for quantities.

Figure 3. (a) Initial guess (homogeneous material distribution). (b) True solution, i.e. the location of the damage.

Figure 4. Reconstruction results where both E(x) (a) and k(x) (b) have been updated independently by the iteration in (16).

Now, the obtained results are combined immediately in the following manner: (45) in every step in (16), where the newly introduced parameter χ serves as an indicator of damage. The constants γ_E and γ_k compensate for the different orders of magnitude of E and k. On the left-hand side of Figure 5, the distribution of the parameter χ is shown, which gives a more precise reconstruction result due to the correlation of the parameters. The sharp boundaries around the damaged zone in the right-hand side of Figure 5(b) in are obtained by projecting the identified values of χ which exceed (or fall below) a certain value to one (or zero). This threshhold is found by indicating positions where the gradient of χ becomes large, i.e. where a remarkable change in the structure of the dam is observed. Numerical computations reveal that the successful identification of material inhomogeneities strongly depends on the location of the inhomogeneity and, in particular, on the distance between damage and measurements. Obviously, a crack close to the tip of the dam is mainly identified via the mechanical quantities (see Figure 6).

Figure 5. Reconstruction results. (a) Distribution of parameter χ(x) according to (45). (b) Sharpened solution.

Figure 6. Location and identification, i.e. the distribution of the parameter χ(x) according to (45), of a crack close to the dam's tip.

The single mechanical measurement along the boundary Γ₆ (see Figure 1) turns out to be indispensable for the identification of the cracked zone in both examples. Cracks close to the upstream side of the dam are observed to be less precisely identifiable than those on the downstream side. This is due to the lack of mechanical measurements under the water. However, these partially dissatisfying results can give valuable hints as to where different sensors can be properly placed along or within existing dam structures. As a consequence, the example may be further investigated to answer questions concerning optimal experimental design for hydromechanically coupled problems in the field of civil engineering.

Finally, results concerning convergence of the implemented scheme when data errors occur in the measurements are given in Figures 7 and 8. In Figure 7, the norm of the difference between the identified and the exact parameters over an increasing noise level is plotted. Figure 8 shows the decreasing number of iterations for increasing noise level (Figure 8(a)). The linear relation between the residual and noise level (Figure 8(b)) is given by the stopping rule in (19). The numerical convergence results are consistent with the theory given in Equation (43).

Figure 7. (a) Error e_E(δ) ≔ ‖E^† − E^{n_*(δ,yδ),δ}‖ over increasing data noise δ. (b) Error e_k(δ) ≔ ‖k^† − k^{n_*(δ,yδ),δ}‖ over increasing data noise δ.

Figure 8. (a) Number of iterations n_*(δ, y^δ) for different noise levels δ. (b) Residual res(δ) = ‖F(p^n,δ) − y^δ‖ over δ.

Aside from the noise level and the distance of the damage to the measurements, the size of the damaged zone strongly influences the quality of the identification results. Obviously, a small crack only slightly influences the fluid and mechanical fields and therefore, there is little to measure that could indicate damage. Figure 9 displays the relation between the size of the damaged area and the quality of the identification.

Figure 9. Relative reconstruction error over size of the damaged zone.

5. Summary and outlook

This work presents a numerical scheme for the identification of irregular material distributions in coupled hydro-mechanical systems. Furthermore, it demonstrates that more information is gained by evaluating data from the two physical fields involved. Future research will be devoted to the description of the fluid-flow in the dam by considering the effects of partially saturated media which incorporate the nonlinear relationships between saturation and suction such as permeability and suction (see 35). Additionally, the procedure will be applied to measurements of arch type gravity dams run by the Ruhrverband Essen in Germany.

As mentioned above, the results of this work can be further explored for an optimal experimental design, offering clues as to where measurement devices should be optimally located and how many of them are required in order to detect damage in all parts of the structure reliably.

Acknowledgements

The author gratefully acknowledges the support of the German Research Foundation (DFG) through the Research Training Group 1462. Furthermore, the author wishes to kindly thank Professor C. Könke and Professor T. Schanz for their motivation and initialization of this research.