Heat conduction analysis by adaptive iterative BEM-FEM coupling procedures

Engineering Analysis with Boundary Elements 73 (2016) 79–94

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Heat conduction analysis by adaptive iterative BEM-FEM coupling procedures D. Soares a,n, L. Godinho b a b

Structural Engineering Department, Federal University of Juiz de Fora, CEP 36036-330 Juiz de Fora, MG, Brazil ISISE, Department of Civil Engineering, University of Coimbra, 3030-788 Coimbra, Portugal

art ic l e i nf o

a b s t r a c t

Article history: Received 9 May 2016 Received in revised form 31 July 2016 Accepted 8 September 2016

This work explores the application of coupled numerical models in thermal conduction analysis, taking into account frequency domain formulations. For this purpose, iterative coupling techniques between the Boundary Element Method (BEM) and the Finite Element Method (FEM) are discussed, also considering adaptive discretization procedures. Two coupling approaches are studied here, both using optimal relaxation parameters to ensure and/or to speed up the convergence of the iterative analysis. Non-matching discretizations on the common interfaces of the different subdomains of the model are allowed, and a single iterative cycle is adopted, incorporating both the adaptive and the coupling processes simultaneously. The methodology is quite flexible, and the possibility of using independent discretizations on each subdomain greatly facilitates the use of enhanced adaptive remeshing; in this context, specialized adaptive techniques for each subdomain/methodology may be applied. The adoption of a single iterative loop also renders a very efficient methodology, avoiding the excessive computational costs of sequentially chained iterative cycles. Numerical examples are presented in the end of the paper, illustrating the good performance of the proposed techniques, as well as their potentialities in engineering applications. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Iterative coupling Adaptive analysis Boundary elements Finite elements Heat conduction Optimal relaxation parameters

1. Introduction The analysis of complex systems may be more effectively handled considering the combination of different numerical methods. In this context, each numerical technique can be applied dealing with the particularities of the model that better fit its positive features. Taking into account combined formulations, a great amount of research works focuses on the coupling of the Boundary Element Method (BEM) and the Finite Element Method (FEM), which are very popular numerical techniques that contribute with complementary beneficial features. In fact, the BEM is quite suitable to handle infinite or semi-infinite media and high variations or discontinuous behaviour, while the FEM is very effectively applied to model complex configurations, in which heterogeneity, anisotropy, nonlinear behaviour etc. may occur. Indeed, it did not take long until some researchers started to combine the BEM and the FEM in order to profit from their respective advantages, trying to evade their disadvantages, and nowadays several works dealing with the BEM-FEM coupling are available [1–20]. Considering transient applications, the first works n

Corresponding author. E-mail addresses: delfi[email protected] (D. Soares), [email protected] (L. Godinho). http://dx.doi.org/10.1016/j.enganabound.2016.09.003 0955-7997/& 2016 Elsevier Ltd. All rights reserved.

dealing with the BEM-FEM coupling were concentrated in the establishment of a coupled system of equations, as reported by manuscripts dealing with time and frequency domain analyses [2– 4]. Later on, iterative coupling algorithms have been proposed, considering once again time [12–15] and frequency [16–19] domain approaches (in this case, an overview is presented by Soares and Godinho [20]). In the iterative coupling approach, each subdomain of the global model is analysed independently (as an uncoupled model) and a successive updating of the variables at the common interfaces is performed, until convergence is achieved. These iterative methodologies exhibit several advantages when compared to standard coupling schemes, as for instance: (i) different sub-domains can be analysed separately, leading to smaller and better-conditioned systems of equations (different solvers, suitable for each sub-domain, may be employed); (ii) only interface routines are required when one wishes to use existing codes to build coupling algorithms (thus, coupled systems may be solved by separate program modules, taking full advantage of specialized features and disciplinary expertise); (iii) matching nodes at common interfaces are not required, greatly improving the flexibility and versatility of the coupled analyses, especially when different discretization methods are considered; (iv) more efficient analyses can be obtained, once the global model can be reduced to several sub-domains with reduced size matrices; etc.

80

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

Iterative coupling methodologies have been extensively applied to analyse elliptic [5–11] and hyperbolic [12–20] models in the last years, and advanced techniques are nowadays available taking into account these configurations. However, iterative approaches have not yet been properly explored considering parabolic models, and the lack of information regarding the performance of iterative coupling procedures applied to these type of problems takes place considering both time and transformed domain analyses. In the present work, adaptive iterative BEM-FEM coupling procedures are proposed, in order to analyse transient heat conduction models. In this case, a frequency domain approach is considered and a single iterative loop is employed to account for the coupling and the adaptive discretization procedures. Adaptive approaches are important to seemly treat complex models, in which solutionbased refinement is necessary to properly represent intricate phenomena. Since in the present coupling methodology a successive updating of computed variables is already settled, adaptive refinement can be introduced into the analyses without further expenditures. As far as the authors are concerned, this is the first work on the topic, and it proposes advanced numerical procedures to analyse a physical phenomenon where coupled BEM-FEM procedures have not yet been properly explored. Transient heat conduction problems can arise in the design and analysis of different structures. Common examples include, for instance, the heat distribution and evolution during the construction of dams or other structures and the simulation of building heat losses and thermal bridges (see, for example, [21]). Closed-form solutions are only known for very simple examples, and thus numerical strategies are frequently necessary to deal with practical and realistic situations. Classical approaches usually tackle this problem directly in the time domain [22–25], making use of timemarching algorithms. Several other works propose addressing the problem in a transformed domain [26–28], mainly using the Laplace transform. However, in that case, difficulties may arise regarding inversion procedures, which may compromise the accuracy of the solution, and special measures need to be adopted. Both approaches have already been successfully applied over the years using the finite element method, the finite difference method and the boundary element method. Mostly in the last decades, a strategy has been proposed and developed trying to avoid the limitations of classical procedures. In this strategy, a Fourier transformation is applied directly to the governing differential equation and the solution is carried out in the frequency domain, similarly to what is regularly considered taking into account hyperbolic models. Later, the time domain solution can be computed by means of the inverse Fourier transformation. This procedure has shown to avoid the accuracy drawbacks of the transformed domain approaches and to be very effective using a variety of numerical methods, including the Boundary Element Method, the Method of Fundamental Solutions, the Kansa’s method and the Meshless Local Petrov Galerkin [29–33]. As it was previously stated, two iterative analyses are considered here, namely the coupling approach and the adaptive refinement. Multiple iterative processes have already been considered taking into account coupled BEM-FEM formulations, mostly considering solid mechanics applications (see, for instance, the work of Elleithy et al. [6], Boumaiza and Aour [10], Soares and Godinho [11] etc.). However, in most of these works, two (or more) iterative loops are considered, being a complete iterative cycle carried out within an iterative step of an external iterative cycle (and so on), which describes a very computationally demanding configuration. The work of Soares and Godinho [11] considering inelastic models is very effective in this sense, illustrating that several iterative analyses may be carried out together, within a unique iterative step, rendering more efficient techniques. Two iterative coupling approaches are discussed here. In these approaches, either temperatures or heat fluxes may be considered prescribed to the BEM or FEM common interfaces. Once these two

approaches are available, their use can be decided in accordance to the characteristics of the model. In fact, it is important to observe that the effectiveness of the proposed coupling formulations is deeply related to the properties of the model, a feature that has been previously observed by Coulier et al. [18], taking into account solid mechanics applications, and that also occurs here. In addition, for both proposed coupling approaches, optimal relaxation parameters are employed to speed up and/or to ensure the convergence of the iterative coupling analysis. As it has been reported [20], frequency domain analyses of hyperbolic models usually give rise to ill-posed problems and, in these cases, the convergence of simple iterative coupling algorithms can either be too slow or unachievable. Thus, in order to deal with illposed problems and ensure convergence of the iterative coupling algorithms, an optimal iterative procedure is adopted here, with optimal relaxation parameters being computed at each iterative step. The paper is organized as follows: first, the governing equations of the transient heat conduction model are briefly presented, and the basic aspects of the BEM and FEM are described; in the sequence, the iterative BEM-FEM coupling algorithms are discussed, and some advanced techniques are further referred, in order to improve the performance of the discussed methodologies; finally, at the end of the paper, numerical applications are considered, illustrating the versatility and good performance of the proposed techniques.

2. Governing equations The basic governing equations related to the transient heat conduction model are given by:

∇⋅(κ (x) ∇τ (t , x) ) − ρ(x) c (x) τ (̇ t , x) = γ (t , x) for x ∈ Ω

(1)

τ (t , x) = τ¯(t , x) for x ∈ Γτ

(2a)

υ(t , x) = ∇τ (t , x)⋅n(x) = υ¯(t , x) for x ∈ Γυ

(2b)

where Eq. (1) and (2) stand for domain ( Ω ) and boundary ( Γ ) equations, respectively, and the boundary of the model is divided into an essential ( Γτ ) and a natural ( Γυ ) boundary, where Γτ ∪ Γυ = Γ and Γτ ∩ Γυ = 0. In these equations, τ and υ stand for the temperature and the heat flux, respectively, and t and x describe the time and the spatial domain of the model. κ , ρ and c stand for the thermal conductivity, density and specific heat, respectively, and γ describes a domain source term. The boundary unit outward normal vector is represented by n, and overdots indicate time derivatives as well as overbars indicate prescribed values. Making use of a Fourier transform, these equations can be rewritten in the frequency domain as follows:

∇⋅(κ (x) ∇τ (ω, x)) − ωi ρ(x)c (x) τ (ω, x) = γ (ω, x) for x ∈ Ω

(3)

τ (ω, x) = τ¯(ω, x) for x ∈ Γτ

(4a)

υ(ω, x) = ∇τ (ω, x)⋅n(x) = υ¯(ω, x) for x ∈ Γυ

(4b)

where i stands for the imaginary unit and ω represents the frequency domain. Since frequency domain analyses are focused here, null initial conditions are considered. In the next subsections, boundary and finite element techniques are briefly presented, describing how the governing Eqs. (3) and (4) can be numerically treated.

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

81

2.1. The boundary element method The boundary element system of equations that describes the focused heat conduction model is given by:

(C + H)τ = Gυ + s

(5)

where τ and υ stand for the temperature and flux vectors, respectively, G and H represent influence matrices, C stands for the location matrix (which is equal to ½I if constant boundary elements are employed, where I is the identity matrix) and vector s accounts for source terms. Matrices G and H can be computed taking into account each boundary element e of the model, and vector s can be evaluated considering each domain integration cell c:

υ N ∫Γ ⌢

H=∪ e

dΓ

(6a)

e

∫Γ ⌢τ N dΓ e

G=∪

(6b)

e

s=∪ c

∫Ω

⌢ τ γ dΩ

Fig. 1. Spatial interpolation procedures: example of an interpolation scheme of vm values in order to obtain v (linear interpolation: v = (djvi + divj)/(dj + di) ).

(6c)

c

In equations (6), ⌢ υ stand for the temperature and the τ and ⌢ flux fundamental solutions, respectively, and N represents the adopted BEM interpolation matrix. Taking into account 2D analysis, the fundamental solutions of the model are given by ⌢ τ (ω, r ) = − aH0(br ) and ⌢ υ (ω, r ) = abH1(br )(∇r⋅n), where a = 1/(4κ ),

b = −iωρc/κ , Hm stands for the Hankel function of the second kind and order m, and r represents the distance between the source and field points. After considering the boundary conditions of the problem (translating all the known variables to the right hand side of Eq. (5), and the unknown fields to the left hand side), the BEM responses of the model can be computed for the given frequency ω. Eq. (5) and (6) only intend to summarily describe the boundary element formulation considered here; for further details on the topic, references [33–35] are suggested.

FEM BEM

L/2

Table 1 Physical properties of the materials (example 1).

The finite element system of equations that represents the heat conduction model is given by:

(7)

where K and M stand for the finite element conductivity and specific heat matrix, respectively, and f represents the excitation vector, accounting for the contribution of domain sources and boundary prescribed fluxes. The terms represented in Eq. (7) may be described as follows, considering each finite element e:

K=∪ e

∫Ω

BT κ B dΩ

(8a)

e

Material

κ

c

ρ

1 (concrete) 2 (steel) 3 (insulator)

1.4 63.9 0.027

880 434 1210

2300 7832 55

temperature field. Once again, Eqs. (7) and (8) only intend to summarily describe the finite element formulation considered here; for further details on the topic, references [36–38] are suggested.

3. Coupled BEM-FEM analysis On the interface between two subdomains, field continuity conditions are defined as F

τ (ω, x) = Bτ (ω, x)

∫Ω e

(9a)

T

M=∪

N ρ c N dΩ

(8b)

e

F

⎛ f=∪⎜ e ⎝

L/2

L

Fig. 2. Sketch of the circular inclusion and of the initial FEM discretization.

2.2. The finite element method

(K + iω M)τ = f

receptor

source

∫Γ

υe

NT υ¯ dΓ −

∫Ω

⎞ NT γ dΩ⎟ ⎠ e

(8c)

where N stands for the adopted FEM interpolation matrix and B stands for the gradient matrix (i.e., B = ∇N ). By re-arranging the system of Eq. (7), taking into account the essential boundary conditions of the problem, the FEM system of equations can be solved, allowing computing the requested

(κ (x)υ(ω, x)) + B(κ (x)υ(ω, x)) = 0

(9b)

where the superscripts B and F indicate if a variable is related to the BEM or to the FEM subdomain, respectively. In the present work, two different coupling procedures are considered, taking into account the boundary conditions at the BEM-FEM common interface: (i) in the first coupling approach (coupling 1), the common interface stands as an essential or Dirichlet boundary for the FEM (i.e., temperatures are prescribed as boundary conditions), and as a natural or Neumann boundary for

82

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94 -3

-3

x 10

1.5

1

1

0.5

0.5

Temperature

Temperature

1.5

0

0

-0.5

-1

x 10

-0.5

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

-1

2.4

0.6

0.8

1

1.2

1.4

1.6

1.8

-5

2

2.2

2.4 -5

x 10

x 10

(a)

(b)

Fig. 3. Temperature at the receptor for the homogeneous model (solid line – analytical solution; circular symbol – real part of the computed results; square symbol – imaginary part of the computed results): (a) coupling 1; (b) coupling 2.

22

Coupling 1 Coupling 2

Iterations

20

the FEM temperatures F τk + 1, at the iterative step kþ1, are evaluated. These FEM temperatures are then considered to compute the fluxes at the BEM nodes at the common interfaces. This is carried out as indicated below: B

υ(ω, B x)k + 1 = [F g(ω, B x)k + 1⋅Bn(B x)](Fκ (B x)/Bκ (B x))

(10)

where B x stands for the BEM node in focus, B n stands for the BEM unit outward normal vector, B κ and F κ stand for the BEM and FEM subdomain conductivities, respectively, and F g stands for the FEM gradient of the temperature, which is computed at the finite element level as function of F τ (i.e., g = Bτ ). In the sequence, these values are applied to the BEM as prescribed boundary conditions, taking into account a relaxation parameter α . Thus, the following vector may be computed:

18

16

B

14 5.0x10

1.0x10

1.5x10

2.0x10

2.5x10

Fig. 4. Total amount of iterative steps for the homogeneous model.

υ¯ k + 1 = (α ) B υ k + 1 + (1 − α ) B υ¯ k

(11)

Once the prescribed variables at the BEM common interfaces are computed, the BEM subdomains can be analysed, allowing to compute the temperatures B τk + 1 at the common interfaces. These values are applied to the FEM as prescribed boundary conditions, following equations (9). Thus, the FEM prescribed values may be obtained considering the following expression:

the BEM (i.e., fluxes are prescribed as boundary conditions); (ii) in the second approach (coupling 2), the inverted relation takes place and the common interface stands as a Neumann boundary for the FEM and as a Dirichlet boundary for the BEM. For both cases, an iterative updating of the variables at the common interface is carried out, taking into account Eq. (9). In the next subsections these two coupling approaches are detailed. In the sequence, some advanced numerical techniques are further discussed, completely describing the proposed adaptive iterative coupling algorithms. It is important to highlight that, in the present work, non-matching BEM and FEM nodes may be considered, which is particularly important when adaptive discretizations are employed. Thus, special procedures to spatially link the FEM and BEM computed values must be taken into account. In this context, mapping and interpolating routines are deeply associated to the algorithms that follow.

where, the FEM temperatures are evaluated taking into account the interpolation of the BEM computed values (in Eq. (12), N stands for interpolating coefficients). These interpolated values can be evaluated taking into account the interpolating functions of the correlated elements or they can be obtained considering enhanced procedures. Fig. 1 illustrates, for instance, a linear interpolation procedure that can be easily applied in order to compute a variable “v” as function of adjacent “vm” variables. Once the FEM prescribed boundary conditions are computed, the iterative cycle is reinitiated, and all the above described procedures are repeated until convergence is achieved.

3.1. Coupling 1

3.2. Coupling 2

In this first algorithm, the FEM subdomain is analysed first, and

F

τ¯(ω, F x)k + 1 =

∑ BN m(F x)Bτ(ω, B x m)k + 1 m

(12)

In this second iterative coupling algorithm, the FEM temperatures

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

83

1.0

Relaxation Parameter

0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 2

4

6

8

10

12

14

16

Iterative Step

(a)

(b)

(c)

(d)

Fig. 5. Computed results for the homogeneous model (coupling 1) considering ω¼ 10  5: (a) final FEM mesh; (b) evolution of the relaxation parameter (circular symbol – real part; square symbol – imaginary part); and computed temperatures at the FEM subdomain – (c) real part; (d) imaginary part. F

τk + 1 are evaluated and applied to the BEM subdomains, taking into account their values at the BEM nodes of the common interfaces, following relations (9). Thus, the temperatures at the BEM nodes of the common interfaces can be computed as follows: B

τ (ω, B x)k + 1 =

∑ FN m(B x)Fτ(ω, F x m)k + 1 m

(13)

where, once again, N stands for interpolating coefficients. In the sequence, these values are employed to compute the BEM prescribed boundary conditions, taking into account a relaxation parameter α , as it is indicated below:

B

τ¯ k + 1 = (α ) B τ k + 1 + (1 − α ) B τ¯ k

(14)

Once the prescribed variables at the BEM common interfaces are evaluated, the BEM subdomains can be analysed, allowing to calculate the fluxes B υ k + 1 at the common interfaces. These values are applied to the FEM as prescribed boundary conditions, following equations (9). Thus, the FEM prescribed values may be obtained considering the following expression: F

⎡ ⎤ υ¯(ω, F x)k + 1 = − ⎢ ∑ BN m(F x)Bυ(ω, B x m)k + 1⎥(Bκ (F x)/ Fκ (F x)) ⎢⎣ m ⎥⎦

(15)

where, once again, the FEM quantities are evaluated taking into

84

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

1.0

Relaxation Parameter

0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 2

4

6

8

10

12

14

16

18

Iterative Step

(a)

(b)

(c)

(d)

Fig. 6. Computed results for the homogeneous model (coupling 2) considering ω¼ 10  5: (a) final FEM mesh; (b) evolution of the relaxation parameter (circular symbol – real part; square symbol – imaginary part); and computed temperatures at the FEM subdomain – (c) real part; (d) imaginary part.

3.3. Adaptive discretization

Table 2 Material distributions in the FEM-BEM coupled models.

Homogeneous Heterogeneous 1 Heterogeneous 2

FEM

BEM

1 (concrete) 2 (steel) 3 (insulator)

1 (concrete) 1 (concrete) 1 (concrete)

account the interpolation of the BEM computed values. Once the FEM prescribed boundary conditions are computed, the iterative cycle is reinitiated, and all the above described procedures are repeated until convergence is achieved.

Taking into account the frameworks of the above discussed iterative coupling algorithms, adaptive discretization techniques can be easily introduced into the analysis, once the proposed algorithms already consider non-matching BEM/FEM nodes. In this context, the discussed algorithms may adapt or not (according to the user’s choice) the FEM meshes, in a step parallel to the BEM subdomains solution, as well as they may adapt or not the BEM meshes, in a step parallel to the FEM subdomains solution. These subdomains adaptive procedures can be carried out taking into

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

x 10

-3

1

0.5

0.5

0

0

-0.5

-0.5

Temperature

Temperature

1

-1 -1.5

-1.5 -2

-2.5

-2.5

-3

-3

0.6

0.8

1

1.2

1.4

1.6

1.8

2

2.2

-3.5

2.4 x 10

-3

-1

-2

-3.5

x 10

0.6

0.8

1

1.2

-4

6

5

4

4

3

2

0

0

1

1.2

1.4

2.2

2.4 -5

1.6

1.8

2

2.2

-1

2.4 x 10

x 10

-4

2

1

0.8

2

3

1

0.6

1.8

(b1)

5

-1

1.6

x 10

Temperature

Temperature

x 10

1.4

-5

(a1) 6

85

0.6

0.8

1

1.2

1.4

1.6

-5

(a2)

1.8

2

2.2

2.4 x 10

-5

(b2)

Fig. 7. Temperature at the receptor for the heterogeneous models 1 and 2 (solid line – analytical solution; circular symbol – real part of the computed results; square symbol – imaginary part of the computed results): (a) coupling 1; (b) coupling 2.

account the computed results of the subdomain, introducing no modifications into the basic format of the discussed coupling algorithms. In addition, following this configuration, just one iterative loop may be considered to deal with both the BEM-FEM coupling and the adaptive discretization (i.e., just a single re-mesh may be considered within an iterative step of the coupling algorithm), resulting in a very efficient approach. In this work, just the FEM meshes are enabled to adapt. This is the case since we believe that this approach provides more efficient analyses for the problems that are usually focused when BEM-FEM coupling procedures are employed; i.e., for problems in which the FEM is employed to discretize an intricate small extension of the model (which is usually where the focus of the analysis lays, demanding good discretizations), and the BEM is employed to discretize the remainder of the model (which is usually of great extension, and the BEM plays the role of a semi-infinite/infinite boundary). In this context, the adaptive refinement of the BEM subdomains does not expressively contribute to enhance the accuracy of the analysis, and it is usually preferable to compute and treat

the BEM matrices just once (i.e., to adopt a non-adaptive approach), improving the efficiency of the solution process. Moreover, creating a boundary mesh is considerably simpler than creating a domain mesh, and elaborating an initial boundary mesh that is already properly refined may not represent a hard task to the user. Thus, the importance of the adaptive process is reduced in the BEM context. To discretize the FEM subdomains, the present work employs linear triangular finite elements, once discretizations considering this type of element are easier to adaptively refine. The adaptive procedure adopted here considers a solve-estimate-mark-refine/ coarsen approach, as described by Chen and Zhang [38]. Following this adaptive methodology, the subdomains discretized by the FEM are analysed and errors are estimated taking into account the computed solutions, at each subdomain; these errors are then employed to mark the elements that will be refined/coarsened, allowing new meshes to be generated maintaining the conformity and shape regularity of the elements. Here, it is considered that the FEM discretization may adapt just at the first η iterative steps of the analysis, where η is provided by the user. Thus, the user can

86

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

B

(a) 80

(17b)

Substituting Eq. (17) into Eq. (16) yields: Coupling 1 Coupling 2

¯ ||2 = ( α 2|| w||2 f (α ) = || α w + (1 − α ) w ¯ ) + (1 − α )2|| w ¯ ||2 ) + 2α(1 − α ) (w, w

60

Iterations

v¯ k = (α ) B v k + (1 − α ) B v¯ k − 1

(18)

where the inner product definition is employed (e.g., (w, w) = || w||2) and new variables, as defined in equations (19), are considered.

40

20

w = Bv k+1 − Bv k

(19a)

¯ = B v¯ k − B v¯ k − 1 w

(19b)

To find the optimal α that minimizes the functional f (α ), Eq. (18) is differentiated with respect to α and the result is set to zero, as described below:

0 5.0x10

1.0x10

1.5x10

2.0x10

2.5x10

¯ ) + (α − 1)|| w ¯ ||2 = 0 α|| w||2 + (1 − 2α ) (w, w

(20)

Re-arranging the terms in Eq. (20), yields:

(b)

¯,w ¯ − w)/|| w ¯ − w||2 α = (w

80

which is an easy to implement expression that provides an optimal value for the relaxation parameter α , at each iterative step. It is important to observe that, by following Eq. (21), the relaxation parameter is computed taking into account just BEM values at the common interfaces. Thus, the computation of the relaxation parameter is not directly influenced by the adaptive discretization of the FEM subdomains.

Coupling 1 Coupling 2 60

Iterations

(21)

40

4. Numerical applications 20

0 5.0x10

1.0x10

1.5x10

2.0x10

2.5x10

Fig. 8. Total amount of iterative steps for the heterogeneous model (a) 1; and (b) 2.

control if no adaptive refinement is to be considered ( η = 0) or if very refined models are to be employed (high η values). 3.4. Optimal relaxation parameter As it has been previously reported [5–20], the effectiveness of iterative coupling techniques is closely related to the introduction of the relaxation parameter, and its value plays a main role in the efficiency of the iterative process. In this context, an unsuitable selection for α may considerably increase the number of iterative steps in the analysis and, in the limit, it may turn convergence unachievable. Thus, employing optimal values for the relaxation parameter is of great importance. In order to evaluate an optimal relaxation parameter, the following square error functional is here minimized:

f (α ) = ||B v¯ k + 1(α ) − B v¯ k(α )||2

(16)

where vector v stands for the temperature or flux vector, according to the coupling algorithm in focus. Taking into account the relaxation of vector v for the k + 1 and k iterations, Eq. (17) may be written, regarding Eqs. (11) or (14): B

v¯ k + 1 = (α ) B v k + 1 + (1 − α ) B v¯ k

(17a)

In this section, two numerical examples are presented, illustrating the performance and potentialities of the discussed techniques. In the first example, a circular inclusion, which is modelled by the FEM, is considered within an infinite domain, which is modelled by the BEM. In this case, different materials are adopted for the inclusion, illustrating the performances of both coupling procedures taking into account the properties of the model. In the second example, a multimaterial wall, which divides a water reservoir, is considered. In this case, the wall is modelled by the FEM, whereas the water subdomains are discretized by the BEM. For both examples, the computed results are compared to analytical answers, whenever possible. For all the analyses that follow, constant boundary elements and linear triangular finite elements are adopted (with η = 10 for the adaptive discretization). For the first step of the iterative analyses, a real value of α = 0.5 is assumed for the relaxation parameter, and Eq. (21) is then followed, providing different complex values for the relaxation parameter, as solution evolves. A tight relative error tolerance of 10  5 is considered for the convergence of the iterative process, which is based on the prescribed BEM values at the common interfaces. In the subsections that follow, SI units are adopted. 4.1. Circular inclusion in an infinite domain This problem consists of a circular inclusion within an infinite domain. A thermal source is applied outside the inclusion, and temperatures are measured at its opposite side, at a receptor. A sketch of the model is depicted in Fig. 2. The initial FEM discretization is also illustrated in the figure, and 40 boundary elements of equal size are employed for the BEM discretization. The geometry of the problem is defined by L = 1. The physical properties of the model are described in Table.1.

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

87

1.0

Relaxation Parameter

0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 2

4

6

8

10

12

14

Iterative Step

(a)

(b)

(c)

(d) 5

Fig. 9. Computed results for the heterogeneous model 1 (coupling 2) considering ω¼ 10 : (a) final FEM mesh; (b) evolution of the relaxation parameter (circular symbol – real part; square symbol – imaginary part); and computed temperatures at the FEM subdomain – (c) real part; (d) imaginary part.

Initially, a homogeneous model is considered and both the inclusion and the infinite domain are characterized by the so-called material 1 (concrete). For this configuration, the computed temperatures at the receptor are depicted in Fig. 3, taking into account different frequencies and coupling algorithms 1 and 2 (analytical answers are also depicted in the figure, for reference). As one can observe, good agreement is obtained between the computed results and the analytical solutions, considering both coupling procedures. However, in this case, as Fig. 3 illustrates, coupling algorithm 2 is able to provide a slightly improved accuracy, taking into account some frequencies. In Fig. 4, the amounts of iterative steps, considering both coupling procedures, are illustrated. As the figure describes, for this homogeneous model, basically the same amount

of iterative steps is required by both coupling procedures, and their efficiencies are equivalent, with coupling algorithm 1 being slightly more efficient. In Fig. 5, the computed temperature field along the FEM subdomain is depicted, taking into account the first coupling algorithm and ω ¼ 10  5. The final computed FEM mesh and the evolution of the relaxation parameter along the iterative analysis are also depicted in the figure. As expected, the FEM mesh is refined in the region closer to the source, in accordance to the distribution of the computed fields (see Fig. 5(a)). Fig. 5(b) illustrates the dynamic behaviour of the computed relaxation parameter, describing its intricate variation along the analysis. This complex evolution is expected, since the evaluation of this parameter is based on

88

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

1.4 1.2

Relaxation Parameter

1.0 0.8 0.6 0.4 0.2 0.0 -0.2 -0.4 2

4

6

8

10

12

Iterative Step

(a)

(b)

(c)

(d)

Fig. 10. Computed results for the heterogeneous model 2 (coupling 1) considering ω ¼10  5: (a) final FEM mesh; (b) evolution of the relaxation parameter (circular symbol – real part; square symbol – imaginary part); and computed temperatures at the FEM subdomain – (c) real part; (d) imaginary part.

residuals, as described by Eqs. (19) and (21). Analogous results are described in Fig. 6, taking into account the second coupling algorithm. As one can observe, the computed temperatures along the FEM subdomain taking into account both coupling algorithms (see Figs. 5(c) and (d) and 6(c) and (d)), are very similar, as well as the obtained final FEM meshes (see Figs. 5(a) and 6(a)). As a second stage of this first example, two heterogeneous models are considered, as described in Table.2. In the so-called heterogeneous model 1, the circular inclusion is made of steel; whereas, in the heterogeneous model 2, the inclusion represents a thermal insulator. In both cases, the surrounding infinite domain is still made of concrete.

Fig. 7 describes the computed temperatures at the receptor, considering these both heterogeneous models, coupling algorithms 1 and 2, and different frequencies. Analytical responses are also depicted in the figure, for reference. Once again, the proposed algorithms provide good results, with slightly improved accuracy being provided by the coupling algorithm 2, when the heterogeneous model 1 is considered. In Fig. 8, the amounts of iterative steps, considering both coupling procedures and both heterogeneous models, are illustrated. In this case (i.e., taking into account the criterion of efficiency), a huge difference of performance is observed between the techniques, depending on the heterogeneous model in focus. As Fig. 8

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

(a)

Table 3 Physical properties of the materials (example 2).

2.0x10

Temperature

1.5x10

Material

κ

c

ρ

1 2 3 4

0.606 1.4 0.72 0.027

4181 880 1860 1210

998 2300 780 55

(water) (concrete) (concrete) (insulator)

1.0x10

Source

5.0x10

0.0 Time

0.0

2.0x10

4.0x10

6.0x10

8.0x10

1.0x10

Time

(b)

Source

4.0x10

3.0x10

Temperature

89

Time

2.0x10

1.0x10

0.0 0.0

2.0x10

4.0x10

6.0x10

8.0x10

1.0x10

Time Fig. 11. Time history of the temperature at the receptor (solid line – analytical solution; dash line – computed results) considering the depicted time behaviour for the source: (a) heterogeneous model 2 and coupling 1; (b) heterogeneous model 1 and coupling 2.

(a) illustrates, for the heterogeneous model 1, the coupling procedure 2 is much more effective, providing good results with a considerably reduced amount of iterative steps. On the other hand, as Fig. 8(b) illustrates, for the heterogeneous model 2, the opposite behaviour occurs. In this case, the coupling procedure 1 is much more effective, providing good results with a considerably reduced amount of iterative steps. In fact, the performance of the coupling

procedures depends on the material proprieties involved in the analysis. As previously remarked taking into account analogous coupling procedures and different physical phenomena [18], the coupling algorithm 1 is supposed to be more effective when the conductivity of the FEM subdomain is considerably lower than that of the BEM subdomain; and the coupling algorithm 2 is supposed to be more effective, otherwise. Alternatively stated, better performance is supposed to be obtained when prescribed temperatures are imposed to the subdomain of considerably lower conductivity. In Fig. 9, the computed temperature field along the FEM subdomain is depicted, taking into account the heterogeneous model 1 and the coupling algorithm 2 (for ω ¼10  5). The final computed FEM mesh and the evolution of the relaxation parameter along the iterative analysis are also depicted in the figure. In Fig. 10, analogous results are presented, considering the heterogeneous model 2 and the coupling algorithm 1. The temperatures in the time domain can be obtained by a numerical inverse fast Fourier transform. In this case, complex frequencies, with a small imaginary part, may be employed, in order to avoid the aliasing phenomenon. Therefore, the definition ωm = (m − iλ )Δω can be adopted, where λ controls the introduction of this imaginary part and Δω stands for the frequency increment. Later, in the time domain, the shift caused by this adoption can be taken into account by applying an exponential window to the response, which is given by the form e λΔωt . Following this procedure, time domain results are depicted in Fig. 11, considering the heterogeneous model 2 and the coupling algorithm 1 and the heterogeneous model 1 and the coupling algorithm 2. The time behaviour of the applied source is also depicted in the figure. In this analysis, 512 complex frequencies are considered, with Δω = 2π⋅10−7 and λ = 0.7. As one can observe in Fig. 11, a very good agreement is observed between the time domain analytical solutions and the computed results. Although very low frequencies are adopted considering the analyses that are carried out here (which is usually the case taking into account heat conduction models [29–33]), it must be highlighted that the proposed iterative coupling algorithms remain stable for higher frequencies. In fact, it has been observed that

source

BEM

Material 1

FEM

Material 2

FEM

source

2H

2H

BEM

Material 1

H

FEM

Material 2

H

Material 3

H

FEM

Material 4

H

FEM

Material 2

H

FEM

Material 2

H

BEM

Material 1

BEM

Material 1

(a)

Material 3

(b)

Fig. 12. Sketch of the multiple layers model (a) without and (b) with the introduction of an insulator (material 4).

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94 0.3

0.3

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

Temperature

Temperature

90

0.05 0

0.05 0

-0.05

-0.05

-0.1

-0.1

-0.15

-0.15

-0.2

1

1.5

2

2.5

3

3.5

4

4.5

-0.2

5 x 10

1

1.5

2

2.5

3

3.5

4

4.5

-6

5 x 10

(a)

-6

(b)

Fig. 13. Temperature at the middle point of the multiple layers model (solid line – analytical solution; circular symbol – real part of the computed results; square symbol – imaginary part of the computed results): (a) coupling 1; (b) coupling 2.

lower amounts of iterative steps are usually obtained when higher frequencies are considered, which is in contrast to what usually occurs considering hyperbolic models [20].

50

Coupling 1 Coupling 2

Iterations

45

4.2. Multi-material wall in a reservoir In this application, a model composed of multiple layers is analysed. A sketch of the model is depicted in Fig. 12. The example describes a multi-material wall, dividing a water reservoir. In this case, the wall is discretized by the FEM (heterogeneous subdomain) and the water is discretized by the BEM (semi-infinite subdomains). The geometry of the problem is defined by H = 0.04 , and the physical properties of the model are described in Table.3. A regular initial mesh is considered for the FEM, with the size of the triangles being characterized by 0.05  0.02 (width x height). For each BEM subdomain, 40 elements of size 0.025 are considered in the central region of the model, which are followed by larger elements, sufficiently extending the discretized domain of analysis. Initially, the configuration described in Fig. 12(a) is considered.

40

35

30 1.0x10

2.0x10

3.0x10

4.0x10

5.0x10

0.3

0.3

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.1

Temperature

Temperature

Fig. 14. Total amount of iterative steps for the multiple layers model.

0.05 0

0.05 0

-0.05

-0.05

-0.1

-0.1

-0.15

-0.15

-0.2 -0.5

-0.4

-0.3

-0.2

-0.1

0 x

(a)

0.1

0.2

0.3

0.4

0.5

-0.2 -0.5

-0.4

-0.3

-0.2

-0.1

0 x

0.1

0.2

0.3

0.4

0.5

(b)

Fig. 15. Temperature at the middle line of the multiple layers model considering ω ¼10  6 (solid line – analytical solution; circular symbol – real part of the computed results; square symbol – imaginary part of the computed results): (a) coupling 1; (b) coupling 2.

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

91

0.2

(a)

0.1 0 -0.1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

(b)

(c)

Fig. 16. Computed results for the multiple layers model (coupling 1) considering ω ¼10  6: (a) final FEM mesh; and computed temperatures at the FEM subdomain – (b) real part; (c) imaginary part (solution with 35 iterative steps).

0.36

Temperature

0.34

0.32

0.30

0.28

0.26

-0.10

-0.05

0.00

0.05

0.10

x Fig. 17. Modulus of the temperature at the middle line of the multiple layers model (zoom at the central region) considering ω¼ 10  6 and coupling 1: solid line – analytical solution; white symbol – adaptive analysis (finishing with 2060 FE); light grey symbol – analysis with fixed 480 FE (initial FE mesh of the adaptive analysis); grey symbol – analysis with fixed 1920 FE; dark grey symbol – analysis with fixed 7680 FE.

In this case, the FEM subdomain is formed by 3 layers, with the outer layers being composed of the so-called material 2, and the inner layer being composed of material 3 (see Table 3). This sequence of layers is employed so that a less conductive core is included to reduce the heat conduction capacity of the wall. The two outer layers of regular concrete correspond to an usual material solution in water reservoirs, which can ensure specific features such as water-tightness. A source is applied at the BEM subdomain, as depicted in Fig. 12. For this configuration, the computed temperatures at the middle point of the wall are depicted in Fig. 13, taking into account different frequencies and coupling algorithms 1 and 2 (analytical responses [39] are also depicted in the figure, for reference). As one can observe, good agreement is obtained between the computed results and the analytical solutions, considering both coupling algorithms. In Fig. 14, the amounts of iterative steps, considering both coupling procedures, are illustrated. As the figure describes, for this model, similar amounts of iterative steps are required by both coupling procedures, with the coupling algorithm 1 being usually a bit more efficient, taking into account the frequencies in focus. In Fig. 15, the accuracies of the proposed techniques are further illustrated, with the temperatures along a horizontal line at the middle height of the wall (central region) being depicted for ω ¼10  6. In Fig. 16, results that are computed along the FEM subdomain are shown, taking into account the coupling algorithm

92

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

0.2

(a)

0.1 0 -0.1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

(b)

(c)

Fig. 18. Computed results for the multiple layers model with an insulator (coupling 1) considering ω ¼10  6: (a) final FEM mesh; and computed temperatures at the FEM subdomain – (b) real part; (c) imaginary part (solution with 35 iterative steps). 10 1.2

0.8

0.1

Error

Relaxation Parameter

1

0.4

0.01 1E-3 1E-4

0.0

1E-5 -0.4 5

10

15

20

25

30

35

1E-6

5

10

Iterative Step

15

20

25

30

35

Iterative Step

(a)

(b) 6

Fig. 19. Computed results for the multiple layers model with an insulator (coupling 1) considering ω¼ 10 : (a) evolution of the relaxation parameter (circular symbol – real part; square symbol – imaginary part); (b) computed relative errors for the convergence of the iterative analysis (a tolerance of 10  5 is adopted).

1 and ω ¼10  6. In this case, the obtained final FEM mesh is depicted in Fig. 16(a), and the obtained temperature fields are illustrated in Fig. 16(b) and (c). As one can observe, the adaptive discretization properly refines the expected critical regions of the

model, allowing a better simulation by the FEM. This is further illustrated in Fig. 17, where results are depicted for the central area of the horizontal line at the middle of the wall, considering the adaptive analysis in focus and analyses taking into account fixed FE

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

93

0.2

(a)

0.1 0 -0.1 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.2

(b)

0.1 0 -0.1 -0.5 0.2

(c)

0.1 0 -0.1 -0.5

0.2

(d)

0.1 0 -0.1 -0.5

Fig. 20. Evolution of the adaptive FEM mesh for the multiple layers model with an insulator (coupling 1) considering ω ¼10

meshes. In this case, regular fixed FE meshes with 480, 1920 and 7680 elements are considered. For the adaptive analysis, the final solution is obtained considering 2016 FE (see Fig. 16(a)). As Fig. 17 describes, the adaptive approach is able to provide accurate results taking into account a considerably lower amount of finite elements. As a second stage of this second example, the model described in Fig. 12(b) is considered. In this case, a thermal insulator is introduced into the wall, as it is indicated in the figure. Taking into account this configuration, Fig. 18 illustrates the computed results along the FEM subdomain, considering coupling algorithm 1. Thus, in Fig. 18(a), the obtained final FEM mesh is depicted and, in Fig. 18(b) and (c), the obtained temperature fields are described, for ω ¼ 10  6. It is interesting to observe the deep effect that the insulator has over the model: in this case, the adaptive refinement mainly occurs taking into account its contour; and, as expected, it blocks the influence of the thermal source, allowing lower temperatures to remain at its opposite side. In Fig. 19, the obtained evolution of the relaxation parameter is depicted, as well as the computed evolution of the relative errors, which describes the convergence of the iterative analysis, is provided. As the figure illustrates, a slightly oscillatory tendency is observed considering the convergence curve. This is expected, since the proposed evaluation of the relaxation parameter allows its value to freely

6

– iterative steps: (a) 2; (b) 4; (c) 6; (d) 8.

vary, which may produce oscillations in the convergence curve, according to the intensities of these computed variations. In fact, this behaviour can be more clearly observed when Fig. 19(a) and (b) are analysed together, illustrating that the oscillations in both curves are related. Finally, in Fig. 20, the evolution of the adaptive analysis is described, depicting the computed FEM meshes at some iterative steps.

5. Conclusions In this work, two iterative BEM-FEM coupling approaches are discussed to analyse frequency domain heat conduction models. The two approaches are defined according to the prescribed boundary conditions at the BEM-FEM common interfaces and, in both cases, the use of iterative coupling schemes allows the problem to be easily analysed considering independent discretizations on each subdomain of the problem. The use of non-matching nodes at the common interfaces of different subdomains renders the whole process quite flexible and of general application, especially considering complex engineering problems. In addition, it allows an easy implementation of adaptive mesh refinement within the iterative coupling process, which is of great practical interest. Both

94

D. Soares, L. Godinho / Engineering Analysis with Boundary Elements 73 (2016) 79–94

the iterative coupling and the adaptive analysis are performed within the same loop, thus avoiding unnecessary multiple nested iterative loops to perform the coupling and the adaptive refinement. It should be noted that, in addition, in order to allow a faster convergence, an optimal relaxation parameter is computed at each step, also improving the performance of the proposed algorithms. The examples presented here show that the two proposed iterativeadaptive procedures exhibit good convergence properties when they are properly applied (i.e., when the use of each methodology takes into account the features of the model), and that they are highly practical and efficient from a computational point of view. These applications also indicate the applicability of these algorithms to adequately tackle complex and realistic problems. It can be stated that the proposed algorithms present several practical advantages when compared with more traditional approaches (such as direct coupling). The various subdomains are analysed separately, leading to smaller and better-conditioned systems of equations, allowing each of them to be handled with maximal efficiency, using specific strategies tailored for each individual method. Completely independent numerical methods can be used for each part of the model, allowing the best method to be chosen for each subdomain. Here, this has been demonstrated using the BEM and the FEM, the first handling infinite or semi-infinite parts of the problem, and the latter dealing with localized and more intricate structures. Furthermore, the use of iterative coupling allows the use of independent codes for each subdomain, making it easy to use existing routines and specific expertise. Non-matching nodes can also be used along the interface, without any specific limitation, leading to a very flexible and general methodology of analysis, and permitting mesh modification on one subdomain without interfering with the other. This greatly facilitates incorporating adaptive mesh refinement algorithms, which may even act within the same iterative loop of the coupling, as demonstrated in the present paper. Overall, the defined procedures can be seen as very powerful tools for the analysis of thermal conduction problems, allowing the proper handling of complex models.

[8] [9]

[10] [11] [12] [13] [14] [15]

[16]

[17]

[18]

[19]

[20] [21] [22]

[23] [24] [25] [26]

[27]

Acknowledgements The financial support by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais) is greatly acknowledged. This work was also financed by FEDER funds through the Competitivity Factors Operational Programme – COMPETE and by national funds through FCT – Foundation for Science and Technology within the scope of the project POCI-01–0145-FEDER-007633.

[28] [29]

[30]

[31]

[32]

References [1] Zienkiewicz OC, Kelly DM, Bettes P. The coupling of the finite element method and boundary solution procedures. Int J Numer Methods Eng 1977;11:355–76. [2] Pavlatos GD, Beskos DE. Dynamic elastoplastic analysis by BEM/FEM. Eng Anal Bound Elem 1994;14:51–63. [3] von Estorff O, Firuziaan M. Coupled BEM/FEM approach for nonlinear soil/ structure interaction. Eng Anal Bound Elem 2000;24:715–25. [4] Amini S, Harris PJ, Wilton DT. Coupled boundary and finite element methods for the solution of the dynamic fluid-structure interaction problem. Berlin: Springer-Verlag; 1992. [5] Lin CC, Lawton EC, Caliendo JA, Anderson LR. An iterative finite element – boundary element algorithm. Comput Struct 1996;39:899–909. [6] Elleithy WM, Tanaka M, Guzik A. Interface relaxation FEM–BEM coupling method for elasto-plastic analysis. Eng Anal Bound Elem 2004;28:849–57. [7] Elleithy WM, Grzhibovskis R. An adaptive domain decomposition coupled

[33]

[34] [35] [36] [37] [38] [39]

finite element–boundary element method for solving problems in elastoplasticity. Int J Numer Methods Eng 2009;79:1019–104. Elleithy WM. Multi-region adaptive finite element-boundary element method for elasto-plastic analysis. Int J Comput Math 2012;89:1525–39. Jahromi HZ, Izzuddin BA, Zdravkovic L. A domain decomposition approach for coupled modelling of nonlinear soil-structure interaction. Comput Methods Appl Mech Eng 2009;198:2738–49. Boumaiza D, Aour B. On the efficiency of the iterative coupling FEM–BEM for solving the elasto-plastic problems. Eng Struct 2014;72:12–25. Soares D, Godinho L. Inelastic 2D analysis by adaptive iterative BEM–FEM coupling procedures. Comput Struct 2015;156:134–48. Soares D, von Estorff O, Mansur WJ. Iterative coupling of BEM and FEM for nonlinear dynamic analyses. Comput Mech 2004;34:67–73. Soares D. An optimised FEM–BEM time-domain iterative coupling algorithm for dynamic analyses. Comput Struct 2008;86:1839–44. Soares D. Fluid-structure interaction analysis by optimised boundary element - finite element coupling procedures. J Sound Vib 2009;322:184–95. Soares D. FEM-BEM iterative coupling procedures to analyze interacting wave propagation models: fluid-fluid, solid-solid and fluid-solid analyses. Couple Syst Mech 2012;1:19–37. Bendali A, Boubendir Y, Fares M. A FETI-like domain decomposition method for coupling finite elements and boundary elements in large-size problems of acoustic scattering. Comput Struct 2007;85:526–35. Soares D, Godinho L. An optimized BEM-FEM iterative coupling algorithm for acoustic-elastodynamic interaction analyses in the frequency domain. Comput Struct 2012;106–107:68–80. Coulier P, François S, Lombaert G, Degrande G. Coupled finite element – hierarchical boundary element methods for dynamic soil–structure interaction in the frequency domain. Int J Numer Methods Eng 2014;97:505–30. Godinho L, Soares D. Frequency domain analysis of interacting acoustic-elastodynamic models taking into account optimized iterative coupling of different numerical methods. Eng Anal Bound Elem 2013;37:1074–88. Soares D, Godinho L. An overview of recent advances in the iterative analysis of coupled models for wave propagation. J Appl Math 2014:426283. Branco FG, Tadeu A, Simões N. Heat conduction across double brick walls via BEM. Build Environ 2004;39:51–8. Chang YP, Kang CS, Chen DJ. The use of fundamental Green functions for solution of problems of heat conduction in anisotropic media. Int J Heat Mass Transf 1973;16:1905–18. Shaw RP. An integral equation approach to diffusion. Int J Heat Mass Transf 1974;17:693–9. Wrobel LC, Brebbia CA. A formulation of the boundary element method for axisymmetric transient heat conduction. Int J Heat Mass Transf 1981;24:843–50. Dargush GF, Banerjee PK. Application of the boundary element method to transient heat conduction. Int J Numer Methods Eng 1991;31:1231–47. Cheng AHD, Abousleiman Y, Badmus T. A Laplace transform BEM for axysymmetric diffusion utilizing pre-tabulated Green’s function. Eng Anal Bound Elem 1992;9:39–46. Zhu SP, Satravaha P, Lu X. Solving the linear diffusion equations with the dual reciprocity methods in Laplace space. Eng Anal Bound Elem 1994;13:1–10. Zhu SP, Satravaha P. An efficient computational method for nonlinear transient heat conduction problems. Appl Math Model 1996;20:513–22. Tadeu A, Simões N, Simões I. Coupling BEM/TBEM and MFS for the simulation of transient conduction heat transfer. Int J Numer Methods Eng 2010;84 (2):179–213. Tadeu A, António J, Godinho L, Simoes N. Boundary element method analyses of transient heat conduction in an unbounded solid layer containing inclusions. Comput Mech 2004;34(2):99–110. Godinho L, Dias-da-Costa. The MLPG (5) for the Analysis of Transient Heat Transfer in the Frequency Domain. CMES: Computer Modeling in Engineering and Sciences, 96(5); 2013, p. 293–316. Godinho L, Branco F. Formulation of Kansa's method in the frequency domain for the analysis of transient heat conduction. Int J Numer Methods Heat Fluid Flow 2014;24(7):1437–53. Godinho L, Tadeu A, Simões N. Study of transient heat conduction in 2.5D domains using the boundary element method. Eng Anal Bound Elem 2004;28:593–606. Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Berlin: Springer-Verlag; 1984. Wrobel LC. The boundary element method: applications in thermo-fluids and acoustics. England: Wiley; 2002. Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method, its basis and fundamentals. 6th ed.. Oxford: Butterworth-Heinemann; 2005. Lewis RW, Nithiarasu P, Seetharamu KN. Fundamentals of the Finite element method for heat and fluid flow. England: Wiley; 2004. Chen L, Zhang CS. AFEM@matlab: a Matlab package of adaptive finite element methods. University of Maryland at College Park; 2006. Simões N, Tadeu A. FUndamental solutions for transient heat transfer by conduction and convection in an unbounded, half-space, slab and layered media in the frequency domain. Eng Anal Bound Elem 2005;29:1130–42.