On adaptive refinement techniques in multi-field

Comput. Methods Appl. Mech. Engrg. 195 (2006) 444–461 www.elsevier.com/locate/cma

On adaptive refinement techniques in multi-field problems including cohesive fracture Bernhard A. Schrefler a

a,*

, Stefano Secchi b, Luciano Simoni

a

Department of Structural and Transportation Engineering, University of Padova, Via Marzolo 9, 35121 Padova, Italy b CNR, Isib, Corso Stati Uniti 4, 35127 Padova, Italy Received 8 March 2004; received in revised form 12 October 2004; accepted 29 October 2004

Abstract This paper presents a generalized finite element formulation which incorporates solid and fluid phases together with a temperature field. The model is developed to obtain time-dependent solutions of complex 2-D cases, such as concrete gravity dams subjected to loading–unloading cycles, non-homogeneous specimens subjected to thermo-mechanical effects, etc. The solid behaviour incorporates a fully coupled cohesive-fracture discrete model, which includes thermal and hydraulic loads and the resulting crack nucleation and propagation is fully described. The evolution of fractures leads to continuous topological changes of the domain and these are handled by systematic local remeshing of the domain and a corresponding change of fluid and thermal boundary conditions. Optimality of the size of automatically generated finite elements is controlled, and the mesh density is adaptively adjusted on the basis of an a posteriori error estimation. For the process zone an element threshold number is introduced to obtain mesh independent results. The presented applications demonstrate the efficiency of the procedure and the importance of mesh refinement in multi-physics problems. 2005 Elsevier B.V. All rights reserved. Keywords: Porous media; Cohesive fracture; Non-homogeneous material; Thermal effects; Hydraulic fracture; Delaunay tessellation

1. Introduction The mechanical behaviour of porous solid materials, e.g. concrete, soils, biological tissues, etc., naturally falls into the category of multi-physics problems, which are usually described by use of a macroscopic approach, as for instance the mixture or the averaging theories [1]. Particular care must be taken in the *

Corresponding author. Tel.: +39 049 827 5611; fax: +39 049 827 5604. E-mail addresses: [email protected] (B.A. Schrefler), [email protected] (S. Secchi), [email protected] (L. Simoni). URL: http://www.dic.unipd.it (L. Simoni).

0045-7825/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.10.014

B.A. Schrefler et al. / Comput. Methods Appl. Mech. Engrg. 195 (2006) 444–461

445

definition of the coupling terms, which describe the effects of the other fields on the analysed one. From the computational point of view, further difficulties arise in the approximation of the field variables, which in general possess completely different behaviour in space, especially close to singularities, and in time. We refer, as an example, to the formation of the fluid lag near the tip of a crack moving in a saturated porous medium. In a linear elastic fracture mechanics context it has been demonstrated that a finite lag is required to ensure coherence of the mathematical solution [2]. A combination of the water mass balance equation and the singular solution for the solid, which is r0.5 for the stress field, implies a logarithmic singularity for the pressure. Under some simplifying assumptions, however, pressure near the tip exhibits an r0.5 behaviour [3] which is similar to some experimental results. The tip cavity is in fact filled with evaporated fracturing fluid under constant pressure, negligibly small when compared with the stress field far from the fracture, which justifies the assumption of zero pressure value at the tip [4]. Whereas the crack tip moves at a speed depending on the mechanical characteristics of the medium, the fluid within the fracture moves depending on permeability, for instance defined by the Poiseuille law [5]. The lag dimensions hence change in time. Similar conclusions can be drawn in the presence of cohesive forces for the solid. The importance of properly modelling the fluid field stems from the fact that pore pressure acts as a body force on the solid skeleton. On the other hand, when the crack takes place in a thermo-elastic medium and is driven by the thermal field, this is practically not influenced by the fracture, while the solid exhibits the usual square root singularity, if linear elastic behaviour is assumed. This is not the case in the presence of cohesive forces, owing that plasticity smoothes the singularity of the stress field. In such cases, when modelling the involved fields, a unique singular approximation can not be used, as in the case of the quarter point elements in linear elastic fracture mechanics. But this usual approximation can also be questioned in the presence of nonhomogeneous solid, owing that the strength of the singularity depends on the mechanical properties of the different solids which form the medium. In similar situations the adaptive technique plays a role of paramount importance to obtain a correct numerical solution. On the base of these and similar observations, we can conclude that the solid field possesses a clear singularity, yet not a definite time scale: in computational fracture mechanics time is usually disregarded, even though time dependence of fracture has been introduced, e.g. by Carpinteri et al. [6]. On the other hand, the temperature and pressure fields do not exhibit singularities, but possess a time scale depending on the specific transport laws. Coupling terms transfer the time scale to the solid. Moreover, assuming the fracture path as unknown, the use of refinement techniques is mandatory to obtain the solution of this kind of problems. Refinement strategies in space have been widely investigated in recent years, see for instance Zienkiewicz and Taylor [7], whereas refinement in time has obtained less attention. A simple and attractive adaptive time-stepping procedure has been proposed by Wiberg and Li [8] limited to linear dynamics problems and could be adapted for first order time differential equations. These strategies could also be applied to multi-physics problems, but the error estimation, which is the fundamental part for refinement, should be revisited and coupling terms should be involved. An adaptive method has been recently presented for thermo-mechanical coupled problems with contact by Rieger and Wriggers [9] by using a time-space-finite element discretization in conjunction with a staggered solution algorithm. In the following of the paper we do not use this approach, due to the fact that coupling between pressure and displacement fields is stronger that between temperature displacement fields. This introduces further numerical difficulties, also shown by the fact that partitioned solutions without iterations produce wrong results. Usual refinement procedures in space will be used, whereas in time several simulations with different time steps will be performed in order to estimate the local error. The discussion of a benchmark test case will put in evidence the capabilities and the limits of the refinement techniques which are usually adopted in the space domain and the needs for further developments will be evidenced. This paper presents a generalized finite element formulation which incorporates solid and fluid phases together with a temperature field [1]. The model is developed to obtain time-dependent solutions of 2-D

446


cases, such as concrete dams subjected to loading–unloading cycles, non-homogeneous specimens subjected to thermo-mechanical effects, etc. A fully coupled cohesive-fracture discrete model, which includes thermal and hydraulic loads, is adopted to analyse crack nucleation and propagation. The evolution of fractures leads to continuous topological changes of the domain. No special approximations are used for the field variables to represent the possible singularities: bilinear interpolations are assumed for all of them in conjunction with h-adaptive algorithms based on a posteriori local error estimations [10]. Spatial discretization is continuously updated as the phenomenon evolves and the domain of definition changes together with boundaries and related conditions. This technique allows to define the initial mesh without special care and it is also able to handle non-homogeneous materials. Optimality of the size of automatically generated finite elements is controlled, and the mesh density is adaptively adjusted avoiding, at the same time, the presence of too distorted elements [11]. The new element size distribution allows for grid enrichment and derefinement and is particularly useful in the presence of evolutionary problems, e.g. in transport phenomena with phase changes located on moving fronts, strain localization and plasticity problems, solid–solid interfaces etc. Cracks may nucleate everywhere depending only on the stress field and propagate along paths and with a velocity of the tip that is unknown a priori. Also the fluid velocity inside the crack is a priori unknown, as the possibility of fluid lag formation. The determination of the crack path, the velocity of the tip propagation and the velocity of the fluid represent an important part of the solution, as the pressure, temperature and stress fields and allows for correct updating of the domain [12]. These difficulties are peculiar for multiphysics problems and refinement techniques must be adapted to such situations. One of the most important aspects considered in the paper deals with the discretization of the process zone, i.e. the area where the cohesive tractions develop. In particular the minimum number of elements located in the process zone is detected (element threshold number) in order to obtain a mesh independent solution. A goal-oriented adaptivity is used, i.e. the adaptivity problem is restricted to the process zone whereas the quantity of interest is assumed far away (for instance an applied load). This allows to construct a refinement for the control of the solution in terms of cohesive forces and trial advancing step of the crack Ds. It will be shown that also time discretization must be considered in order to obtain mesh independent results. 2. The field equations Governing equations of the analysed problems are briefly recalled. For a more detailed presentation see Lewis and Schrefler [1] and Simoni and Secchi [12]. Within the framework of BiotÕs theory, non-isothermal, quasi-static conditions, small displacements and displacement gradients are assumed. The mechanical behaviour of the solid is dependent on the effective stress r0ij defined, following Biot and Willis as a ð1Þ r0ij ¼ cijrs ðers drs eT Þ apdij ; eT ¼ T ; 3 ers being the total strain tensor, p the fluid pressure, dij the Kronecker symbol, a BiotÕs coefficient, which accounts for small volumetric strain due to pressure. eT is the strain associated to temperature T changes, according to cubic expansion coefficient a. As usual, a Green-elastic material is assumed for the solid at a distance from the process zone, being cijrs the elastic coefficients dependent on the strain energy function W. The linear momentum balance for the mixture (solid plus water), in weak form, hence containing the natural boundary conditions, may be written as: Z Z Z Z Z Z a deij cijrs eij dX qdui gi dX deij dui ti dC dui ci dC0 ¼ 0; adij p dX deij cijrs dij T dX 3 X X X X Ce C0 ð2Þ


447

where X is the domain of the initial boundary value problem, Ce is the external boundary and C 0 the boundary of the fracture and process zone (Fig. 1). deij is the strain associated with virtual displacement dui,q the density of the mixture, gi the gravity acceleration vector, ti the traction on boundary Ce and ci the cohesive tractions on the process zone as defined in the following. Forced conditions fix the field variables along the constrained boundary and completely define the problem. The fracturing material in the process zone undergoes mixed mode crack opening. This is usually the case of a crack moving along an interface separating two solid components. Whereas the crack path in a homogeneous medium may be governed by the principal stress direction, the interface has an orientation that is generally different. The mixed cohesive mechanical model involves the simultaneous activation of normal and tangential displacement discontinuity with respect to the crack and corresponding tractions. The interaction between the two cohesive mechanisms is treated as in Margolin [13] and Dienes [14], which generalizes the Hilleborg [15] model for mode I cohesive behaviour. Before presenting the mass balance equation for water, we have to discuss with a certain care the definition of permeability tensor. DarcyÕs law with constant absolute permeability is assumed for the fluid fully saturated medium surrounding the fracture. Within the crack the Poiseuille or cubic law is assumed to be valid. This has been stated for a long time in the case of laminar flow through open fractures and its validity has been confirmed in the case of closed fractures where the surfaces are in contact and the aperture decreases under applied stress, the sides remaining parallel. Permeability is not dependent on the rock type or stress history, but is defined by crack aperture only in the form [5] k ij ¼

1 w2 ; f 12

ð3Þ

w being the fracture aperture. f a coefficient is in the range 1.04–1.65 and takes into account the deviation from the ideal parallel surfaces conditions. In the following, this parameter is assumed as constant and equal to 1.0. Incorporating DarcyÕs law, the weak form of the mass balance equation for water in all the domain except for the fracture zone is given by Z Z n op oT k ij an s þ avi;i ½ð dp þ a nÞa þ nb p;j þ qw gj dX dX ðdpÞ;i K w ot ot lw Ks X X Z Z dpqw dC þ dp qw dC0 ¼ 0; ð4Þ þ Ce

C0

Fig. 1. Geometry of the problem.

448


where X is the domain of the fluid and solid fields (Fig. 1), dp is a continuous pressure distribution satisfying boundary conditions, n the porosity, Ks and Kw the bulk modulus for solid and liquid phases, b the thermal expansion coefficient of water, vsi the velocity vector of the solid phase, kij the permeability tensor of the medium, lw the dynamic viscosity of water, qw its density and qw the imposed flux on the external boundary. In the last term of Eq. (4) qw represents the water leakage flux along the fracture toward the surrounding medium. This term is defined along the entire fracture, i.e. the open part and the process zone. Incorporating fluid transport law into the water mass balance equation within the crack results in Z

2 Z Z n op ow oT w þ nb dp p;j þ qw gj dX þ dpqw dC0 ¼ 0; dX ðdpÞ;i K w ot ot ot 12lw X X C0

ð5Þ

which represents the fluid flow equation along the fracture. In this equation, X and C 0 represent the domain and the boundary of the fracture (Fig. 1). It should be noted that the last term, representing the leakage flux into the surrounding porous medium across the fracture borders, is of paramount importance in hydraulic fracturing techniques. In the present formulation, this term does not require special assumptions, as in Boone and Ingraffea [16], for instance. It can be represented by means of DarcyÕs law using the permeability of the surrounding medium and pressure gradient generated by the application of water pressure on the fracture lips. Except for the compressibility term and thermal effects, this equation is also presented in Ref. [16], even though it is there treated in a completely different way from the numerical point of view. Heat transfer is addressed next. When mechanical terms are neglected, internal energy depends on temperature only and is related to heat capacity of the mixture at constant volume Cv. Volume heat sources are however retained, owing that they can represent very different phenomena, e.g. heat production due to hydration of concrete or to plastic work. Hence they represent coupling effects between stress and thermal fields. Source terms may also arise along the boundary and represent frictional effects. In particular, the rate at which heat is generated at the contact surfaces is r0 ¼ t kvk;

ð6Þ

where t is the contact traction and kvk is the jump in velocity across the contact. This effect can be incorporated in the weak form of the energy balance, which hence takes the form [18] Z Z Z Z Z Z 0 0 conv _ qðC v T Þh dX þ r h dC þ q h dC þ qh dC ¼ q div h dX þ sh dX ð7Þ X

C0

C0

C0

X

X

being h an admissible virtual temperature, qconv and q the convective and imposed heat flux normal to the boundary. FourierÕs law is used as constitutive assumption for heat flux (kij being the effective thermal conductivity tensor), and NewtonÕs law to represent convective flux (being h the convective heat transfer coefficient and T1 the temperature in the far field of the undisturbed surroundings), qi ¼ kij T ;j

qconv ¼ hðT T 1 Þ. i

ð8Þ

Except for the convective and boundary fluxes, this equation is similar to the one used by Camacho and Ortiz [17]. The governing equations (2), (4), (5) and (7) are firstly discretized by means of the Galerkin procedure, then solved simultaneously to obtain the displacement and pressure and temperature fields together with the fracture path. The topology of the domain X and boundary change with the evolution of the fracture phenomenon. In particular, the fracture path, the position of the process zone and the cohesive forces are unknown and must be determined during the analysis.


449

3. Generation of the mesh In the following applications an improved implementation of 2-D Delaunay triangulation is used leading to unstructured planar meshes defined in simply or multiply connected regions of any arbitrary shape [11]. The algorithms belonging to the Delaunay–Voronoi class offer the advantage of efficiency together with a sound mathematical basis. The computational cost of the used algorithm is ranging from O(N1.5) to O(N log N), but, avoiding some not strictly necessary quality checks, can be reduced to almost O(N). Further, it pursues flexibility, reduction of user intervention, quality of the mesh and computational efficiency, independently of the complexity of the domain. To obtain high computational performances, the capabilities of an object-oriented language are exploited and a modified edge-based data structure is built for storing mesh topology information and handling meshing operations. A series of newly defined primitives and higher-level topological operators perform the necessary operations, speed up the procedure and limit memory requirements. The defined primitives and operators, together with an edge-based data structure almost completely avoid the need of computational loops in meshing operations, resulting in a highly efficient generation procedure. The only user intervention consists in a schematic description of the boundary of the domain (polygonal contours and/or analytical curves) and the location of singularities, for instance sources or sinks. Each geometrical entity is associated with a function controlling the spacing of the nodes to be generated (spacing function). Automatic refinements deal with a priori node grading along the boundaries and with the management of mesh and geometric singularities. Once new nodes are generated, the spacing function is interpolated from the existing nodes. The final discretization depends only on the initial description of the boundary (assigned nodes and spacing function) and there is absolutely no guarantee that with arbitrary input data a good quality mesh can be obtained. To improve this situation, a special algorithm is defined which corrects the boundary description when highly distorted elements originate from the assigned data, running from the initial steps of the generation. Particularly interesting for non homogeneous materials and fracture problems are some capabilities of the generator such as the node grading along the boundaries, management of mesh and geometric singularities and a strict preservation of the shape of the singularity and/or the interface. In the case of non-homogeneous materials, each homogeneous region is meshed independently from the others, using completely independent spacing function values with the unique requirement of assuming the same value along the interfaces for the neighbouring zones in order to guarantee mesh consistency and approximation continuity. Further, it is necessary that the triangulation completely represents the boundary of each region and the complete domain. It should be remembered that Delaunay triangulation is a unique geometric construction for an assigned set of points, hence it does not guarantee the consistency of a generic contour defined through a subset of points of the initial set. The interested reader is referred to Secchi and Simoni [11] for more details. Evolution of discrete fracture results in a continuous change in the topology of the domain together with the updating of boundary conditions. Whereas changes in topology are dealt with by the remeshing of the domain, mechanical and kinematical boundary conditions are treated by means of the Lagrange multipliers technique. To model the open crack and process zone, finite elements are used that for the solid phase present only the cohesive forces and account for the complete behaviour for the fluid phase. As a consequence, these elements allow for a straightforward integration of Eq. (5) in the fracture domain, without requiring special treatment of the continuity equation along the open fractures. Finite elements located along the fracture path can be viewed as special non-linear contact elements, capable of representing mixed mode fracture behaviour, as well as integrating the fluid continuity equation along the crack path. More details in the representation of crack nucleation and propagation can be found in Ref. [12]. When geometric or mesh singularities are located near the boundary, generating high gradients of the spacing requirements (Fig. 2), the algorithm updates the grading along the boundary and improves the geometric description, maintaining the topology of the domain (solution independent mesh improvements). The

450


Fig. 2. Mesh with (bottom) and without (top) solution independent multi-constrained refinement.

same procedure which controls the grading along the boundary allows for quality control of the complete resulting grid. In this sense the subdivisions produced, although not optimized as finite element meshes, are certainly very appropriate from the geometric point of view and represent an optimal starting point for adaptive algorithms in numerical solutions. Good meshes are hence obtained by means of a unique process, without any correction phase, even in the presence of very irregular boundaries and geometric constraints. Once the triangulation has been obtained, it can be converted into a quadrilateral mesh, which is sometimes preferred in finite element applications. The quality of the transformed mesh is ensured from the beginning on, by optimising the aspect ratio for the generated quadrilaterals and by performing logical checks to choose between the different possibilities offered by the conversion. 4. Refinement strategy In the frame of h-refinement method in space and in spite of its cost, a total or partial remeshing technique is adopted in the following. For all variables, linear approximation of the field is used. To reduce computation time, the successive refinement/derefinement operations can be limited to suitable sub-domains containing the singularities zones (Fig. 1). The dimensions of these areas and local density of discretization can be decided on the basis of the solution at the previous steps. The main reason for adopting this type of remeshing is to preserve a good mesh, which is necessary especially for modelling the transport mechanisms, even though some problems may arise for transferring data from a mesh to another. The link between the refinement and the mesh generator is the spacing function, which is point-wise updated according to the a posteriori calculated error. The spacing function and the ensuing node distribution is then regularized by interpolation when the new nodes are inserted. In this way the resulting mesh always present a smooth distribution of element dimensions. The general scheme of the refinement procedure is presented in Fig. 3. It is important to remark that the first step is performed by the mesh generator using the multi-constrained algorithm [11], when a geometrical singularity is detected. Then the Zhu–Zienkiewicz [10] technique is adopted, which relies on an a posteriori recovery-based error estimator. The error energy norm is calculated locally over a patch constituted by the elements, usually six, surrounding each node of the actual mesh. This error is related to the maximum permissible error. A weighting parameter is hence obtained by which the spacing function in the central node of each patch is multiplied. Note that in the used mesh generator the spacing function is defined point-wise, whereas in standard adaptive procedures the generic element is directly handled. The used approach has the


451

Fig. 3. General refinement procedure.

advantage of producing regular and graded meshes. All field variables are involved in the adaptive procedure, however, in our applications, the stricter requirements come from the solid field. Particular care has to be used in handling the discretization of the cohesive zone and the crack tip movement in order to avoid mesh dependency of the solution. Zhou and Molinari [19] suggest three basic rules: (i) small mesh size, (ii) mesh size as uniform as possible, (iii) random orientation of elements in order to avoid preferred directions. The adopted strategy guarantees all these requirements. In fact, in correspondence of the tip of the fracture we locate a point element source with strength controlled by the user or by the discretization error; uniformity of the mesh in the area surrounding the crack is assured by the interpolation of the spacing function and by the remeshing; finally, the random distribution of elements is a characteristic of Delaunay tessellation, without using highly distorted elements. Even though the above presented mesh adaptivity algorithm is always active, for relative energy norm percentage errors typical of engineering we have found that the most pressing requirements are those of a proper representation of cohesive tractions. The same requirements allows for a good representation of the fluid field in the process zone and the fluid lag. The refinement procedure of the process zone is as follows: the length of the process zone is a priori estimated for the assigned material properties and a geometry similar to the case at hand (in the application BarenblattÕs theory [20] has been used). Once this length is obtained, the analyst can choose the number of elements to discretize it (element threshold number). This allows to construct a goal-oriented refinement (Fig. 4) and to control the trial advancing step Ds. As far as adaptivity in time is concerned, the importance of refinement will be assessed by repeating the complete transient analysis using reduced time steps. This is due to a lack of a suitable error estimator for non-linear coupled problems. However, this rough choice allows the assessment of the effects of space/time discretizations. 5. Numerical applications Three numerical applications are presented in order to demonstrate capabilities and deficiencies of the existing adaptive techniques when applied to multi-physics problems. For a better understanding of the results, solid-fluid and solid-temperature coupling are separately addressed.

452


Fig. 4. Goal-oriented refinement procedure for the process zone.

The first application deals with an hydraulically driven fracture due to a fluid pumped at constant flow rate. Assuming 2-D conditions (plane strain), simplified analytical solutions have been obtained and a more general numerical model has been presented by Boone and Ingraffea [16]. Fig. 5 presents the geometry of the problem together with the finite element discretization. A notch with a sharp tip is present along the symmetry axis of the analysed area. Material properties are presented in Table 1. This test is not presented for comparison purposes, but only to discuss the problems of refinement. The main characteristic of this test case is the possibility to change the time scale of the coupled problem by fixing different water inflow rates at the crack mouth. The effects of combined spatial/temporal discretizations is clearly seen in Fig. 6, where the crack length is presented versus time for different tip advancement, Ds, and time step increments, Dt. We can conclude that the more correct time history (case E) is obtained by simultaneously reducing these two parameters, whereas the reduction of only one discretization parameters leads to errors (about ±20%) even using a small tip advancement, if compared to the crack length. At the same time, the importance of the element threshold number is evident for the choice of Ds (crack length, according to BarenblattÕs hypotheses is about 0.9 m). We can also conclude that crack tip velocity is very mesh sensitive. Hence the element number threshold must be satisfied to obtain mesh independent results. A lower number of elements result in wrong crack tip velocity and the development of fluid lag may be missed. This in fact is the result of the interplay between crack velocity and permeability (through continuity equation) and is of great importance because it determines if there is negative pressure in the process zone, hence determines different body forces. As can be seen from Fig. 7, where time histories of the crack mouth opening and fluid pressure are presented, not all parameters are sensitive to crack tip advancement. For instance pressure at the mouth is nearly independent of Ds (Fig. 7b) and of Dt (not shown in the figure due to the overlapping of the responses). We only observe a more regular curve for the lower values of Ds. The crack mouth opening


453

Fig. 5. Problem geometry for water injection benchmark and overall discretization.

Table 1 Material properties for water injected test case Permeability coefficient Shear modulus PoissonÕs ratio Bulk modulus of solid Bulk modulus of fluid Porosity Fluid viscosity Max allowable stress (mode I) Critical width

2 · 105 m2/(MPa s) 6000 MPa 0.2 36,000 MPa 3000 MPa 0.19 109 MPa s 1.1 MPa 0.26 mm

exhibits an intermediate behaviour: its time history slightly depends on Ds, whereas is independent of Dt (results not shown in the figure, owing that the response curves are practically overlapping). The second application deals with the benchmark exercise A2 proposed by ICOLD [21]. The benchmark consists in the evaluation of failure conditions as a consequence of an overtopping wave acting on a concrete gravity dam. The geometry of the dam is shown in Fig. 8 together with initial, boundary conditions and an intermediate discretization. Differently from the original benchmark, the dam concrete foundation is also considered, which has been assumed homogeneous with the dam body. In such a situation, the crack path is unknown. On the contrary, when a rock foundation is present, the crack naturally develops at the interface between dam and foundation.

454


Fig. 6. Crack length time history.

Fig. 7. Crack mouth opening displacement and mouth pressure time histories.

As far as initial conditions for water pressure are concerned, it is assumed that during building operations and before filling up the reservoir, pressure can dissipate in all the dam body. As consequence zero initial pore pressure is assumed in the simulation. A more realistic assumption is the hypothesis of partial saturation of the concrete, which would require a further extension of the present mathematical model. Applied loads are the dam self-weight and the hydrostatic pressure due to water in the reservoir growing from zero to the overtopping level h (which is higher than the dam). The material data for the concrete are those assigned in the benchmark, whereas for permeability the value of 1012 cm/s has been assumed. This value could suggest the hypothesis of an impermeable material. This limit case can be analysed by the present model locating the diffusion phenomenon in restricted areas near the wetted side of the dam and along the crack sides. Such a condition is easily handled by the used mesh generator, but has not been applied in the following. The necessity of adaptivity during the solution is demonstrated in Fig. 8 where two different crack paths are obtained depending on the assumed permeability of the fluid within the fracture. Also the proper representation of the cohesive forces requires a fine mesh in the area of the process zone. This is shown in Fig. 9, which represents the process zone when the fracture length is 3.5 m and corresponds to an intermediate step of the analysis when the water level is 78 m. Finally, the formation of the fluid lag is studied. The lag is dependent on the different velocities of propagation of the crack tip and the one of


455

Fig. 8. Problem geometry for ICOLD benchmark and calculated crack position.

Fig. 9. Zoom near the fracture for maximum principal stress contour.

seepage inside the fracture, hence the simulation of this feature requires a simultaneous correct representation of the solid and fluid field. Zhu–ZienkiewiczÕs [10] adaptive strategy for gradient dependent quantities (seepage velocity) and goal-oriented refinement for crack velocity are in this case very useful. Fig. 10 depicts the fluid lag (water pressure is compression positive). Some important conclusions can be drawn from this application: • The mechanical behaviour of the solid skeleton strongly depends on the characteristic permeability of the fluid within the crack. Crack paths are in fact different, as result of the different stress fields.

456


Fig. 10. Zoom for pressure distribution within the crack and fluid lag.

• The fluid lag is responsible of the differences in the stress field in the process zone, hence correct modelling of the fluid and solid fields is mandatory. • Crack path can not be forecast, hence the traditional use of special/interface elements to simulate fracture propagation in large structures is prevented. The alternative to the successive remeshing is the use of cumbersome discretizations of the areas interested by fracture, but also this strategy is not viable in the case of dams. Further, the used technique for the analysis of the nucleation of the fracture does not require the presence of an initial notch and requires a very limited amount of information to be initially defined. The third application is related to temperature driven fracture propagation. Reference is made to the laboratory experiment [22] for a three-point bending test performed on a bi-material specimen subjected to a thermo-mechanical loading. The geometry is presented in Fig. 11. The sample is made by bonding with methacrylate adhesive a part of aluminium 6061 and another of polymethylmethacrylate (PMMA) and presents a notch with a sharp tip of 1 mm width and 30 mm height shifted 3 mm from the interface in the PMMA zone. The two materials present very different properties (Table 2), so that, when the system is subjected to heat, stresses arise near the interface as a result of the mismatch in thermal expansion. A variation of PMMA YoungÕs modulus with temperature has been experimentally observed and values are shown in Table 3 were measured. Fracture energy has been obtained on the basis of the stress intensity factors declared [22] for the two limiting cases (T = 25 C, 60 C), whereas the ultimate stress r0 is assumed of 25 MPa at the temperature of 60 C. Linear variation with temperature of both these parameters has been adopted. In numerical simulation the PMMA is assumed as a linear elastic material. In the experiment, operating at a room temperature of 25 C, a load was applied 3 mm from the interface in the PMMA zone (Fig. 11) in order to trigger the fracture process. The loading rate was very low, so that quasi-static conditions could be assumed. The crack path was obtained and stresses near the crack tip in the


457

Fig. 11. Geometry of the three-point bending test for a bi-material specimen (length in mm).

Table 2 Mechanical and thermal properties of the materials at 25 C

Al6061 PMMA

YoungÕs modulus (GPa)

Poisson coefficient

Thermal exp. coeff. (C1)

Thermal conductivity (W/mK)

Fracture energy (MPa m1)

r0 (MPa)

69 3.338

0.33 0.35

2.36 · 105 7.45 · 105

167 0.2

320.0

50.0

Table 3 Dependence of PMMA YoungÕs modulus on temperature

YoungÕs modulus (GPa)

25 C

35 C

45 C

60 C

3.4

3.4

2.6

2.05

PMMA were measured using a shearing interferometer. The experiment was then repeated when the temperature of the aluminium was 60 C in steady state conditions. To reach these conditions, a cartridge heater (Q in Fig. 11) was inserted into the aluminium part near the external vertical side. The variation in time of the PMMA temperature was checked before the fracture test, which was performed when steady state conditions were reached [18]. In the two experiments, the difference in crack propagation trajectory is remarkable as shown in Fig. 12a and b. In particular the crack path is closer to the interface when the temperature is higher. A highly

Fig. 12. Crack paths: (a) T = 25 C, (b) T = 60 C experimental (reproduced from Ref. [2]), (c) numerical. Case A: uniform temperature (25 C); Case B: thermal load with E, r0, drcr varying with temperature; Case C: thermal load with E = E (25 C), r0 = r0 (25 C), drcr = drcr (25 C).

458


Table 4 Comparison for crack initiation angle: Case A: uniform temperature (25 C); Case B: thermal load with E, r0, drcr varying with temperature; Case C: thermal load with E = E (25 C), r0 = r0 (25 C), drcr = drcr (25 C) Case

A B C

Crack initiation angle [] Experimental

Numerical

24 13 –

25 12 25

representative parameter is the crack initiation angle, which experimental and numerical values are compared in Table 4. Further details for the numerical simulation can be found in Ref. [18]. Qualitative comparison of the present solution with the experimental results is performed on the basis of crack path. Fig. 12c presents calculated crack paths in three cases. Case A considers uniform temperature (25 C) in conjunction with a mechanical load of 835 N. In case B the effects of the thermal load due to heating and vertical load of 650 N are analysed, assuming the PMMA YoungÕs modulus, limiting stress r0 and critical opening drcr varying with temperature. Case C is equal to the previous one, except for parameters E, r0 and drc which are assumed constant and corresponding to the ones at room temperature. The obtained crack path configurations are of the same nature of the experimental ones and compare quite well with them (Fig. 12a and b), in particular showing that the path is closer to the interface when the temperature is higher. The agreement is corroborated by the crack initiation angle values listed in Table 4: numerical and experimental results are in good accordance. The same happens for the temperature field as shown in Fig. 13, where temperature distribution in PMMA versus distance from the interface is presented. A steep gradient near the interface appears, which is correctly caught by the numerical solution. Previously presented results are obtained by using the standard refinement techniques both for the displacement and temperature field. The good accordance with the laboratory experiment demonstrates the sufficiency of this approach. The steady state temperature conditions for aluminium in which the laboratory experiment was performed result in an almost steady temperature in PMMA. From the numerical point of view, this fact reduces the importance and the difficulty of the time integration. To analyse the importance of the goal-oriented h-refinement, the discretization of the process zone is now analysed in detail and the element threshold number is determined. Fig. 14 presents the contour

Fig. 13. Temperature distribution of PMMA versus distance from the interface.


Fig. 14. Contour map of maximum principal stress and cohesive forces at fracture tip (Case B).

Fig. 15. External force versus vertical displacement and mesh size.

459

460


map of the maximum principal stress and cohesive force distribution in an initial stage of fracture propagation. The calculated length of the process zone is very limited. The maximum length is nearly 0.8 mm at the beginning and decreases during the evolution of the phenomenon. This value is in good accordance with BarenblattÕs [20] cohesive theory, which yields for PMMA ‘¼

pEGc ffi 0.75 mm. 8ð1 m2 Þr20

ð9Þ

To obtain a mesh independent solution it is necessary to properly discretize the process zone, hence severe limits are introduced depending on the approximation of the used elements. An heuristic assessment has been made on the influence of the mesh size in the process zone: using linear elements of decreasing size, the value of the force P (Fig. 11), corresponding to an applied vertical displacement on the same point, is calculated. Results are summarized in Fig. 15. As can be seen the peak of the external load and the softening branch result independent of discretization when the process zone is subdivided into at least five elements with edges of 0.15 mm or smaller. This situation is handled by the mesh generator simply by locating an element source in correspondence of the crack tip. Its weight may be a priori stated (using Eq. (9)) and/ or can be a posteriori updated during the adaptive remeshing procedure once the real length of the process zone has been determined. It is interesting to remark that only the use of and adaptive procedure can handle such type of problems. The unique alternative is a very fine mesh in large areas, resulting in a huge numerical model. On the contrary the use of interface elements is prevented by the fact that crack path is a priori unknown.

6. Conclusions The paper discusses the application of the existing adaptive refinement techniques to multi-field problems including cohesive fracture. Capabilities and shortcomings of this application are identified, the latter depending on the presence of coupling terms. The importance of a proper simultaneous space/time discretization is analysed in order to obtain correct results, even though some quantities are nearly independent of very fine discretizations. The most important conclusion, however, regards the existence of an element threshold number to get independent mesh results, when using cohesive fracture models. This number is obtained by means of a goal-oriented adaptive refinement of the process zone, in which the quantity of interest is chosen far away from the process zone. The importance of properly determining the crack tip advancement and velocity in multi-physics problems when at least one field is time dependent clearly appears from this paper.

Acknowledgements This research was partially supported by the Italian Ministero dellÕIstruzione, dellÕUniversita` e della Ricerca (PRIN 2002087915_004) and by European Network NW-IALAD.

References [1] R.W. Lewis, B.A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, Wiley, 1998. [2] D. Garagash, E. Detournay, The tip region of a fluid-driven fracture in an elastic medium, J. Appl. Mech. 67 (2000) 183–192.


461

[3] N.C. Huang, S.G. Russel, Hydraulic fracturing of a saturated porous medium—II: special cases, Theor. Appl. Fract. Mech. 4 (1985) 215–222. [4] B.J. Carter, J. Desroches, A.R. Ingraffea, P.A. Wawrzynek, Simulating fully 3-D hydraulic fracturing, in: M. Zaman, J. Booker, G. Gioda (Eds.), Modeling in Geomechanics, Wiley, 2000, pp. 525–567. [5] P.A. Witherspoon, J.S.Y. Wang, K. Iwai, J.E. Gale, Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res. 16 (1980) 1016–1024. [6] A. Carpinteri, S. Valente, F.P. Zhou, G. Ferrara, G. Melchiorri, Crack propagation in concrete specimen subjected to sustained loads, in: F.H. Wittmann (Ed.), Proceedings of the Second International Conference on Fracture Mechanics of Concrete Structures, FRA.M.CO.S. II, 1995, pp. 1315–1328. [7] O.C. Zienkiewicz, R.L. TaylorThe Finite Element Method, vol. 1–3, Butterworth-Heinemann, Oxford, 2000. [8] N.E. Wieberg, X.D. Li, A postprocessed error estimate and an adaptive procedure for the semidiscrete finite elements method in dynamic analysis, Int. J. Numer. Methods Eng. 37 (1994) 3585–3603. [9] A. Rieger, P. Wriggers, Adaptive methods for thermomechanical coupled contact problems, Int. J. Numer. Methods Engrg. 59 (2004) 871–894. [10] J.Z. Zhu, O.C. Zienkiewicz, Adaptive techniques in the finite element method, Commun. Appl. Numer. Methods 4 (1988) 197– 204. [11] S. Secchi, L. Simoni, An improved procedure for 2-D unstructured Delaunay mesh generation, Adv. Engrg. Software 34 (2003) 217–234. [12] L. Simoni, S. Secchi, Cohesive fracture mechanics for a multi-phase porous medium, Engrg. Comput. 20 (2003) 675–698. [13] L.G. Margolin, A generalized Griffith criterion for crack propagation, Engrg. Fract. Mech. 19 (1984) 539–543. [14] J.K. Dienes, Comments on a generalized Griffith criteria for crack propagation, Engrg. Fract. Mech. 23 (1986) 615–617. [15] A. Hilleborg, M. Modeer, P.E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, Cement Concrete Res. 6 (1976) 773–782. [16] T.J. Boone, A.R. Ingraffea, A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media, Int. J. Numer. Anal. Methods Geomech. 14 (1990) 27–47. [17] G.T. Camacho, M. Ortiz, Computational modelling of impact damage in brittle materials, Int. J. Solids Struct. 33 (1996) 2899– 2938. [18] S. Secchi, L. Simoni, B.A. Schrefler, Cohesive fracture growth in a thermoelastic bi-material medium, Comput. Struct. 82 (2004) 1875–1887. [19] F. Zhou, J.F. Molinari, Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency, Int. J. Numer. Methods Engrg. 59 (2004) 1–24. [20] G.I. Barenblatt, The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses axially-symmetric cracks, J. Appl. Math. Mech. 23 (1959) 622–636. [21] ICOLD, Fifth International Benchmark Workshop on Numerical Analysis of dams, Them A2, Denver, Colorado. [22] J.-S. Bae, S. Krishnaswamy, Subinterfacial cracks in bimaterial systems subjected to mechanical and thermal loading, Engrg. Fract. Mech. 68 (2001) 1081–1094.