- Email: [email protected]
Iterative simulation of 3D heat diffusion in a medium with multiple cracks
Engineering Analysis with Boundary Elements 41 (2014) 10–17
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Iterative simulation of 3D heat diffusion in a medium with multiple cracks C. Serra a,n, A. Tadeu b, N. Simões b a b
ITeCons—Instituto de Investigação e Desenvolvimento Tecnológico em Ciências da Construção, Rua Pedro Hispano s/n., 3030-289 Coimbra, Portugal Department of Civil Engineering, University of Coimbra, Pólo II, Rua Luís Reis Santos, 3030-788 Coimbra, Portugal
art ic l e i nf o
a b s t r a c t
Article history: Received 19 August 2013 Accepted 28 December 2013 Available online 25 January 2014
This paper presents an iterative three-dimensional (3D) normal-derivative equation model (TBEM) to simulate 3D heat diffusion generated by a point heat source in the presence of 3D cracks embedded in an unbounded spatially uniform solid medium. The method is intended to reduce the processing time (CPU time) needed to compute 3D heat diffusion using a TBEM formulation. In the proposed formulation each inclusion is modelled individually and successively. The first crack is submitted to an incident heat field and produces a disturbance. Each one of the cracks analyzed next is reached by a heat field generated by the previous one, which is seen as an incident field. The iterative process is stopped when the heat field disturbance generated by each inclusion is negligible. The final solution is the sum of all the contributions (disturbances in the heat field). Performance of the iterative approach proposed in this study is evaluated by comparing results generated using the full 3D TBEM and using the iterative model, in terms of temperature results, CPU time required for a given frequency, as well as number of iterations. The applicability of the proposed method is illustrated via a numerical example of heat field in time domain computation. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Point heat source 3D heat diffusion Numerical modelling Normal-derivative equation Iterative TBEM Multi-inclusions
1. Introduction Numerical modelling of heat diffusion in media with embedded inclusions has been used to evaluate the influence that defects have on temperature distribution and heat fluxes inside the domain, which has proved to be useful in several areas. In particular, numerical models can be useful in establishing the limits of the effectiveness of using infrared-thermography (IRT) in defect detection studies, as stated by Maldague [1]. In such studies, the numerical modelling of inclusions with different geometries can indicate defect detectability and determine the optimum observation time window of the best subsurface defect visibility without the expense of making and testing specimens. Solving transient heat transfer problems is useful for the detection and quantitative characterization of the properties of subsurface defects, but the 3D nature of subsurface defects, plus the need to simulate heat transfer and diffusion phenomena in transient regime, presents challenges for researchers [2,3]. Conventionally, numerical modelling requires either the discretization of the domain of the problem, as in finite elements methods (FEM) [4] and finite differences methods (FDM) [5,6], or
n
Corresponding author. Tel.: þ 351 910304560. E-mail addresses: [email protected], [email protected] (C. Serra), [email protected] (A. Tadeu), [email protected] (N. Simões). 0955-7997/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.12.010
the discretization of the boundary in the case of the boundary element method (BEM) [7,8]. Researchers have recently been developing meshless formulations which require neither domain nor boundary discretization, such as the method of fundamental solutions (MFS) [9]. Methods based on domain discretization are however better suited to deal with a bounded domain, since they require its full discretization. Furthermore, in the case of multiple inclusions the space between them also needs to be discretized, which requires special attention. This becomes unfeasible if there is a very large number of inclusions embedded in an unbounded domain. Of the available numerical methods for homogeneous unbounded or semi-infinite systems modelling, the BEM is the technique that automatically satisfies far field conditions and therefore only requires the discretization of the boundaries of the inclusions. FDM and FEM techniques lead to sparse systems of equations, while the BEM results in fully populated systems of equations. Another drawback of the BEM is that it can only be applied to more general geometries and media when the relevant fundamental solutions are known. In addition, it is well known that when using the BEM the boundary integrals may become singular or nearly singular, depending on the distance between the source point and the node being integrated. Also, when the thickness of the heterogeneity being modelled tends towards zero, as in the case of delaminations, cracks or thin defects, the conventional direct BEM degenerates and becomes inaccurate, and is no longer a valid basis for numerical modelling. Among the
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
techniques that have been proposed to overcome this is the dual boundary element method (DBEM) or the normal-derivative integral equation (TBEM), which leads to hypersingular integrals. The correct integration of singular and hypersingular integrals is one of the big challenges of these techniques. A number of approaches have been proposed to deal with hypersingular integrals that arise in DBEM [10]. Solutions for specific 2D problems can be found in Cruse [11], Sládek and Sládek [12], Prosper [13], Prosper and Kausel [14], Amado Mendes and Tadeu [15]. Tadeu et al. [17] described an analytical evaluation of the singular and hypersingular integrals that appear in 3D boundary element formulations for heat diffusion, in the frequency domain. More recently, Matsumoto et al. [16] presented a solution for exterior acoustic problems governed by the Helmholtz equation using the normal derivative Burton–Miller formulation which leads to hypersingular integrals, by using constant triangular elements to discretize the boundary. Nevertheless, the BEM is still regarded as one of the most suitable tools for modelling 3D heat diffusion generated by heat sources in an unbounded spatially uniform solid medium. However, these simulations are usually associated with a high computational effort. Various BEM formulations aimed at reducing it have been proposed. Ma et al. [18] used a BEM formulation to study transient heat conduction in 3D solids with fiber inclusions. Jablonski [19] proposed the analytical evaluation of the BEM surface integrals by means of 3D Laplace and Poisson equations. Qin et al. [20] implemented changes to the conventional distance transformation technique to evaluate nearly singular integrals on 3D boundary elements, including planar and curved surface elements and very irregular elements of slender shape. Iterative solvers or multi-domain methods have also been suggested by authors to solve/avoid systems in 3D BEM problems with very large meshes. As described by Škerget and Rek [21] and Ramšak et al. [22], by using a multi-domain method the system matrix becomes sparse, thus considerably decreasing the amount of memory needed. By analyzing subdomains separately, where independent discretization can be considered for each subdomain and suitable solvers can be used for their individual systems of equations, smaller and better conditioned systems of equations can be obtained. One disadvantage of 3D multi-domain BEM stems from the difficulties of applying interface conditions between subdomains with a highly sparse system matrix. To overcome this, Ramšak and Škerget have extended their previous work using discretization of mixed elements [21,22] to 3D problems, and proposed a 3D BEM formulation for very large meshes using a multi-domain method where each element is itself a subdomain [23]. Valente and Pina [24] have explored iterative techniques based on conjugate gradient type methods as an alternative to the direct solution techniques for large scale three-dimensional problems. They concluded that these methods are competitive for BEM algebraic systems of equations, especially if used with an appropriate preconditioner [25]. Researchers such as Marburg and Schneider [26], Ylä-Oijala and Järvenpää [27] and Alia et al. [28] have proposed iterative approaches for acoustics problems. In this paper we describe an iterative approach to simulate 3D heat diffusion in the presence of multi-inclusions using the BEM formulated in the frequency domain. Null thickness inclusions are dealt with by means of a normal-derivative equation formulation. Because only one inclusion is solved at each step, the matrix storage requirements are reduced, as a smaller system of equations is generated. In the first iteration each inclusion is solved considering the prescribed boundary conditions and disregarding the other inclusions. In each of the subsequent iterations only the heat field disturbance generated by the inclusions modelled previously (viewed as an incident field) is considered. Since the coefficient matrixes remain the same throughout the process, the systems of equations are only solved during the first iteration.
11
The iterative process is stopped when the new heat field generated by each inclusion is sufficiently small. This is achieved by defining a convergence criterion which consists of comparing the responses obtained in specific receivers from each iteration. Performance of the proposed iterative 3D TBEM is studied by comparing the results generated by the full 3D normal-derivative equation formulation with those obtained using the iterative model. The number of iterations and the CPU time taken to compute the numerical responses are used to evaluate the computational efficiency of the proposed iterative formulation. In the sections that follow, first the problem is defined and the iterative boundary element formulation used is presented (an iterative approach to the normal-derivative integral equations— 3D TBEM—formulated in the frequency domain). Analytical solutions are used to solve the hypersingular integrals that appear in the 3D TBEM formulation when the element being integrated is the loaded element. The performance of the proposed iterative method is studied. The method employed to obtain time-domain responses from frequency-domain calculations is also described. Finally, a numerical application is presented to illustrate the applicability and usefulness of the proposed approach.
2. Numerical formulation 2.1. Problem definition In this paper the three-dimensional (3D) heat diffusion generated by a point heat source in the presence of 3D cracks is modelled. The cracks are embedded in an unbounded spatially uniform solid medium of density ρ, thermal conductivity λ and specific heat c. To illustrate the numerical formulation used, consider two 3D cracks buried in a uniform medium, with surfaces S1 and S2 , subjected to a point heat source O, as shown in Fig. 1. The cracks are assumed to be of null thickness, and null heat fluxes are prescribed along their surfaces. The harmonic point heat source O placed at xs ¼ ðxs ; ys ; zs Þ generates an incident heat field in x ¼ ðx; y; zÞ which is expressed by pffiffiffiffiffiffiffiffiffiffiffiffiffiffi P e i ðiω=KÞr0 ; ð1Þ T inc ðx; xs ; ωÞ ¼ 2λr 0 in which T inc is the incident heat field at x when the point heat source is located at xs and oscillates at a frequency of ω, K is the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi thermal diffusivity defined by λ=ρ c, i ¼ 1, r 0 ¼ ðx xs Þ2 þ ðy ys Þ2 þ ðz zs Þ2 and P is the amplitude of the heat source. 2.2. The iterative 3D normal derivative integral equation (iterative 3D TBEM) This section describes the iterative 3D TBEM formulation used to simulate 3D heat diffusion in the presence of multiple cracks, generated by a point heat source.
Fig. 1. 3D view of the geometry of the problem.
12
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
we get a T ð0Þ ðx0 ; ωÞ ¼
Z S1
Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds þ T inc ðx0 ; nn2 ; xs ; ωÞ:
ð5Þ
Green0 s functions H are defined by applying the gradient operator to H. These can be seen as the derivatives of the former Green0 s functions to obtain heat fluxes. In these equations, nn2 is the unit outward, normal to the boundary at the collocation points x0 . In this equation, the factor a is null for piecewise planar boundary elements. The required 3D Green0 s functions for an unbounded space are now defined as Hðx; nn1 ; nn2 ; x0 ; ωÞ ¼
∂H ∂x ∂H ∂y ∂H ∂z þ þ ; ∂x ∂nn2 ∂y ∂nn2 ∂z ∂nn2
ð6Þ
with
n h i h io ∂r 2 ∂x ∂r ∂r ∂y ∂r ∂r ∂z ∂x ; A ∂x ∂nn1 þ ∂x ∂y ∂nn1 þ ∂x ∂z ∂nn1 þ B ∂nn1
2 h i ∂y ∂y ∂H 1 ∂r ∂r ∂x ∂r ∂r ∂r ∂z ; ∂nn1 þ ∂y ∂z ∂nn1 þ B ∂nn1 ∂y ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π A ∂x ∂y ∂nn1 þ ∂y n h i h io 2 ∂H 1 ∂r ∂r ∂x ∂r ∂r ∂y ∂r ∂z ∂z ; ∂z ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π A ∂x ∂z ∂nn1 þ ∂y ∂z ∂nn1 þ ∂z ∂nn1 þ B ∂nn1 ∂H 1 ∂x ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π
A¼
Fig. 2. Iteration 0—Step 1: (a) geometry of the problem; and (b) discretization of the first crack: nodal points and boundary elements.
Iteration 0—Step 1: The incident heat field reaches the first crack while the second crack is assumed to be absent (see Fig. 2a)
The temperature T at any point x of the spatial domain can be calculated by the diffusion equation,
∂2 ∂2 ∂2 þ 2þ 2 2 ∂x ∂y ∂z
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi in which x ¼ ðx; y; zÞ and kc ¼ iω=K . Since the cracks being modelled have null thickness, a normalderivative integral equation is used to overcome the limitations of the classical boundary element method. The normal derivative integral equation can be derived by applying the gradient operator to the boundary integral equation, b T ð0Þ ðx0 ; ωÞ ¼
Z S1
Hðx; nn1 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds þ T inc ðx0 ; xs ; ωÞ:
ð3Þ
where H are the fundamental solutions (Green0 s functions) for the heat fluxes at x, due to a virtual point heat load at x0 ¼ ðx0 ; y0 ; z0 Þ and nn1 is the unit outward normal along the boundary of the crack at x. T inc is the incident heat field at x0 when the point heat source is located at xs . The factor b is a constant defined by the shape of the boundary which takes the value of 1/2 if x0 A S1 and 1 otherwise. The superscript used in T ð0Þ indicates the number of the iteration (0). Green0 s function for heat flux in an unbounded medium, in Cartesian coordinates, is given by eikc r ð ikc r 1Þ ∂r ; ð4Þ ∂nn1 4λπr 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with r ¼ ðx x0 Þ2 þ ðy y0 Þ2 þ ðz z0 Þ2 . By applying the gradient operator to Eq. (3) and assuming the inclusion is loaded with dipole heat sources (dynamic doublets), Hðx; nn1 ; x0 ; ωÞ ¼
and
ikc r
B ¼ ikc er2
ikc r
e r3 :
In Eq. (5) the incident heat field is computed as Peikc r0 ð ikc r 0 1Þ ∂r 0 ∂x ∂r 0 ∂y ∂r 0 ∂z T inc ðx0 ; nn2 ; xs ; ωÞ ¼ þ þ : ∂x ∂nn2 ∂y ∂nn2 ∂z ∂nn2 2r 0 2 ð7Þ The solution is found by solving Eq. (5), which requires the discretization of the crack S1 into N 1 planar boundary elements, with one nodal point at the center of each element (see Fig. 2b). This leads to a system of ½N 1 N 1 equations (B T ð0Þ ¼ T ð0Þ ). B is the inc R matrix that compiles the integration of S1 Hðx; nn1 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds along each boundary element, for each virtual load placed along the nodal points. This system of equations can expressed as kl
Tðx; ωÞ þ ðkc Þ2 Tðx; ωÞ ¼ 0;
ikc r ikc r k2c e ikc r þ 3ikc er2 þ 3e r3 r
ð0Þk
½H ½T ð0Þl ¼ ½T inc ;
ð8Þ
where k; l ¼ 1 to N 1 (l identifies the element being integrated and k R kl corresponds to the element being loaded), H ¼ Al Hðxl ; nn1 ; nn2 ; ð0Þk
xk ; ωÞdAl , T ð0Þl is the nodal temperature at the iteration ð0Þ, T inc is the incident temperature field at the iteration ð0Þand Al is the area of each boundary element. The required integrations in these equations are performed using a Gaussian quadrature scheme when the element to be integrated is not the loaded element. For the loaded element, the resulting hypersingular integrands of the Green0 s functions are calculated analytically [17]. The solution to the system of Eq. (8) gives the nodal heat T ð0Þ along the boundary of S1 , which allows the heat diffusion field to be defined anywhere in the field xrec as Z T 01 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ð0Þ ðx; ωÞ ds: ð9Þ S1
In this equation, the subscripts in T 01 ðxrec ; ωÞ indicate the number of the iteration (0) and the crack which produces it (1). Iteration 0—Step 2: The second crack is reached by the incident heat field and by the field generated by the first crack in Step 1 (see Fig. 3a). The second null thickness crack is modelled in the same way as the first one. S2 is discretized into N 2 planar boundary elements with one nodal point at the center of each element. Now the heat field generated by the crack S1 is taken into account. The heat field generated by the first crack defined in Step 1 is viewed as an
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
13
Fig. 4. Iteration k—Step 1.
Fig. 3. Iteration 0—Step 2: (a) geometry of the problem; and (b) discretization of the second crack: nodal points and boundary elements.
incident heat field that strikes the second crack R S1 Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞds, leading to Z a T ð0Þ ðx0 ; ωÞ ¼ Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds
which leads to a T ðkÞ ðx0 ; ωÞ ¼
S2
ð0Þ
þ T inc ðx0 ; nn2 ; xs ; ωÞ þ T 12 ðx0 ; nn2 ; ωÞ:
This leads to a system of ½N 2 N2 equations (C T ð0Þ ¼ T ð0Þ ) inc which allows the unknown nodal amplitudes to be defined. The system is expressed as ð0Þk
½H ½T ð0Þl ¼ ½T inc ;
ð11Þ ð0Þ
where k; l ¼ 1 to N 2 and T ð0Þ inc ðx0 ; nn2 ; ωÞ ¼ T inc ðx0 ; nn2 ; xs ; ωÞ þ T 12 ðx0 ; nn2 ; ωÞ. The disturbance in the heat field at xrec can then be expressed as Z T 02 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ð0Þ ðx; ωÞ ds: ð12Þ
Z S1
Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞ ds
ðk 1Þ
ð10Þ
Once again the a factor is null for piecewise planar boundary elements.
kl
Fig. 5. Iteration k—Step 2.
ð0Þ T 12 ðx0 ; nn2 ; ωÞ ¼
þ T 21
ðx0 ; nn2 ; xn_ext ; ωÞ
ð15Þ
The solution of this equation leads to a system of ½N 1 N 1 ðkÞ equations (B T ðkÞ ¼ T inc ), similar to the one calculated previously, where only the constant matrix needs to be modified. In this equation (k) identifies the number of the iteration. If in iteration 0 the system has been solved by defining its inverse matrix B 1 , the new solution will not require the system to be solved, as ðkÞ
T ðkÞ ¼ B 1 T inc . The heat field disturbance at the receiver xrec can now be computed as Z T k1 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ðkÞ ðx; ωÞ ds:
ð16Þ
S1
S2
At the end of this iteration the total heat field recorded at the receiver will be M
Tðxrec ; ωÞ ¼ T inc ðxrec ; xs ; ωÞ þ ∑ T 0 m ðxrec ; ωÞ; m¼1
ð13Þ
where M is the number of inclusions (in this case M ¼ 2). Iteration k—Step 1: The first crack is reached only by the field generated by the second crack (in the conditions defined by Step 2 of iteration k 1) (see Fig. 4). At this stage the incident heat field corresponds to the field generated by the second inclusion in the previous iteration and is given by Z ðk 1Þ T 21 ðx0 ; nn2 ; ωÞ ¼ H ðx; nn1 ; nn2 ; x0 ; ωÞ T ðk 1Þ ðx; ωÞ ds; ð14Þ S2
Iteration k—Step 2: The second crack is reached only by the heat field generated by the first crack (subjected to the incident heat field in Step 1) (see Fig. 5). The heat field generated by the first crack in the previous step is assumed to be the only heat field which reaches the second R ðkÞ crack, T 12 ðx0 ; nn2 ; ωÞ ¼ S1 Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞds, leading to Z ðkÞ a T ðkÞ ðx0 ; ωÞ ¼ Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞ ds þ T 12 ðx0 ; nn2 ; ωÞ: S2
ð17Þ ), This leads to a system of ½N2 N 2 equations (C T ðkÞ ¼ T ðkÞ inc similar to the one defined before, in Eq. (15), where the constant ðkÞ
matrix needs to be replaced by T ðkÞ inc ðx0 ; nn2 ; ωÞ ¼ T 12 ðx0 ; nn2 ; ωÞ.
14
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
Values of T ðkÞ can be obtained since T ðkÞ ¼ C 1 T ðkÞ . The new heat inc field produced by this inclusion recorded at xrec is expressed by Z ð18Þ T k2 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ðkÞ ðx; ωÞ ds: S2
At the end of iteration k the total heat field recorded at the receiver will be k
Tðxrec ; ωÞ ¼ T inc ðxrec ; xs ; ωÞ þ ∑
M
∑ T km ðxrec ; ωÞ:
k¼0m¼1
ð19Þ
The iterative process continues until the contribution of the last new field to the total heat field at a specific receiver reaches a predefined threshold. The proposed iterative process requires only the individual linear system of equations related to each inclusion to be solved. Given the example used to illustrate the algorithm procedure, for two inclusions, this requires solving two individual systems of ½N 1 N 1 and ½N 2 N 2 equations only once. The full modelling process would require solving a system of ½ðN1 þ N 2 Þ ðN 1 þ N 2 Þ equations. The proposed iterative approach is more relevant as the number of cracks increases. The size of the system of equations used by the full model becomes greater when compared with the system associated with each inclusion (as used in the proposed iterative model), which raises the interest in using this approach. 2.3. Heat field in the time domain The heat field in the time domain is determined by applying a numerical inverse fast Fourier transform in the frequency domain.
Aliasing phenomena are dealt with by introducing complex frequencies with a small imaginary part, taking the form ωc ¼ ω iη (where η ¼ 0:7Δω and Δω is the frequency step). This shift is subsequently taken into account in the time domain by applying an exponential window (eη t ) to the response. The source can have any time variation and the frequency domain solution can be determined by applying a time Fourier transformation to a frequency range from f¼ 0.0 Hz up to high frequencies. However, we do not need to compute the highest frequencies in the range since the heat response falls rapidly at higher frequencies. The frequency f¼0.0 Hz corresponds to the static response. The use of complex frequencies allows this response to be computed since the argument of the functions in the integral equations is iη, i.e., other than zero.
3. Results and discussion 3.1. Case studies The proposed model is used to simulate heat diffusion around four 3D null thickness cracks, represented by the planes in Fig. 6. The inclusions are embedded in a uniform unbounded solid medium of density ρ¼2300 kg m 3, thermal conductivity λ¼ 1.40 W m 1 1C 1 and specific heat c¼880 J kg 1 1C 1. The point heat source O is located at the origin ðxs ; ys ; zs Þ ¼ ð0 m ; 0 m ; 0 mÞ. Temperature computations are recorded on two fine grids of receivers placed over 0.5 0.4 m2 planes: Grid G1 which is placed away from the source along the z axis, parallel to the xy plane, at z¼0.08 m and grid G2 which is placed parallel to the xz plane at y¼ 0.0 m, as seen in Fig. 6. Each grid contains a total of 4141 receivers spaced at equal intervals of Δx¼0.005 m, Δy¼0.01 m and Δz¼ 0.007 m. As shown in Fig. 7, four 0.2 0.2 m2 cracks are placed parallel to the yz plane, away from the source O along the x axis, 0.1 m apart from each other. Each plane is discretized into 800 elements: 40 along the z axis by 20 along the y direction. 3.2. Simulation results
Fig. 6. 3D view of the systems modelled.
In this section, the performance of the proposed iterative 3D TBEM is studied. This is done by comparing the results obtained for modelling heat diffusion in the presence of 3D cracks using the full 3D normal-derivative equation formulation (3D TBEM) with those given by using the iterative approach. All computations where performed using a 2.30 GHz Intels Core™ i5-2410M CPU with 4.00 GB RAM. The time taken to compute the generated heat field (CPU time) is compared using three distinct approaches: (I) using the full 3D TBEM to calculate heat diffusion in the cracked medium; (II) defining two sub-domains by pairing two cracks and using the iterative 3D TBEM for each subdomain; and (III) using the iterative 3D TBEM and considering each crack separately.
Fig. 7. 2D view of the geometry of the systems modelled: (a) yz plane; (b) xy plane and (c) xz plane.
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
The number of iterations needed to solve the problem using the proposed iterative model was also recorded. The iteration loop is stopped and the number of iterations (NI ) is registered when the temperature at each individual receiver does not change by more than a predefined threshold (ε) when compared to the former
15
iteration, according to the following equation:
NI M
NI M NI 1 M
T ðx ; ωÞ T ðx ; ωÞ ∑ ∑ km rec
= ∑ ∑ T iter
∑ ∑ km rec
k ¼ 0 m ¼ 1
k ¼ 0 m ¼ 1 k¼0m¼1
Fig. 8. Temperature responses, CPU time and number of iterations recorded in G2.
ðx ; ωÞ
r ε: m rec
ð20Þ
16
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
Fig. 9. Snapshots of heat field in time domain taken at (a) t ¼ 5.0 h, (b) t ¼10.6 h and (c) t¼ 12.8 h.
Additionally, transient temperature results are presented in order to illustrate heat diffusion in the time domain simulation in the cracked medium. 3.2.1. Performance of the iterative 3D TBEM Three excitation frequencies have been selected to illustrate the main findings: f ¼0.0 Hz, f ¼0.5 10 6 Hz and f¼ 6.4 10 5 Hz. Fig. 8 shows the real and imaginary temperature responses recorded in grid G2 for each of the approaches, along with the required CPU time and number of iterations. The results clearly indicate that, for any frequency, using an iterative method greatly reduces CPU time, when compared to the full 3D TBEM method. The performance of the iterative method is slightly better when each crack is considered individually (approach III) as CPU time is further reduced. The results show that the number of iterations (and the CPU time) varies for each receiver. Fewer iterations are required in zones which are less affected by disturbances in the heat field, which, for approach III (iterative 3D TBEM considering each crack separately), is the zones behind the cracks placed at x ¼0.2 m and x ¼0.4 m. In the case of the two sub-domains approach (II), it is the zone before the first crack and after the last crack. In both iterative approaches, the number of iterations diminishes as the frequency increases, however CPU time increases slightly. 3.2.2. Heat field in time domain Transient temperature computations are performed in the frequency range ½0:0; 1:024 10 3 Hz with a frequency increment of Δf¼0.5 10 6 Hz, which results in a time window of 1/(0.5 10 6) s. The imaginary part of the frequency is given by η ¼ 0:7 2 π Δf . The host medium is assumed to be initially at 20.0 1C. The source starts emitting energy at instant t¼0.5 h and continues for 9.5 h.
The heat source time dependence is assumed to be rectangular with an amplitude of P¼6.27. This amplitude is defined so that a maximum temperature increase of 15.0 1C is registered by a receiver located in grid G2 at ðx; y; zÞ ¼ ð0:1 m; 0 m; 0 mÞ when the medium is not cracked. In order to illustrate the main results, snapshot images taken at different instants are shown in Fig. 9. Fig. 9a clearly illustrates the effect that the first crack has on temperature diffusion. By instant t¼10.6 h (see Fig. 9b) the heat source has stopped emitting energy and the temperature at ðx; y; zÞ ¼ ð0:1 m; 0 m; 0 mÞ, immediately in front of the center of the first crack, reaches its peak. The snapshot in Fig. 9c shows the instant when the temperature is maximum in the zone immediately before the second crack.
4. Conclusions This paper presents an iterative model to simulate 3D heat diffusion in an unbounded spatially uniform solid medium with embedded 3D inclusions. Heat diffusion is generated by a point heat source and the resulting heat field in the presence of multiple cracks is computed using a normal-derivative integral equation (3D TBEM) formulation. An iterative approach to the 3D TBEM formulation has been described. The method to obtain time-domain responses from frequency-domain calculations has also been described. The performance of the proposed iterative 3D TBEM was assessed by comparing results generated using the full 3D normal-derivative equation formulation (3D TBEM) and those using the iterative model. Two different iterative approaches were considered: solving each crack separately or considering sub-domains by pairing cracks and solving two at each turn. By solving each crack separately, only one inclusion is computed at a time, at each step of the iterative approach. Employing this method required solving a
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
lower number of small systems of equations rather than one large system, and led to a significant decrease in CPU time (for both iterative approaches). For a given frequency, the CPU time and the number of iterations required for each receiver varied in both iterative approaches. The number of iterations required fell as the analyzed frequency increased, though the average CPU time increased slightly. Considering each crack separately generated only slightly better results than when sub-domains of pairs of inclusions were considered. It is concluded that using the proposed iterative model greatly reduces the computational effort needed to simulate 3D heat diffusion generated by a point heat source in an unbounded uniform medium with multiple cracks. In order to further demonstrate the applicability of the method, a numerical application to transient temperature calculations has also been presented. Acknowledgments The research work presented herein was supported by the FEDER funds through the Operational Programme for Competitiveness Factors—COMPETE and by national funds through the FCT —Portuguese Foundation for Science and Technology, under research project PTDC/ECM/114189/2009 and was supported in part by the doctoral grant SFRH/BD/91686/2012 attributed by the FCT. This work has been framed under the Initiative Energy for Sustainability of the University of Coimbra and supported by the Energy and Mobility for Sustainable Regions—EMSURE Project (CENTRO-07-0224-FEDER-002004). References [1] Maldague X. Introduction to NDT by active infrared thermography. Mater Eval 2002;6:1060–73. [2] Tartarini P, Corticelli M, Tarozzi L. Dropwise cooling: experimental tests by infrared thermography and numerical simulations. Appl Therm Eng 2009;29:1391–7. [3] Rodríguez F, Nicolau V. Inverse heat transfer approach for IR image reconstruction: application to thermal non-destructive evaluation. Appl Therm Eng 2012;33–34:109–18. [4] Bathe KJ. Numerical methods in finite element analysis. New Jersey: PrenticeHall; 1976. [5] Ozisik MN. Finite difference methods in heat transfer. USA: CRC Press Inc.; 1994. [6] Juncu Gh. Unsteady conjugate forced convection heat/mass transfer from a finite flat plate. Int J Therm Sci 2008;47:972–84.
17
[7] Brebbia CA, Telles JC, Wrobel LC. Boundary elements techniques: theory and applications in engineering. Berlin-New York: Springer-Verlag; 1984. [8] Ochiai Y. Steady heat conduction analysis in orthotropic bodies by triplereciprocity BEM. Comput Model Eng Sci 2001;2:435–46. [9] Šarler B. Towards a mesh-free computation of transport phenomena. Eng Anal Bound Elem 2002;26:731–8. [10] Rudolphi TJ. The use of simple solutions in the regularisation of hypersingular boundary integral equations. Math Comput Model 1991;15:269–78. [11] Cruse TA. Boundary element analysis in computational fracture mechanics. Dordrecht: Kluwer Academic Publishers; 1988. [12] Sládek V, Sládek J. Transient elastodynamics three-dimensional problems in cracked bodies. Appl Math Model 1984;8:2–10. [13] Prosper D. Modeling and detection of delaminations in laminated plates (Ph.D. thesis). Cambridge: MIT; 2001. [14] Prosper D, Kause E. Wave scattering by cracks in laminated media. In: CD: advances in computational engineering and sciences, proceedings of the international conference on computational engineering and science ICES0 01. Atluri SN, Nishioka T, Kikuchi M, (editors), Tech Science Press; 2001. [15] Amado Mendes P, Tadeu A. Wave propagation in the presence of empty cracks in an elastic medium. Comput Mech 2006;38:183–99. [16] Matsumoto T, Zheng C, Harada S, Takahashi T. Explicit evaluation of hypersingular boundary integral equation for 3-D Helmholtz equation discretized with constant triangular element. J Comput Sci Technol 2010;4:194–206. [17] Tadeu A, Prata J, Simões N. Closed form integration of singular and hypersingular integrals in 3D BEM formulations for heat conduction. Math Probl Eng, 2012; Article ID 647038, 21 p, http://dx.doi.org/10.1155/2012/647038. [18] Ma F, Chatterjee J, Henry DP, Banerjee PK. Transient heat conduction analysis of 3D solids with fibre inclusions using the boundary element method. Int J Numer Methods Eng 2008;73:1113–36. [19] Jablonski P. Integral and geometrical means in the analytical evaluation of the BEM integrals for a 3D Laplace equation. Eng Anal Bound Elem 2010;34:264–73. [20] Xianyun Q, Zhang J, Xie G, Zhou F, Li G. A general algorithm for the numerical evaluation of nearly singular integrals on 3D boundary element. J Comput Appl Math 2011;235:4174–86. [21] Škerget L, Rek Z. Boundary-domain integral method using a velocity–vorticity formulation. Eng Anal Bound Elem 1995;15:359–70. [22] Ramšak M, Škerget L, Hriberšk M, Zunič Z. A multidomain boundary element method for unsteady laminar flow using stream function—vorticity equations. Eng Anal Bound Elem 2005;29:1–14. [23] Ramšak M, Škerget L. 3D multidomain BEM for solving the Laplace equation. Eng Anal Bound Elem 2007;31:528–38. [24] Valente FP, Pina HL. Iterative techniques for 3-D boundary element method systems of equations. Eng Anal Bound Elem 2001;25:423–9. [25] Van der Vorst HA, Dekker K. Conjugate gradient type methods and preconditioning. J Comput Appl Math 1988;24:73–87. [26] Marburg S, Schneider S. Performance of iterative solvers for acoustic problems. Part I. Solvers and effect of diagonal preconditioning. Eng Anal Bound Elem 2003;27:727–50. [27] Ylä-Oijala P, Järvenpää S. Iterative solution of high-order boundary element method for acoustic impedance boundary value problems. J Sound Vib 2006;291:824–43. [28] Alia A, Sadok H, Souli M. CMRH method as iterative solver for boundary element acoustic systems. Eng Anal Bound Elem 2012;36:346–50.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Iterative simulation of 3D heat diffusion in a medium with multiple cracks C. Serra a,n, A. Tadeu b, N. Simões b a b
ITeCons—Instituto de Investigação e Desenvolvimento Tecnológico em Ciências da Construção, Rua Pedro Hispano s/n., 3030-289 Coimbra, Portugal Department of Civil Engineering, University of Coimbra, Pólo II, Rua Luís Reis Santos, 3030-788 Coimbra, Portugal
art ic l e i nf o
a b s t r a c t
Article history: Received 19 August 2013 Accepted 28 December 2013 Available online 25 January 2014
This paper presents an iterative three-dimensional (3D) normal-derivative equation model (TBEM) to simulate 3D heat diffusion generated by a point heat source in the presence of 3D cracks embedded in an unbounded spatially uniform solid medium. The method is intended to reduce the processing time (CPU time) needed to compute 3D heat diffusion using a TBEM formulation. In the proposed formulation each inclusion is modelled individually and successively. The first crack is submitted to an incident heat field and produces a disturbance. Each one of the cracks analyzed next is reached by a heat field generated by the previous one, which is seen as an incident field. The iterative process is stopped when the heat field disturbance generated by each inclusion is negligible. The final solution is the sum of all the contributions (disturbances in the heat field). Performance of the iterative approach proposed in this study is evaluated by comparing results generated using the full 3D TBEM and using the iterative model, in terms of temperature results, CPU time required for a given frequency, as well as number of iterations. The applicability of the proposed method is illustrated via a numerical example of heat field in time domain computation. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Point heat source 3D heat diffusion Numerical modelling Normal-derivative equation Iterative TBEM Multi-inclusions
1. Introduction Numerical modelling of heat diffusion in media with embedded inclusions has been used to evaluate the influence that defects have on temperature distribution and heat fluxes inside the domain, which has proved to be useful in several areas. In particular, numerical models can be useful in establishing the limits of the effectiveness of using infrared-thermography (IRT) in defect detection studies, as stated by Maldague [1]. In such studies, the numerical modelling of inclusions with different geometries can indicate defect detectability and determine the optimum observation time window of the best subsurface defect visibility without the expense of making and testing specimens. Solving transient heat transfer problems is useful for the detection and quantitative characterization of the properties of subsurface defects, but the 3D nature of subsurface defects, plus the need to simulate heat transfer and diffusion phenomena in transient regime, presents challenges for researchers [2,3]. Conventionally, numerical modelling requires either the discretization of the domain of the problem, as in finite elements methods (FEM) [4] and finite differences methods (FDM) [5,6], or
n
Corresponding author. Tel.: þ 351 910304560. E-mail addresses: [email protected], [email protected] (C. Serra), [email protected] (A. Tadeu), [email protected] (N. Simões). 0955-7997/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.12.010
the discretization of the boundary in the case of the boundary element method (BEM) [7,8]. Researchers have recently been developing meshless formulations which require neither domain nor boundary discretization, such as the method of fundamental solutions (MFS) [9]. Methods based on domain discretization are however better suited to deal with a bounded domain, since they require its full discretization. Furthermore, in the case of multiple inclusions the space between them also needs to be discretized, which requires special attention. This becomes unfeasible if there is a very large number of inclusions embedded in an unbounded domain. Of the available numerical methods for homogeneous unbounded or semi-infinite systems modelling, the BEM is the technique that automatically satisfies far field conditions and therefore only requires the discretization of the boundaries of the inclusions. FDM and FEM techniques lead to sparse systems of equations, while the BEM results in fully populated systems of equations. Another drawback of the BEM is that it can only be applied to more general geometries and media when the relevant fundamental solutions are known. In addition, it is well known that when using the BEM the boundary integrals may become singular or nearly singular, depending on the distance between the source point and the node being integrated. Also, when the thickness of the heterogeneity being modelled tends towards zero, as in the case of delaminations, cracks or thin defects, the conventional direct BEM degenerates and becomes inaccurate, and is no longer a valid basis for numerical modelling. Among the
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
techniques that have been proposed to overcome this is the dual boundary element method (DBEM) or the normal-derivative integral equation (TBEM), which leads to hypersingular integrals. The correct integration of singular and hypersingular integrals is one of the big challenges of these techniques. A number of approaches have been proposed to deal with hypersingular integrals that arise in DBEM [10]. Solutions for specific 2D problems can be found in Cruse [11], Sládek and Sládek [12], Prosper [13], Prosper and Kausel [14], Amado Mendes and Tadeu [15]. Tadeu et al. [17] described an analytical evaluation of the singular and hypersingular integrals that appear in 3D boundary element formulations for heat diffusion, in the frequency domain. More recently, Matsumoto et al. [16] presented a solution for exterior acoustic problems governed by the Helmholtz equation using the normal derivative Burton–Miller formulation which leads to hypersingular integrals, by using constant triangular elements to discretize the boundary. Nevertheless, the BEM is still regarded as one of the most suitable tools for modelling 3D heat diffusion generated by heat sources in an unbounded spatially uniform solid medium. However, these simulations are usually associated with a high computational effort. Various BEM formulations aimed at reducing it have been proposed. Ma et al. [18] used a BEM formulation to study transient heat conduction in 3D solids with fiber inclusions. Jablonski [19] proposed the analytical evaluation of the BEM surface integrals by means of 3D Laplace and Poisson equations. Qin et al. [20] implemented changes to the conventional distance transformation technique to evaluate nearly singular integrals on 3D boundary elements, including planar and curved surface elements and very irregular elements of slender shape. Iterative solvers or multi-domain methods have also been suggested by authors to solve/avoid systems in 3D BEM problems with very large meshes. As described by Škerget and Rek [21] and Ramšak et al. [22], by using a multi-domain method the system matrix becomes sparse, thus considerably decreasing the amount of memory needed. By analyzing subdomains separately, where independent discretization can be considered for each subdomain and suitable solvers can be used for their individual systems of equations, smaller and better conditioned systems of equations can be obtained. One disadvantage of 3D multi-domain BEM stems from the difficulties of applying interface conditions between subdomains with a highly sparse system matrix. To overcome this, Ramšak and Škerget have extended their previous work using discretization of mixed elements [21,22] to 3D problems, and proposed a 3D BEM formulation for very large meshes using a multi-domain method where each element is itself a subdomain [23]. Valente and Pina [24] have explored iterative techniques based on conjugate gradient type methods as an alternative to the direct solution techniques for large scale three-dimensional problems. They concluded that these methods are competitive for BEM algebraic systems of equations, especially if used with an appropriate preconditioner [25]. Researchers such as Marburg and Schneider [26], Ylä-Oijala and Järvenpää [27] and Alia et al. [28] have proposed iterative approaches for acoustics problems. In this paper we describe an iterative approach to simulate 3D heat diffusion in the presence of multi-inclusions using the BEM formulated in the frequency domain. Null thickness inclusions are dealt with by means of a normal-derivative equation formulation. Because only one inclusion is solved at each step, the matrix storage requirements are reduced, as a smaller system of equations is generated. In the first iteration each inclusion is solved considering the prescribed boundary conditions and disregarding the other inclusions. In each of the subsequent iterations only the heat field disturbance generated by the inclusions modelled previously (viewed as an incident field) is considered. Since the coefficient matrixes remain the same throughout the process, the systems of equations are only solved during the first iteration.
11
The iterative process is stopped when the new heat field generated by each inclusion is sufficiently small. This is achieved by defining a convergence criterion which consists of comparing the responses obtained in specific receivers from each iteration. Performance of the proposed iterative 3D TBEM is studied by comparing the results generated by the full 3D normal-derivative equation formulation with those obtained using the iterative model. The number of iterations and the CPU time taken to compute the numerical responses are used to evaluate the computational efficiency of the proposed iterative formulation. In the sections that follow, first the problem is defined and the iterative boundary element formulation used is presented (an iterative approach to the normal-derivative integral equations— 3D TBEM—formulated in the frequency domain). Analytical solutions are used to solve the hypersingular integrals that appear in the 3D TBEM formulation when the element being integrated is the loaded element. The performance of the proposed iterative method is studied. The method employed to obtain time-domain responses from frequency-domain calculations is also described. Finally, a numerical application is presented to illustrate the applicability and usefulness of the proposed approach.
2. Numerical formulation 2.1. Problem definition In this paper the three-dimensional (3D) heat diffusion generated by a point heat source in the presence of 3D cracks is modelled. The cracks are embedded in an unbounded spatially uniform solid medium of density ρ, thermal conductivity λ and specific heat c. To illustrate the numerical formulation used, consider two 3D cracks buried in a uniform medium, with surfaces S1 and S2 , subjected to a point heat source O, as shown in Fig. 1. The cracks are assumed to be of null thickness, and null heat fluxes are prescribed along their surfaces. The harmonic point heat source O placed at xs ¼ ðxs ; ys ; zs Þ generates an incident heat field in x ¼ ðx; y; zÞ which is expressed by pffiffiffiffiffiffiffiffiffiffiffiffiffiffi P e i ðiω=KÞr0 ; ð1Þ T inc ðx; xs ; ωÞ ¼ 2λr 0 in which T inc is the incident heat field at x when the point heat source is located at xs and oscillates at a frequency of ω, K is the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi thermal diffusivity defined by λ=ρ c, i ¼ 1, r 0 ¼ ðx xs Þ2 þ ðy ys Þ2 þ ðz zs Þ2 and P is the amplitude of the heat source. 2.2. The iterative 3D normal derivative integral equation (iterative 3D TBEM) This section describes the iterative 3D TBEM formulation used to simulate 3D heat diffusion in the presence of multiple cracks, generated by a point heat source.
Fig. 1. 3D view of the geometry of the problem.
12
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
we get a T ð0Þ ðx0 ; ωÞ ¼
Z S1
Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds þ T inc ðx0 ; nn2 ; xs ; ωÞ:
ð5Þ
Green0 s functions H are defined by applying the gradient operator to H. These can be seen as the derivatives of the former Green0 s functions to obtain heat fluxes. In these equations, nn2 is the unit outward, normal to the boundary at the collocation points x0 . In this equation, the factor a is null for piecewise planar boundary elements. The required 3D Green0 s functions for an unbounded space are now defined as Hðx; nn1 ; nn2 ; x0 ; ωÞ ¼
∂H ∂x ∂H ∂y ∂H ∂z þ þ ; ∂x ∂nn2 ∂y ∂nn2 ∂z ∂nn2
ð6Þ
with
n h i h io ∂r 2 ∂x ∂r ∂r ∂y ∂r ∂r ∂z ∂x ; A ∂x ∂nn1 þ ∂x ∂y ∂nn1 þ ∂x ∂z ∂nn1 þ B ∂nn1
2 h i ∂y ∂y ∂H 1 ∂r ∂r ∂x ∂r ∂r ∂r ∂z ; ∂nn1 þ ∂y ∂z ∂nn1 þ B ∂nn1 ∂y ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π A ∂x ∂y ∂nn1 þ ∂y n h i h io 2 ∂H 1 ∂r ∂r ∂x ∂r ∂r ∂y ∂r ∂z ∂z ; ∂z ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π A ∂x ∂z ∂nn1 þ ∂y ∂z ∂nn1 þ ∂z ∂nn1 þ B ∂nn1 ∂H 1 ∂x ðx; nn1 ; nn2 ; x0 ; ωÞ ¼ 4π
A¼
Fig. 2. Iteration 0—Step 1: (a) geometry of the problem; and (b) discretization of the first crack: nodal points and boundary elements.
Iteration 0—Step 1: The incident heat field reaches the first crack while the second crack is assumed to be absent (see Fig. 2a)
The temperature T at any point x of the spatial domain can be calculated by the diffusion equation,
∂2 ∂2 ∂2 þ 2þ 2 2 ∂x ∂y ∂z
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi in which x ¼ ðx; y; zÞ and kc ¼ iω=K . Since the cracks being modelled have null thickness, a normalderivative integral equation is used to overcome the limitations of the classical boundary element method. The normal derivative integral equation can be derived by applying the gradient operator to the boundary integral equation, b T ð0Þ ðx0 ; ωÞ ¼
Z S1
Hðx; nn1 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds þ T inc ðx0 ; xs ; ωÞ:
ð3Þ
where H are the fundamental solutions (Green0 s functions) for the heat fluxes at x, due to a virtual point heat load at x0 ¼ ðx0 ; y0 ; z0 Þ and nn1 is the unit outward normal along the boundary of the crack at x. T inc is the incident heat field at x0 when the point heat source is located at xs . The factor b is a constant defined by the shape of the boundary which takes the value of 1/2 if x0 A S1 and 1 otherwise. The superscript used in T ð0Þ indicates the number of the iteration (0). Green0 s function for heat flux in an unbounded medium, in Cartesian coordinates, is given by eikc r ð ikc r 1Þ ∂r ; ð4Þ ∂nn1 4λπr 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with r ¼ ðx x0 Þ2 þ ðy y0 Þ2 þ ðz z0 Þ2 . By applying the gradient operator to Eq. (3) and assuming the inclusion is loaded with dipole heat sources (dynamic doublets), Hðx; nn1 ; x0 ; ωÞ ¼
and
ikc r
B ¼ ikc er2
ikc r
e r3 :
In Eq. (5) the incident heat field is computed as Peikc r0 ð ikc r 0 1Þ ∂r 0 ∂x ∂r 0 ∂y ∂r 0 ∂z T inc ðx0 ; nn2 ; xs ; ωÞ ¼ þ þ : ∂x ∂nn2 ∂y ∂nn2 ∂z ∂nn2 2r 0 2 ð7Þ The solution is found by solving Eq. (5), which requires the discretization of the crack S1 into N 1 planar boundary elements, with one nodal point at the center of each element (see Fig. 2b). This leads to a system of ½N 1 N 1 equations (B T ð0Þ ¼ T ð0Þ ). B is the inc R matrix that compiles the integration of S1 Hðx; nn1 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds along each boundary element, for each virtual load placed along the nodal points. This system of equations can expressed as kl
Tðx; ωÞ þ ðkc Þ2 Tðx; ωÞ ¼ 0;
ikc r ikc r k2c e ikc r þ 3ikc er2 þ 3e r3 r
ð0Þk
½H ½T ð0Þl ¼ ½T inc ;
ð8Þ
where k; l ¼ 1 to N 1 (l identifies the element being integrated and k R kl corresponds to the element being loaded), H ¼ Al Hðxl ; nn1 ; nn2 ; ð0Þk
xk ; ωÞdAl , T ð0Þl is the nodal temperature at the iteration ð0Þ, T inc is the incident temperature field at the iteration ð0Þand Al is the area of each boundary element. The required integrations in these equations are performed using a Gaussian quadrature scheme when the element to be integrated is not the loaded element. For the loaded element, the resulting hypersingular integrands of the Green0 s functions are calculated analytically [17]. The solution to the system of Eq. (8) gives the nodal heat T ð0Þ along the boundary of S1 , which allows the heat diffusion field to be defined anywhere in the field xrec as Z T 01 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ð0Þ ðx; ωÞ ds: ð9Þ S1
In this equation, the subscripts in T 01 ðxrec ; ωÞ indicate the number of the iteration (0) and the crack which produces it (1). Iteration 0—Step 2: The second crack is reached by the incident heat field and by the field generated by the first crack in Step 1 (see Fig. 3a). The second null thickness crack is modelled in the same way as the first one. S2 is discretized into N 2 planar boundary elements with one nodal point at the center of each element. Now the heat field generated by the crack S1 is taken into account. The heat field generated by the first crack defined in Step 1 is viewed as an
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
13
Fig. 4. Iteration k—Step 1.
Fig. 3. Iteration 0—Step 2: (a) geometry of the problem; and (b) discretization of the second crack: nodal points and boundary elements.
incident heat field that strikes the second crack R S1 Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞds, leading to Z a T ð0Þ ðx0 ; ωÞ ¼ Hðx; nn1 ; nn2 ; x0 ; ωÞ T ð0Þ ðx; ωÞ ds
which leads to a T ðkÞ ðx0 ; ωÞ ¼
S2
ð0Þ
þ T inc ðx0 ; nn2 ; xs ; ωÞ þ T 12 ðx0 ; nn2 ; ωÞ:
This leads to a system of ½N 2 N2 equations (C T ð0Þ ¼ T ð0Þ ) inc which allows the unknown nodal amplitudes to be defined. The system is expressed as ð0Þk
½H ½T ð0Þl ¼ ½T inc ;
ð11Þ ð0Þ
where k; l ¼ 1 to N 2 and T ð0Þ inc ðx0 ; nn2 ; ωÞ ¼ T inc ðx0 ; nn2 ; xs ; ωÞ þ T 12 ðx0 ; nn2 ; ωÞ. The disturbance in the heat field at xrec can then be expressed as Z T 02 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ð0Þ ðx; ωÞ ds: ð12Þ
Z S1
Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞ ds
ðk 1Þ
ð10Þ
Once again the a factor is null for piecewise planar boundary elements.
kl
Fig. 5. Iteration k—Step 2.
ð0Þ T 12 ðx0 ; nn2 ; ωÞ ¼
þ T 21
ðx0 ; nn2 ; xn_ext ; ωÞ
ð15Þ
The solution of this equation leads to a system of ½N 1 N 1 ðkÞ equations (B T ðkÞ ¼ T inc ), similar to the one calculated previously, where only the constant matrix needs to be modified. In this equation (k) identifies the number of the iteration. If in iteration 0 the system has been solved by defining its inverse matrix B 1 , the new solution will not require the system to be solved, as ðkÞ
T ðkÞ ¼ B 1 T inc . The heat field disturbance at the receiver xrec can now be computed as Z T k1 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ðkÞ ðx; ωÞ ds:
ð16Þ
S1
S2
At the end of this iteration the total heat field recorded at the receiver will be M
Tðxrec ; ωÞ ¼ T inc ðxrec ; xs ; ωÞ þ ∑ T 0 m ðxrec ; ωÞ; m¼1
ð13Þ
where M is the number of inclusions (in this case M ¼ 2). Iteration k—Step 1: The first crack is reached only by the field generated by the second crack (in the conditions defined by Step 2 of iteration k 1) (see Fig. 4). At this stage the incident heat field corresponds to the field generated by the second inclusion in the previous iteration and is given by Z ðk 1Þ T 21 ðx0 ; nn2 ; ωÞ ¼ H ðx; nn1 ; nn2 ; x0 ; ωÞ T ðk 1Þ ðx; ωÞ ds; ð14Þ S2
Iteration k—Step 2: The second crack is reached only by the heat field generated by the first crack (subjected to the incident heat field in Step 1) (see Fig. 5). The heat field generated by the first crack in the previous step is assumed to be the only heat field which reaches the second R ðkÞ crack, T 12 ðx0 ; nn2 ; ωÞ ¼ S1 Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞds, leading to Z ðkÞ a T ðkÞ ðx0 ; ωÞ ¼ Hðx; nn1 ; nn2 ; x0 ; ωÞ T ðkÞ ðx; ωÞ ds þ T 12 ðx0 ; nn2 ; ωÞ: S2
ð17Þ ), This leads to a system of ½N2 N 2 equations (C T ðkÞ ¼ T ðkÞ inc similar to the one defined before, in Eq. (15), where the constant ðkÞ
matrix needs to be replaced by T ðkÞ inc ðx0 ; nn2 ; ωÞ ¼ T 12 ðx0 ; nn2 ; ωÞ.
14
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
Values of T ðkÞ can be obtained since T ðkÞ ¼ C 1 T ðkÞ . The new heat inc field produced by this inclusion recorded at xrec is expressed by Z ð18Þ T k2 ðxrec ; ωÞ ¼ Hðx; nn1 ; xrec ; ωÞ T ðkÞ ðx; ωÞ ds: S2
At the end of iteration k the total heat field recorded at the receiver will be k
Tðxrec ; ωÞ ¼ T inc ðxrec ; xs ; ωÞ þ ∑
M
∑ T km ðxrec ; ωÞ:
k¼0m¼1
ð19Þ
The iterative process continues until the contribution of the last new field to the total heat field at a specific receiver reaches a predefined threshold. The proposed iterative process requires only the individual linear system of equations related to each inclusion to be solved. Given the example used to illustrate the algorithm procedure, for two inclusions, this requires solving two individual systems of ½N 1 N 1 and ½N 2 N 2 equations only once. The full modelling process would require solving a system of ½ðN1 þ N 2 Þ ðN 1 þ N 2 Þ equations. The proposed iterative approach is more relevant as the number of cracks increases. The size of the system of equations used by the full model becomes greater when compared with the system associated with each inclusion (as used in the proposed iterative model), which raises the interest in using this approach. 2.3. Heat field in the time domain The heat field in the time domain is determined by applying a numerical inverse fast Fourier transform in the frequency domain.
Aliasing phenomena are dealt with by introducing complex frequencies with a small imaginary part, taking the form ωc ¼ ω iη (where η ¼ 0:7Δω and Δω is the frequency step). This shift is subsequently taken into account in the time domain by applying an exponential window (eη t ) to the response. The source can have any time variation and the frequency domain solution can be determined by applying a time Fourier transformation to a frequency range from f¼ 0.0 Hz up to high frequencies. However, we do not need to compute the highest frequencies in the range since the heat response falls rapidly at higher frequencies. The frequency f¼0.0 Hz corresponds to the static response. The use of complex frequencies allows this response to be computed since the argument of the functions in the integral equations is iη, i.e., other than zero.
3. Results and discussion 3.1. Case studies The proposed model is used to simulate heat diffusion around four 3D null thickness cracks, represented by the planes in Fig. 6. The inclusions are embedded in a uniform unbounded solid medium of density ρ¼2300 kg m 3, thermal conductivity λ¼ 1.40 W m 1 1C 1 and specific heat c¼880 J kg 1 1C 1. The point heat source O is located at the origin ðxs ; ys ; zs Þ ¼ ð0 m ; 0 m ; 0 mÞ. Temperature computations are recorded on two fine grids of receivers placed over 0.5 0.4 m2 planes: Grid G1 which is placed away from the source along the z axis, parallel to the xy plane, at z¼0.08 m and grid G2 which is placed parallel to the xz plane at y¼ 0.0 m, as seen in Fig. 6. Each grid contains a total of 4141 receivers spaced at equal intervals of Δx¼0.005 m, Δy¼0.01 m and Δz¼ 0.007 m. As shown in Fig. 7, four 0.2 0.2 m2 cracks are placed parallel to the yz plane, away from the source O along the x axis, 0.1 m apart from each other. Each plane is discretized into 800 elements: 40 along the z axis by 20 along the y direction. 3.2. Simulation results
Fig. 6. 3D view of the systems modelled.
In this section, the performance of the proposed iterative 3D TBEM is studied. This is done by comparing the results obtained for modelling heat diffusion in the presence of 3D cracks using the full 3D normal-derivative equation formulation (3D TBEM) with those given by using the iterative approach. All computations where performed using a 2.30 GHz Intels Core™ i5-2410M CPU with 4.00 GB RAM. The time taken to compute the generated heat field (CPU time) is compared using three distinct approaches: (I) using the full 3D TBEM to calculate heat diffusion in the cracked medium; (II) defining two sub-domains by pairing two cracks and using the iterative 3D TBEM for each subdomain; and (III) using the iterative 3D TBEM and considering each crack separately.
Fig. 7. 2D view of the geometry of the systems modelled: (a) yz plane; (b) xy plane and (c) xz plane.
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
The number of iterations needed to solve the problem using the proposed iterative model was also recorded. The iteration loop is stopped and the number of iterations (NI ) is registered when the temperature at each individual receiver does not change by more than a predefined threshold (ε) when compared to the former
15
iteration, according to the following equation:
NI M
NI M NI 1 M
T ðx ; ωÞ T ðx ; ωÞ ∑ ∑ km rec
= ∑ ∑ T iter
∑ ∑ km rec
k ¼ 0 m ¼ 1
k ¼ 0 m ¼ 1 k¼0m¼1
Fig. 8. Temperature responses, CPU time and number of iterations recorded in G2.
ðx ; ωÞ
r ε: m rec
ð20Þ
16
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
Fig. 9. Snapshots of heat field in time domain taken at (a) t ¼ 5.0 h, (b) t ¼10.6 h and (c) t¼ 12.8 h.
Additionally, transient temperature results are presented in order to illustrate heat diffusion in the time domain simulation in the cracked medium. 3.2.1. Performance of the iterative 3D TBEM Three excitation frequencies have been selected to illustrate the main findings: f ¼0.0 Hz, f ¼0.5 10 6 Hz and f¼ 6.4 10 5 Hz. Fig. 8 shows the real and imaginary temperature responses recorded in grid G2 for each of the approaches, along with the required CPU time and number of iterations. The results clearly indicate that, for any frequency, using an iterative method greatly reduces CPU time, when compared to the full 3D TBEM method. The performance of the iterative method is slightly better when each crack is considered individually (approach III) as CPU time is further reduced. The results show that the number of iterations (and the CPU time) varies for each receiver. Fewer iterations are required in zones which are less affected by disturbances in the heat field, which, for approach III (iterative 3D TBEM considering each crack separately), is the zones behind the cracks placed at x ¼0.2 m and x ¼0.4 m. In the case of the two sub-domains approach (II), it is the zone before the first crack and after the last crack. In both iterative approaches, the number of iterations diminishes as the frequency increases, however CPU time increases slightly. 3.2.2. Heat field in time domain Transient temperature computations are performed in the frequency range ½0:0; 1:024 10 3 Hz with a frequency increment of Δf¼0.5 10 6 Hz, which results in a time window of 1/(0.5 10 6) s. The imaginary part of the frequency is given by η ¼ 0:7 2 π Δf . The host medium is assumed to be initially at 20.0 1C. The source starts emitting energy at instant t¼0.5 h and continues for 9.5 h.
The heat source time dependence is assumed to be rectangular with an amplitude of P¼6.27. This amplitude is defined so that a maximum temperature increase of 15.0 1C is registered by a receiver located in grid G2 at ðx; y; zÞ ¼ ð0:1 m; 0 m; 0 mÞ when the medium is not cracked. In order to illustrate the main results, snapshot images taken at different instants are shown in Fig. 9. Fig. 9a clearly illustrates the effect that the first crack has on temperature diffusion. By instant t¼10.6 h (see Fig. 9b) the heat source has stopped emitting energy and the temperature at ðx; y; zÞ ¼ ð0:1 m; 0 m; 0 mÞ, immediately in front of the center of the first crack, reaches its peak. The snapshot in Fig. 9c shows the instant when the temperature is maximum in the zone immediately before the second crack.
4. Conclusions This paper presents an iterative model to simulate 3D heat diffusion in an unbounded spatially uniform solid medium with embedded 3D inclusions. Heat diffusion is generated by a point heat source and the resulting heat field in the presence of multiple cracks is computed using a normal-derivative integral equation (3D TBEM) formulation. An iterative approach to the 3D TBEM formulation has been described. The method to obtain time-domain responses from frequency-domain calculations has also been described. The performance of the proposed iterative 3D TBEM was assessed by comparing results generated using the full 3D normal-derivative equation formulation (3D TBEM) and those using the iterative model. Two different iterative approaches were considered: solving each crack separately or considering sub-domains by pairing cracks and solving two at each turn. By solving each crack separately, only one inclusion is computed at a time, at each step of the iterative approach. Employing this method required solving a
C. Serra et al. / Engineering Analysis with Boundary Elements 41 (2014) 10–17
lower number of small systems of equations rather than one large system, and led to a significant decrease in CPU time (for both iterative approaches). For a given frequency, the CPU time and the number of iterations required for each receiver varied in both iterative approaches. The number of iterations required fell as the analyzed frequency increased, though the average CPU time increased slightly. Considering each crack separately generated only slightly better results than when sub-domains of pairs of inclusions were considered. It is concluded that using the proposed iterative model greatly reduces the computational effort needed to simulate 3D heat diffusion generated by a point heat source in an unbounded uniform medium with multiple cracks. In order to further demonstrate the applicability of the method, a numerical application to transient temperature calculations has also been presented. Acknowledgments The research work presented herein was supported by the FEDER funds through the Operational Programme for Competitiveness Factors—COMPETE and by national funds through the FCT —Portuguese Foundation for Science and Technology, under research project PTDC/ECM/114189/2009 and was supported in part by the doctoral grant SFRH/BD/91686/2012 attributed by the FCT. This work has been framed under the Initiative Energy for Sustainability of the University of Coimbra and supported by the Energy and Mobility for Sustainable Regions—EMSURE Project (CENTRO-07-0224-FEDER-002004). References [1] Maldague X. Introduction to NDT by active infrared thermography. Mater Eval 2002;6:1060–73. [2] Tartarini P, Corticelli M, Tarozzi L. Dropwise cooling: experimental tests by infrared thermography and numerical simulations. Appl Therm Eng 2009;29:1391–7. [3] Rodríguez F, Nicolau V. Inverse heat transfer approach for IR image reconstruction: application to thermal non-destructive evaluation. Appl Therm Eng 2012;33–34:109–18. [4] Bathe KJ. Numerical methods in finite element analysis. New Jersey: PrenticeHall; 1976. [5] Ozisik MN. Finite difference methods in heat transfer. USA: CRC Press Inc.; 1994. [6] Juncu Gh. Unsteady conjugate forced convection heat/mass transfer from a finite flat plate. Int J Therm Sci 2008;47:972–84.
17
[7] Brebbia CA, Telles JC, Wrobel LC. Boundary elements techniques: theory and applications in engineering. Berlin-New York: Springer-Verlag; 1984. [8] Ochiai Y. Steady heat conduction analysis in orthotropic bodies by triplereciprocity BEM. Comput Model Eng Sci 2001;2:435–46. [9] Šarler B. Towards a mesh-free computation of transport phenomena. Eng Anal Bound Elem 2002;26:731–8. [10] Rudolphi TJ. The use of simple solutions in the regularisation of hypersingular boundary integral equations. Math Comput Model 1991;15:269–78. [11] Cruse TA. Boundary element analysis in computational fracture mechanics. Dordrecht: Kluwer Academic Publishers; 1988. [12] Sládek V, Sládek J. Transient elastodynamics three-dimensional problems in cracked bodies. Appl Math Model 1984;8:2–10. [13] Prosper D. Modeling and detection of delaminations in laminated plates (Ph.D. thesis). Cambridge: MIT; 2001. [14] Prosper D, Kause E. Wave scattering by cracks in laminated media. In: CD: advances in computational engineering and sciences, proceedings of the international conference on computational engineering and science ICES0 01. Atluri SN, Nishioka T, Kikuchi M, (editors), Tech Science Press; 2001. [15] Amado Mendes P, Tadeu A. Wave propagation in the presence of empty cracks in an elastic medium. Comput Mech 2006;38:183–99. [16] Matsumoto T, Zheng C, Harada S, Takahashi T. Explicit evaluation of hypersingular boundary integral equation for 3-D Helmholtz equation discretized with constant triangular element. J Comput Sci Technol 2010;4:194–206. [17] Tadeu A, Prata J, Simões N. Closed form integration of singular and hypersingular integrals in 3D BEM formulations for heat conduction. Math Probl Eng, 2012; Article ID 647038, 21 p, http://dx.doi.org/10.1155/2012/647038. [18] Ma F, Chatterjee J, Henry DP, Banerjee PK. Transient heat conduction analysis of 3D solids with fibre inclusions using the boundary element method. Int J Numer Methods Eng 2008;73:1113–36. [19] Jablonski P. Integral and geometrical means in the analytical evaluation of the BEM integrals for a 3D Laplace equation. Eng Anal Bound Elem 2010;34:264–73. [20] Xianyun Q, Zhang J, Xie G, Zhou F, Li G. A general algorithm for the numerical evaluation of nearly singular integrals on 3D boundary element. J Comput Appl Math 2011;235:4174–86. [21] Škerget L, Rek Z. Boundary-domain integral method using a velocity–vorticity formulation. Eng Anal Bound Elem 1995;15:359–70. [22] Ramšak M, Škerget L, Hriberšk M, Zunič Z. A multidomain boundary element method for unsteady laminar flow using stream function—vorticity equations. Eng Anal Bound Elem 2005;29:1–14. [23] Ramšak M, Škerget L. 3D multidomain BEM for solving the Laplace equation. Eng Anal Bound Elem 2007;31:528–38. [24] Valente FP, Pina HL. Iterative techniques for 3-D boundary element method systems of equations. Eng Anal Bound Elem 2001;25:423–9. [25] Van der Vorst HA, Dekker K. Conjugate gradient type methods and preconditioning. J Comput Appl Math 1988;24:73–87. [26] Marburg S, Schneider S. Performance of iterative solvers for acoustic problems. Part I. Solvers and effect of diagonal preconditioning. Eng Anal Bound Elem 2003;27:727–50. [27] Ylä-Oijala P, Järvenpää S. Iterative solution of high-order boundary element method for acoustic impedance boundary value problems. J Sound Vib 2006;291:824–43. [28] Alia A, Sadok H, Souli M. CMRH method as iterative solver for boundary element acoustic systems. Eng Anal Bound Elem 2012;36:346–50.