# An Overview of Recent Advances in the Iterative Analysis of Coupled Models for Wave Propagation

Authors: D. Soares; L. Godinho
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Wave propagation problems can be solved using a variety of methods. However, in many cases, the joint use of different numerical procedures to model different parts of the problem may be advisable and strategies to perform the coupling between them must be developed. Many works have been published on this subject, addressing the case of electromagnetic, acoustic, or elastic waves and making use of different strategies to perform this coupling. Both direct and iterative approaches can be used, and they may exhibit specific advantages and disadvantages. This work focuses on the use of iterative coupling schemes for the analysis of wave propagation problems, presenting an overview of the application of iterative procedures to perform the coupling between different methods. Both frequency- and time-domain analyses are addressed, and problems involving acoustic, mechanical, and electromagnetic wave propagation problems are illustrated.


---

## Body

## 1. Introduction

The analysis of wave propagation, either involving electromagnetic, acoustic, or elastic waves, has been widely studied by researchers using different strategies and methodologies, as can be seen, for example, in [1–10], among many others. In many cases, the interaction between different types of media, such as fluid-solid or soil-structure interaction problems, poses significant challenges that can hardly be tackled by means of a single numerical method, requiring the joint use of different procedures to model different parts of the problem. Indeed, taking into consideration the specificities and particular features of distinct numerical methods, their combined use, as coupled or hybrid models, has been proposed by many authors, in order to explore the individual advantages of each technique.In acoustic and elastodynamic problems, coupled models, including, for example, the joint use of the boundary element method (BEM) and the method of fundamental solutions (MFS) [11] or of the BEM and the meshless Kansa’s method [12], have been successfully applied. Similarly, when modelling dynamic fluid-structure and soil-structure interactions, wave propagation in elastic media with heterogeneities, or the transmission of ground-borne vibration, coupled models using the finite element method (FEM) and the BEM have been extensively documented in the literature [13–19], mostly using the FEM to model the structure and the BEM to model the hosting infinite or semiinfinite medium. Although these approaches can be quite useful in addressing many engineering problems, they mostly correspond to standard direct coupling methodologies and thus exhibit well-known limitations. Indeed, directly coupling distinct methods involves assembling a single system matrix, accounting for the contributions of each method and for the required coupling interface conditions, which frequently becomes poorly conditioned due to the different nature of the methods. Since this system is formed from the contributions of distinct methods, it is also usually not possible to make use of their individual advantages in terms of optimized solvers or memory storage (e.g., in BEM-FEM the final system will no longer be banded and symmetric, etc.). In addition to this limitation, by forming a single system of equations, a very large problem usually arises, leading to increased computational efforts and thus to a loss of performance.All these limitations have justified the appearance of iterative algorithms to obtain accurate solutions in a more efficient manner. Perhaps one of the first iterative techniques to be developed for general problems is the well-known Schwarz alternating strategy [20, 21], in which the domain of analysis is partitioned in overlapping subdomains, and the solution is found by successively iterating along these subdomains until convergence is reached. This classical and simple to implement algorithm has been applied to many problem types, including potential problems [22] or electromagnetic wave propagation problems [23]. However, for the case of acoustic problems or elastic wave propagation problems, formulated in the frequency domain, the special oscillatory structure of the solution leads to severe convergence problems when using such classic approaches, and more sophisticated and difficult to implement strategies must be defined.In recent years, more elaborate iterative domain decomposition techniques have been proposed and discussed in order to analyze a wide range of problems, providing good results especially in terms of flexibility and efficiency. Mostly, these techniques have been applied to nontransient applications, and they usually consider the analysis of coupled models, taking into account the interaction of different discretization methods, physical phenomena, and so forth. In fact, for complex models, iterative domain decomposition techniques are recommended, usually providing a better approach for the analysis. Indeed, a proper numerical simulation is hardly achieved by a single numerical technique in those cases, mostly because complex and quite different phenomena interact, requiring particularized advanced expertise, and/or large scale problems are involved, demanding high computational efforts.Nowadays, several works are available discussing iterative nonoverlapping partitioned analysis. Taking into account elliptic problems, Rice et al. [24] presented a quite complete discussion, considering several interface relaxation procedures and comparing formulations and performances. As a matter of fact, most of the publications on the topic are focused on elliptic models, few being devoted to hyperbolic problems. Taking into account computational mechanics, one of the first publications on the topic was presented by Lin et al. [25], which discussed a relaxed iterative procedure to couple the FEM and the BEM, considering linear static analyses. Similar approaches have been presented later on, considering potential and mechanical static linear analyses [26, 27]. In the works of Elleithy et al. [28, 29], concerning mechanical static and potential problems, the authors propose that the domain of the original problem is subdivided into subdomains, each of them modeled by the finite element or boundary element methods; the coupling between the different subdomains is performed using smoothing operators on the interdomain boundaries. Their strategy allows separate computations for the BEM and FEM subdomains, with successive update of the boundary conditions at the interfaces being performed until convergence is achieved. In [30–32], similar approaches for the analysis of different linear problems using domain decomposition techniques were also presented. Further developments of these strategies to nonlinear analysis in solid mechanics can also be found in the works of Elleithy et al. [33], using an interface relaxation finite element-boundary element coupling method (FEM-BEM coupling) for elasto-plastic analysis, or Jahromi et al. [34], who established a coupling procedure based on a sequential iterative Dirichlet-Neumann coupling algorithm for nonlinear soil-structure interaction. It must be noted that the described works refer to nontransient problems, either linear or nonlinear, and no application to wave propagation analysis is focused on in these works.Taking into account time-domain wave propagation models, the first work on the topic seems to have been presented by Soares et al. [35], who described a relaxed FEM-BEM iterative coupling procedure to analyze dynamic nonlinear problems, considering different time discretizations within each sub-domain of the model. Later on, this technique has been further developed to analyze other wave propagation models, including acoustic, elastic, and electromagnetic wave propagation or solid-fluid interaction, taking into account several different numerical procedures using the FEM and the BEM [36–45] or the meshless local Petrov-Galerkin method [46]. Most of these works are focused on the iterative coupling of different numerical discretization techniques, and a review considering the iterative coupling of the FEM and the BEM, taking into account some wave propagation models in computational mechanics, has been presented in [47]. The coupling of acoustic and mechanic wave propagation models, on the other hand, has been reviewed in [48], taking into account different domain decomposition techniques and considering several numerical discretization techniques.In the analysis of wave propagation using frequency-domain formulations, iterative coupling procedures can be found in the literature, mostly considering acoustic-acoustic and acoustic-elastodynamic coupling [49–54]. As it has been reported, frequency-domain wave propagation analyses usually give rise to ill-posed problems and, in these cases, the convergence of the iterative coupling algorithm can be either too slow or unachievable. This is the case in acoustic-acoustic, acoustic-elastodynamic, and elastodynamic-elastodynamic interacting models and, as discussed in this work, convergence can be hardly achieved if no special procedure is considered, especially if higher frequencies are focused on. As referred in the literature, in order to deal with this ill-posed problem and ensure convergence of the iterative coupling algorithm, special techniques, such as the adoption of optimal relaxation parameters, must be considered.In this work, time- and frequency-domain analyses of wave propagation models are reviewed, taking into account relaxed iterative coupling procedures. In this context, several wave propagation models (such as electromagnetic, acoustic, mechanic) are considered, and several numerical procedures (such as the finite element method, the boundary element method, and meshless methods) are employed to discretize the model. In the iterative coupling approach, each sub-domain of the global model is analyzed independently (as an uncoupled model) and a successive renewal 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, for instance,(i)
different subdomains 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)
matching time steps at common interfaces are not required (in time-domain analysis), allowing optimal temporal discretizations within each sub-domain, improving accuracy and stability aspects;
(v)
nonlinear analyses (as well as other iterative-based analyses) may be carried out in the same iterative loop of the iterative coupling, not introducing a relevant extra computational effort for the model;
(vi)
more efficient analyses can be obtained, once the global model can be reduced to several subdomains with reduced size matrices.As a matter of fact, Gauzellino et al. [55] compared the iterative domain decomposition and global solution taking into account three-dimensional Helmholtz problems. Their numerical results show that iterative domain decomposition methods perform far better than global methods. In addition, they observed that iterative domain decomposition methods involving small subdomains work better than those with subdomains involving a large number of elements. Similar results have been obtained by Soares et al. [51], taking into account two-dimensional Helmholtz problems.To give a detailed overview of the recent developments found in many of the referred works, the remainder of this paper will address a number of application examples concerning different phenomena and methods. First, the governing equations related to wave propagation models are generically and briefly presented. In the sequence, an efficient iterative coupling technique is described, including the mathematical derivation of the optimized relaxation methodology. Some numerical applications are finally presented, illustrating the accuracy, performance, and potentialities of the discussed procedures, taking into account different wave propagation models and discretization techniques.


## 2. Governing Equations

Wave propagation phenomena may be generically described by the following time/frequency-domain governing equations:(1a)
c
0
(
x
,
t
)
+
c
1
(
x
)
u
¨
(
x
,
t
)
+
c
2
(
x
)
u
˙
(
x
,
t
)
+
c
3
(
x
)
∂
f
(
u
(
x
,
t
)
)
=
0
,
(1b)
c
-
0
(
x
,
ω
)
-
ω
2
c
1
(
x
)
u
-
(
x
,
ω
)
+
i
ω
c
2
(
x
)
u
-
(
x
,
ω
)
+
c
3
(
x
)
∂
f
(
u
-
(
x
,
ω
)
)
=
0which can be further generalized in order to consider more complex behavior, such as time varying coefficients (c
l
(
x
,
t
), l
=
1,2
,
3), nonlinearities (c
l
(
x
,
u
(
x
,
t
)
), l
=
0,1
,
2,3; etc.) Equation (1a) stands for the time domain governing equation, whereas (1b) stands for its frequency-domain counterpart (overbars indicate frequency-domain values). In these equations, u represents the incognita field, which can be scalar, vectorial, and so forth, according to the physical model in focus. c
i stands for a general coefficient representation, which can as well be a scalar, a tensor, and so forth. Overdots stand for time derivatives, whereas ∂
f indicates a spatial derivative operator. The complex number is denoted by i and the time, frequency, and space domains are represented by t, ω, and x, respectively (in this case, x
∈
Ω, where Ω is the spatial domain of the model).The boundary conditions (x
∈
Γ, where Γ is the boundary of the model) may be generically described as (for simplicity, from this point onwards, overbars are no longer used to indicate frequency-domain values and ς stands for t or ω, according to the case of analysis)
(2)
f
(
u
(
x
,
ς
)
,
v
(
x
,
ς
)
)
=
c
(
x
,
ς
)
,
where, once again, c stands for known terms. In (2), f stands for a generic function, representing the combination of its arguments. The variable v, which may be considered prescribed at the boundary of the model, is a function of u, and it is usually expressed considering some normal projection (normal to the boundary) of the spatial derivatives of u (i.e., v
=
∂
f
n
(
u
)).To completely define the model, initial conditions (which are usually adopted null in frequency-domain analyses) must also be defined. In this case, a generic representation can be given byf
(
u
(
x
,
t
=
0
)
,
u
˙
(
x
,
t
=
0
)
)
=
c
(
x
), where notation analogous to that of (2) is considered.Taking into account coupled models in which different domains interact by a common interface, interface conditions must be stated, indicating how the domains interact. This can be generically expressed as(3)
f
1
(
u
(
x
-
,
ς
)
,
v
(
x
-
,
ς
)
)
=
f
2
(
u
(
x
+
,
ς
)
,
v
(
x
+
,
ς
)
)
,
where x
∈
Γ
I, x
-
∈
Γ
I
∪
Ω
1, x
+
∈
Γ
I
∪
Ω
2, and Γ
I is the common interface between domains Ω
1 and Ω
2. In (3), functions f
1 and f
2 describe how the interaction between the coupled domains takes place by relating their boundary values on the common interface.


## 3. Iterative Coupling Analysis

In order to enable the coupling between sectioned domains of a global model, an iterative procedure is employed here, which performs a successive renewal of the relevant variables at the common interfaces. This approach is based on the imposition of prescribed boundary conditions, properly evaluated, at the interfaces of the sectioned domains, allowing each domain of the global model to be analyzed separately. Since the sectioned domains are analyzed separately, the relevant systems of equations are formed independently, before the iterative process starts (in the case of linear analyses), and are kept constant along the iterative process, rendering a very efficient procedure. The separate treatment of the sectioned domains allows independent discretizations to be considered on each domain, without any special requirement of matching nodes along the common interfaces. Moreover, in the case of time-domain analysis, different time-steps may also be considered for each domain. Thus, the coupling algorithm can be presented for a generic case, in which the interface nodes may not match, and the interface time instants are disconnected, allowing exploiting the benefits of the iterative coupling formulation.To ensure and/or to speed up convergence, a relaxation parameterλ is introduced in the iterative coupling algorithm. The effectiveness of the iterative process is strongly related to the selection of this relaxation parameter, since an inappropriate selection for λ can significantly increase the number of iterations in the analysis or, even worse, make convergence unfeasible. As it has been reported [49, 51], frequency-domain analyses 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. In order to deal with ill-posed problems and ensure convergence of the iterative coupling algorithm, an optimal iterative procedure is adopted here, with optimal relaxation parameters being computed at each iterative step. As it is illustrated in the next section, the introduction of these optimal relaxation parameters allows the iterative coupling technique to be very effective, especially in the frequency domain, ensuring convergence at a low number of iterative steps.


### 3.1. Iterative Algorithm

Initially, in thekth iterative step of the coupled analysis of domains 1 and 2, the so-called domain 1 is analyzed and the variables u or v at the common interfaces of the domain are computed, taking into account prescribed values of v  or u at these common interfaces. These prescribed values of v or u are provided from the previous iterative step (in the first iterative step, null or previous time-step values may be considered). Once the variables u or v are computed, they are applied to evaluate the boundary conditions that are prescribed at the common interfaces of domain 2, as described by (3). Taking into account these prescribed u or v boundary conditions, the so-called domain 2 is analyzed and the variables v  or  u at the common interfaces of the domain are computed. Then, the computed v or u values are applied to evaluate the boundary conditions that are prescribed at the common interfaces of domain 1, reinitiating the iterative cycle. A sketch of this cycle is depicted in Figure 1.Figure 1
Sketch of the iterative coupling algorithm.As previously discussed, relaxation parameters must be considered in order to ensure and/or to speed up the convergence of the iterative process. Thus, the values that are computed after the analysis of the sectioned domain may be combined with its previous iterative step counterpart, relaxing the computation of the actual iterative step value. Mathematically, this can be represented as follows:(4)
y
(
k
+
1
)
=
(
λ
)
y
(
k
+
λ
)
+
(
1
-
λ
)
y
(
k
)
,
where λ is the adopted relaxation parameter and y stands for u or v, according to the case of analysis; one should note that y
(
k
+
λ
) is the value computed at the end of the iterative step, before the application of the relaxation parameter.A proper selection forλ at each iterative step is extremely important for the effectiveness of the iterative coupling procedure. In order to obtain an easy to implement, efficient, and effective expression for the relaxation parameter computation, optimal λ values are deduced in Section 3.2.


### 3.2. Optimal Relaxation Parameter

In order to evaluate an optimal relaxation parameter, the following square error functional is minimized here:(5)
ε
(
λ
)
=
∥
Y
(
k
+
1
)
(
λ
)
-
Y
(
k
)
(
λ
)
∥
2
,
where Y stands for a vector whose entries are u or v values, computed at the common interfaces.Taking into account the relaxation of the field values for the (k
+
1) and (k) iterations, (6a) and (6b) may be written, based on the definition in (4):(6a)
Y
(
k
+
1
)
=
(
λ
)
Y
(
k
+
λ
)
+
(
1
-
λ
)
Y
(
k
)
,
(6b)
Y
(
k
)
=
(
λ
)
Y
(
k
+
λ
-
1
)
+
(
1
-
λ
)
Y
(
k
-
1
)
.Substituting (6a) and  (6b) into (5) yields
(7)
ε
(
λ
)
=
∥
(
λ
)
W
(
k
+
λ
)
+
(
1
-
λ
)
W
(
k
)
∥
2
=
(
λ
2
)
∥
W
(
k
+
λ
)
∥
2
+
2
λ
(
1
-
λ
)
(
W
(
k
+
λ
)
,
W
(
k
)
)
+
(
1
-
λ
)
2
∥
W
(
k
)
∥
2
,
where the inner product definition is employed (e.g., (
W
,
W
)
=
∥
W
∥
2) and new variables, as defined in the following, are considered:
(8)
W
(
k
+
λ
)
=
Y
(
k
+
λ
)
-
Y
(
k
+
λ
-
1
)To find the optimalλ that minimizes the functional ε
(
λ
), (7) is differentiated with respect to λ and the result is set to zero, described as follows:
(9)
(
λ
)
∥
W
(
k
+
λ
)
∥
2
+
(
1
-
2
λ
)
(
W
(
k
+
λ
)
,
W
(
k
)
)
+
(
λ
-
1
)
∥
W
(
k
)
∥
2
=
0
.Rearranging the terms in (9) yields
(10)
λ
=
(
W
(
k
)
,
W
(
k
)
-
W
(
k
+
λ
)
)
∥
W
(
k
)
-
W
(
k
+
λ
)
∥
2
which is an easy to implement expression that provides an optimal value for the relaxation parameter λ, at each iterative step. This expression requires a low computational cost, when compared to other alternatives that can be found in the literature (see, e.g., [28, 29]) and it provides very good results, as it has been reported taking into account different physical models and domain analyses [43, 44, 51–54]. The iterative process is relatively insensitive to the value of the relaxation parameter adopted for the first iterative step and λ
=
0.5 can be considered in this case, for instance.


### 3.3. Interface Compatibility

As previously discussed, independent spatial (and temporal, in time-domain analysis) discretizations may be considered for each domain of the model, not requiring matching nodes (or equal time steps) at the common interfaces. Thus, special procedures must be employed to ensure the interface spatial (and temporal) compatibility. In order to do so, interpolation and extrapolation procedures are considered here. These procedures can be generically described by(11a)
y
(
x
i
,
ς
)
=
∑
j
=
1
J
α
j
y
(
x
j
,
ς
)
,
(11b)
y
(
x
,
t
n
)
=
β
0
y
(
x
,
t
m
)
+
∑
j
=
1
J
β
j
y
(
x
,
(
t
-
j
Δ
t
)
m
/
n
)
,where (11a) stands for spatial interpolations and (11b) stands for time interpolations/extrapolations (α
j and β
j stand for spatial interpolation coefficients and time interpolation/extrapolation coefficients, respectively, where Δ
t represents the time step). In Figure 2, simple sketches for the spatial and temporal interpolation procedures are depicted, taking into account linear interpolations.(a) Sketch for a spatial interpolation of nodal values on the interface:y
1
+
=
I
(
y
1
-
,
y
2
-
), y
2
+
=
I
(
y
1
-
,
y
2
-
), y
3
+
=
I
(
y
2
-
,
y
3
-
), and so forth; (b) sketch for a temporal interpolation of time-step values on the interface: y
+
(
t
+
)
=
I
(
y
-
(
t
-
)
,
y
-
(
t
-
-
Δ
t
-
)
), and so forth, where I stands for a linear interpolation function.
(a)
(b)Although time interpolations usually can be carried out without further difficulties, time extrapolations may give rise to instabilities if not properly elaborated. Thus, extrapolations should be performed in consonance with the field approximations being adopted within each time step and with the time discretization procedures being considered in the analysis, in order to formulate a consistent procedure. Once a consistent methodology is elaborated, time interpolation/extrapolation procedures can be employed with confidence, as referred in the literature [47, 48] and illustrated in the next section. One should notice that usually different optimal (optimal in terms of accuracy, stability and efficiency) time steps are required when taking into account different numerical methods, spatial discretizations, material properties, physical phenomena, and so forth. Thus, in some cases, considering different time steps within each domain of a coupled model is of maximal importance to allow the effectiveness of the analysis.Using space(/time) interpolation(/extrapolation) procedures, optimal modeling of each sectioned domain may be achieved, which is very important in what concerns flexibility, efficiency, accuracy, and stability aspects.


## 3.1. Iterative Algorithm

Initially, in thekth iterative step of the coupled analysis of domains 1 and 2, the so-called domain 1 is analyzed and the variables u or v at the common interfaces of the domain are computed, taking into account prescribed values of v  or u at these common interfaces. These prescribed values of v or u are provided from the previous iterative step (in the first iterative step, null or previous time-step values may be considered). Once the variables u or v are computed, they are applied to evaluate the boundary conditions that are prescribed at the common interfaces of domain 2, as described by (3). Taking into account these prescribed u or v boundary conditions, the so-called domain 2 is analyzed and the variables v  or  u at the common interfaces of the domain are computed. Then, the computed v or u values are applied to evaluate the boundary conditions that are prescribed at the common interfaces of domain 1, reinitiating the iterative cycle. A sketch of this cycle is depicted in Figure 1.Figure 1
Sketch of the iterative coupling algorithm.As previously discussed, relaxation parameters must be considered in order to ensure and/or to speed up the convergence of the iterative process. Thus, the values that are computed after the analysis of the sectioned domain may be combined with its previous iterative step counterpart, relaxing the computation of the actual iterative step value. Mathematically, this can be represented as follows:(4)
y
(
k
+
1
)
=
(
λ
)
y
(
k
+
λ
)
+
(
1
-
λ
)
y
(
k
)
,
where λ is the adopted relaxation parameter and y stands for u or v, according to the case of analysis; one should note that y
(
k
+
λ
) is the value computed at the end of the iterative step, before the application of the relaxation parameter.A proper selection forλ at each iterative step is extremely important for the effectiveness of the iterative coupling procedure. In order to obtain an easy to implement, efficient, and effective expression for the relaxation parameter computation, optimal λ values are deduced in Section 3.2.


## 3.2. Optimal Relaxation Parameter

In order to evaluate an optimal relaxation parameter, the following square error functional is minimized here:(5)
ε
(
λ
)
=
∥
Y
(
k
+
1
)
(
λ
)
-
Y
(
k
)
(
λ
)
∥
2
,
where Y stands for a vector whose entries are u or v values, computed at the common interfaces.Taking into account the relaxation of the field values for the (k
+
1) and (k) iterations, (6a) and (6b) may be written, based on the definition in (4):(6a)
Y
(
k
+
1
)
=
(
λ
)
Y
(
k
+
λ
)
+
(
1
-
λ
)
Y
(
k
)
,
(6b)
Y
(
k
)
=
(
λ
)
Y
(
k
+
λ
-
1
)
+
(
1
-
λ
)
Y
(
k
-
1
)
.Substituting (6a) and  (6b) into (5) yields
(7)
ε
(
λ
)
=
∥
(
λ
)
W
(
k
+
λ
)
+
(
1
-
λ
)
W
(
k
)
∥
2
=
(
λ
2
)
∥
W
(
k
+
λ
)
∥
2
+
2
λ
(
1
-
λ
)
(
W
(
k
+
λ
)
,
W
(
k
)
)
+
(
1
-
λ
)
2
∥
W
(
k
)
∥
2
,
where the inner product definition is employed (e.g., (
W
,
W
)
=
∥
W
∥
2) and new variables, as defined in the following, are considered:
(8)
W
(
k
+
λ
)
=
Y
(
k
+
λ
)
-
Y
(
k
+
λ
-
1
)To find the optimalλ that minimizes the functional ε
(
λ
), (7) is differentiated with respect to λ and the result is set to zero, described as follows:
(9)
(
λ
)
∥
W
(
k
+
λ
)
∥
2
+
(
1
-
2
λ
)
(
W
(
k
+
λ
)
,
W
(
k
)
)
+
(
λ
-
1
)
∥
W
(
k
)
∥
2
=
0
.Rearranging the terms in (9) yields
(10)
λ
=
(
W
(
k
)
,
W
(
k
)
-
W
(
k
+
λ
)
)
∥
W
(
k
)
-
W
(
k
+
λ
)
∥
2
which is an easy to implement expression that provides an optimal value for the relaxation parameter λ, at each iterative step. This expression requires a low computational cost, when compared to other alternatives that can be found in the literature (see, e.g., [28, 29]) and it provides very good results, as it has been reported taking into account different physical models and domain analyses [43, 44, 51–54]. The iterative process is relatively insensitive to the value of the relaxation parameter adopted for the first iterative step and λ
=
0.5 can be considered in this case, for instance.


## 3.3. Interface Compatibility

As previously discussed, independent spatial (and temporal, in time-domain analysis) discretizations may be considered for each domain of the model, not requiring matching nodes (or equal time steps) at the common interfaces. Thus, special procedures must be employed to ensure the interface spatial (and temporal) compatibility. In order to do so, interpolation and extrapolation procedures are considered here. These procedures can be generically described by(11a)
y
(
x
i
,
ς
)
=
∑
j
=
1
J
α
j
y
(
x
j
,
ς
)
,
(11b)
y
(
x
,
t
n
)
=
β
0
y
(
x
,
t
m
)
+
∑
j
=
1
J
β
j
y
(
x
,
(
t
-
j
Δ
t
)
m
/
n
)
,where (11a) stands for spatial interpolations and (11b) stands for time interpolations/extrapolations (α
j and β
j stand for spatial interpolation coefficients and time interpolation/extrapolation coefficients, respectively, where Δ
t represents the time step). In Figure 2, simple sketches for the spatial and temporal interpolation procedures are depicted, taking into account linear interpolations.(a) Sketch for a spatial interpolation of nodal values on the interface:y
1
+
=
I
(
y
1
-
,
y
2
-
), y
2
+
=
I
(
y
1
-
,
y
2
-
), y
3
+
=
I
(
y
2
-
,
y
3
-
), and so forth; (b) sketch for a temporal interpolation of time-step values on the interface: y
+
(
t
+
)
=
I
(
y
-
(
t
-
)
,
y
-
(
t
-
-
Δ
t
-
)
), and so forth, where I stands for a linear interpolation function.
(a)
(b)Although time interpolations usually can be carried out without further difficulties, time extrapolations may give rise to instabilities if not properly elaborated. Thus, extrapolations should be performed in consonance with the field approximations being adopted within each time step and with the time discretization procedures being considered in the analysis, in order to formulate a consistent procedure. Once a consistent methodology is elaborated, time interpolation/extrapolation procedures can be employed with confidence, as referred in the literature [47, 48] and illustrated in the next section. One should notice that usually different optimal (optimal in terms of accuracy, stability and efficiency) time steps are required when taking into account different numerical methods, spatial discretizations, material properties, physical phenomena, and so forth. Thus, in some cases, considering different time steps within each domain of a coupled model is of maximal importance to allow the effectiveness of the analysis.Using space(/time) interpolation(/extrapolation) procedures, optimal modeling of each sectioned domain may be achieved, which is very important in what concerns flexibility, efficiency, accuracy, and stability aspects.


## 4. Numerical Applications

In this section, the general procedures previously discussed are particularized and briefly detailed, taking into account different physical models and discretization techniques. Thus, the discussed iterative coupling methodology is applied considering a wide range of wave propagation models and numerical methods, richly illustrating its performance and potentialities.In this context, time- and frequency-domain analyses are carried out here, and electromagnetic, acoustic, and mechanical wave propagation phenomena (as well as their interactions) are discussed in the applications that follow. Moreover, different numerical techniques (such as the finite element method, the boundary element method, and meshless methods) are applied to discretize the different domains of the model, illustrating the versatility and generality of the discussed iterative method.


### 4.1. Electromagnetic Waves

In electromagnetic models, vectorial wave equations describe the electric and the magnetic field evolution [56, 57]. In this case, (1a) can be rewritten as (in this subsection, time-domain analyses are focused on):(12a)
∇
×
(
μ
(
x
)
-
1
∇
×
E
(
x
,
t
)
)
+
ε
(
x
)
E
¨
(
x
,
t
)
=
-
J
˙
(
x
,
t
)
,
(12b)
∇
×
(
ε
(
x
)
-
1
∇
×
H
(
x
,
t
)
)
+
μ
(
x
)
H
¨
(
x
,
t
)
=
∇
×
(
ε
(
x
)
-
1
J
(
x
,
t
)
)
,and (3) can be rewritten as(13a)
n
(
x
)
×
(
E
(
x
+
,
t
)
-
E
(
x
-
,
t
)
)
=
0
,
(13b)
(
D
(
x
+
,
t
)
-
D
(
x
-
,
t
)
)
·
n
(
x
)
=
ρ
(
x
,
t
)
,
(13c)
(
B
(
x
+
,
t
)
-
B
(
x
-
,
t
)
)
·
n
(
x
)
=
0
,
(13d)
n
(
x
)
×
(
H
(
x
+
,
t
)
-
H
(
x
-
,
t
)
)
=
J
(
x
,
t
)
,where E and H are the electric and magnetic field intensity vectors, respectively; D and B represent the electric and magnetic flux densities, respectively; and J and ρ stand for the electric current and electric charge density, respectively. The parameters ε and μ denote, respectively, the permittivity and permeability of the medium and its wave propagation velocity is specified as c
=
(
ε
μ
)
-
1
/
2. n is the normal vector, from domain 1 to domain 2. Equations (13a) and (13b) state that the tangential component of E is continuous across the interface and that the normal component of D has a step of surface charge on the interface surface, respectively. Equations (13c) and (13d) state that the normal component of B is continuous across the interface and that the tangential component of H is continuous across the interface if there is no surface current present, respectively.In the present application, the electromagnetic fields surrounding infinitely long wires are studied [41]. Two cases of analysis are focused here, namely, (a) case 1, where one wire is considered; (b) case 2, where two wires are employed. For both cases, the wires are carrying time-dependent currents (i.e., I
(
t
)
=
t or I
(
t
)
=
t
2) and they are located along the adopted z-axis. A sketch of the model is depicted in Figure 3.Sketch of the electromagnetic models and adopted FEM/BEM spatial discretizations: (a) case 1, one wire; (b) case 2, two wires.
(a)
(b)The spatial and temporal evolution of the electric field intensity vector is analyzed here taking into account a finite element method (FEM)—boundary element method (BEM) coupled formulation. In this context, the FEM is applied to model the region close to the wires, whereas the BEM simulates the remaining infinity domain. As it is well known, the BEM employs fundamental solutions which fulfill the radiation condition. Thus, this formulation is very suitable to perform infinite domain analysis, once reflected waves from infinity are avoided [58].The adopted spatial discretization is also described in Figure3. In this case, 2344 linear triangular finite elements and 80 linear boundary elements are employed in the analyses (see references [57, 58] for more details regarding the FEM and the BEM applied to electromagnetic analyses). The radius of the FEM-BEM interface is defined by R
=
1
m and matching nodes are considered at the interface. For temporal discretization, the selected time step is given by Δ
t
=
5
·
1
0
-
11
s for both domains. The physical properties of the medium (air) are μ
=
1.2566
·
1
0
-
6
H
/
m and ε
=
8.8544
·
1
0
-
12
F
/
m.Figure4 shows the modulus of the electric field intensity obtained at points A and B (see Figure 3) considering the iterative coupling methodology. Analytical time histories [58] are also depicted in Figure 4, highlighting the good accuracy of the numerical results. In Figure 5, charts are displayed, indicating the percentage of occurrence of different relaxation parameter values (evaluated according to expression (10)), in each analysis. As can be observed, for all considered cases, optimal relaxation parameters are mostly in the interval 0.7
≤
λ
≤
0.8. In fact, an optimal relaxation parameter selection is extremely case dependent. It is function of the physical properties of the model, geometric aspects, adopted spatial and temporal discretizations, and so forth. Equation (10) provides a simple expression to evaluate this complex parameter.Time history results for the electric field intensity at points A and B consideringI
(
t
)
=
t and (a) case 1 and (b) case 2; I
(
t
)
=
t
2 and (c) case 1 and (d) case 2.
(a)
(b)
(c)
(d)Percentage of occurrence of different relaxation parameter values during the analysis, consideringI
(
t
)
=
t and (a) case 1 and (b) case 2; I
(
t
)
=
t
2 and (c) case 1 and (d) case 2.
(a)
(b)
(c)
(d)In order to illustrate the effectiveness of the methodology when considering different time discretizations for different domains, Figure6 depicts results that are computed considering Δ
t
=
2.5
·
1
0
-
11
s for the FEM and Δ
t
=
2.0
·
1
0
-
10
s for the BEM (i.e., a difference of 8 times between the time steps). For simplicity, results are presented considering just the first case of analysis, that is, case 1 and I
(
t
)
=
t. As one can observe in Figure 6(a), good results are still obtained taking into account the iterative formulation, in spite of the existing time disconnections at the interface. In Figure 6(b), the evolution of the relaxation parameter is depicted, taking into account this last configuration. As one can observe, in this case, optimal relaxation parameter values are between 0.7 and 1.0 and mostly concentrate on the interval (
0.9,1.0
). In fact, it is expected that these values get closer to 1.0 when smaller time steps are considered. In the present analysis, an average number of 4.92 iterations per time step is obtained (taking into account 800 FEM time steps), which is a relatively low number, illustrating the good performance of the technique (it must be remarked that a tight tolerance criterion was adopted for the convergence of the iterative analysis).Results considering different time steps for each domain: (a) electric field intensity at points A and B; (b) optimal relaxation parameters for each iterative step.
(a)
(b)


### 4.2. Acoustic Waves

In acoustic models, a scalar wave equation describes the acoustic pressure field evolution [1]. In this case, (1b) can be rewritten as (in this subsection, frequency-domain analyses are focused)
(14)
∇
·
(
κ
(
x
)
∇
p
(
x
,
ω
)
)
+
ω
2
ρ
(
x
)
p
(
x
,
ω
)
=
γ
(
x
,
ω
)
and (3) can be rewritten as(15a)
(
p
(
x
+
,
ω
)
-
p
(
x
-
,
ω
)
)
=
0
,
(15b)
(
q
(
x
+
,
ω
)
-
q
(
x
-
,
ω
)
)
=
g
(
x
,
ω
)
,where p is the hydrodynamic pressure and γ and g stand for domain and surface sources, respectively. The parameters ρ and κ denote, respectively, the mass density and compressibility of the medium and its wave propagation velocity is specified as υ
=
(
κ
/
ρ
)
1
/
2. The hydrodynamic fluxes on the interfaces are represented by q, and they are defined by q
=
κ
∇
p
·
n, where n is the normal vector, from domain 1 to domain 2. Equation (15a) states that the pressure is continuous across the interface, whereas (15b) states that the flux is continuous across the interface if there is no surface source.The advantages of using iterative coupling procedures are revealed when more complex configurations are analyzed. In this subsection, the case of a heterogeneous domain, composed of a homogeneous fluid incorporating multiple circular inclusions with different properties, is analyzed. For this purpose, consider the host medium to allow the propagation of sound with a velocity of 1500 m/s, and this medium is excited by a line source located atx
s
=
-
5
.
0 m and y
s
=
0.0 m. Within this fluid, consider the presence of 8 circular inclusions; all of them are with unit radius and filled with a different fluid, allowing sound waves to travel at 3000 m/s, as depicted in Figure 7.(a) Sketch for the heterogeneous medium with multiple subregions; (b) boundary and domain point distribution, considering the spatial discretization of an inclusion and adjacent fluid.
(a)
(b)The above-described system has been analyzed taking into account the proposed iterative coupling procedure making use of the Kansa’s method (KM) to model all the inclusions and of the method of fundamental solutions (MFS) to model the host fluid (see references [12, 59–61] for more details regarding the KM and the MFS applied to acoustic analyses). One should note that, since the real source is positioned at the outer region, the iterative process is initialized with the analysis of the MFS model, considering prescribed Neumann boundary conditions at the common interfaces. Once the boundary pressures for the outer region are computed, these values are transferred to the closed regions by imposing Dirichlet boundary conditions, incorporating information about the influence of each inclusion on the remaining heterogeneities. Then, each KM subregion is analysed independently and the internal boundary values (normal fluxes) are evaluated autonomously for each inclusion. The iterative procedure then goes further, including the calculation of the relaxation parameter at each iterative step, as well as the correction of boundary variables, until convergence is achieved.To model the system, each MFS boundary is discretized by 55 points. 331 KM domain points are equally distributed within each inclusion, and 66 KM boundary points (around 31 points per wavelength) are used (see Figure7(b)). The complexity of the model hinders the definition of a closed form solution; thus, the results are checked against a numerical model which performs the direct (i.e., noniterative) coupling between both methods. In that model, 66 boundary points are used in the MFS to define the boundary of each inclusion, and 66 and 331 KM boundary and domain points, respectively, are adopted for the discretization of each inclusion (analogously to the iterative coupling procedure). Figure 8 compares the responses computed by the iterative and the direct coupling methodologies. Results are depicted along the boundary of the 8th inclusion, for excitation frequencies of 400 Hz (Figure 8(a)) and 1000 Hz (Figure 8(b)). As can be observed in the figure, there is a perfect match between both approaches, with the iterative procedure clearly converging to the correct solution.Hydrodynamic pressures along the common interface of the 8th inclusion for (a)ω
=
400 Hz and (b) ω
=
1000 Hz (—, real-iterative; - - - -, Imag-iterative; ○, Real-direct; □, Imag-direct).
(a)
(b)It is important, at this point, to highlight the differences in the computational times of the direct and of the iterative coupling approaches. For the present model configuration, the direct coupling approach had to deal simultaneously with 528 boundary points and a total of 2648 internal points (i.e., considering a coupled matrix of dimension 3704), which is implied in 373.89 s of CPU time in a Matlab implementation (being this CPU time independent of the frequency in focus). For the iterative coupling approach, using 55 boundary points for the MFS and 66 boundary points for the KM, it was possible to obtain analogous results considering 12.06 s of CPU time for the frequency of 200 Hz and 32.32 s for the frequency of 1000 Hz (i.e., 3.23% and 8.64% of the computational cost of the direct coupling methodology, resp.). Even if the same number of boundary points is used in the iterative coupling approach for the MFS (i.e., 66 points), the final CPU time would just increase up to 35.71 s (9.55% of the computational cost of the direct coupling methodology). These results are summarized in Table1, where the number of iterations and the CPU time are presented for the first scenarios (i.e., 55 boundary points for the MFS) and for frequencies between 50 Hz and 1000 Hz. The values described in Table 1 further confirm that the difference in calculation times between the iterative and the direct coupling approach is striking and reveal an excellent gain in performance favouring the iterative coupling technique. It is important to understand that this gain is strongly related to the possibility of dealing with smaller-sized matrices when using the iterative coupling procedure. Moreover, it is possible to invert (or triangularize, etc.) the relevant matrices only at the first iterative step and then proceed with the calculations using the inverted matrices (or forward/back substituting, etc.). As a consequence, after the first iteration, only matrix vector multiplication operations are required, and very high savings in what concerns computational time are achieved. In fact, considering that the number of operations required for matrix inversion can be assumed to be of the order of N
3 (N being the matrix size), a simple calculation allows concluding that, for the current model, the relative cost of inverting the eight KM matrices (each one being a square matrix with 397 × 397 entries) and the MFS matrix (with 440 × 440 entries) is less than 2% of the cost of inverting a larger 3704 × 3704 matrix, as required for the direct coupling strategy. Similar conclusions can be obtained considering other solver procedures, such as matrix triangularizations, demonstrating that a considerably less expensive methodology is obtained if the different subdomains are analysed separately (even considering an eventual high number of iterative steps in the iterative analysis).Table 1
Total number of iterations and relative CPU time (iterative/direct coupling) for the acoustic model.
Frequency (Hz)
Iterations
Relative CPU time (%)
50
14
2.09
100
18
2.43
150
29
3.19
200
31
3.23
250
26
2.78
300
40
3.86
350
31
3.13
400
32
3.20
450
30
3.02
500
49
4.49
550
48
4.40
600
45
4.18
650
68
5.94
700
34
3.28
750
39
3.66
800
110
9.20
850
98
8.28
900
78
6.75
950
104
8.83
1000
100
8.64Analyzing the difference in computational times between the two analyzed frequencies (i.e., 200 Hz and 1000 Hz) also reveals a significant difference between them. This difference is related to the number of iterations required for convergence, which was higher when the excitation frequency of 1000 Hz was considered. The plot in Figure9 indicates the number of iterations required for convergence along a range of frequencies between 10 Hz and 1000 Hz, using a constant number of boundary (55 for the MFS and 66 for the KM) and internal points (331 for the KM). As expected, the number of iterations increases with the frequency. It is interesting to note that the maximum necessary number of iterations occurred for a frequency of 990 Hz, requiring 170 iterations and a CPU time of 54.80 s to converge, which is less than 15% of the CPU time required by the direct coupling for the same frequency.Figure 9
Total number of iterations considering different frequencies and 55 collocation points for the MFS and 66 boundary points for the KM, at each circular inclusion.In Figure10, the wavefield produced within and around the inclusions is illustrated for excitation frequencies of 600 Hz and 1000 Hz. As expected, as the frequency increases, the multiple inclusions generate progressively more complex wave fields, with the interaction between them becoming very significant for the higher frequency. Observation of these results also reveals a strong shadow effect produced by the inclusions, with much lower amplitudes being registered in the region behind the inclusions placed further away from the source. This effect is even more pronounced for the higher frequency. Interestingly, for both frequencies, the space between the two lines of inclusions works as a guiding path, along which the sound energy travels with less attenuation.3D plots of the sound field for frequencies of (a) 600 Hz and (b) 1000 Hz.
(a)
(b)


### 4.3. Mechanical Waves

In dynamic models, a vectorial wave equation describes the displacement field evolution [1]. In this case, considering linear behaviour, (1a) can be rewritten as (in this subsection, time-domain analyses are focused)
(16)
∇
×
(
μ
(
x
)
∇
×
u
(
x
,
t
)
)
-
∇
(
(
η
(
x
)
+
2
μ
(
x
)
)
∇
·
u
(
x
,
t
)
)
+
ρ
(
x
)
u
¨
(
x
,
t
)
=
f
(
x
,
t
)
,
and (3) can be rewritten as:(17a)
(
u
(
x
+
,
t
)
-
u
(
x
-
,
t
)
)
=
0
,
(17b)
(
σ
(
x
+
,
t
)
-
σ
(
x
-
,
t
)
)
n
(
x
)
=
τ
(
x
,
t
)
,where u is the displacement vector and f and τ stand for domain and surface forces, respectively. The terms η and μ denote the so-called Lamé parameters, and ρ is the mass density of the medium. In this case, the wave propagation velocities are specified as c
s
=
(
μ
/
ρ
)
1
/
2 (shear wave) and c
d
=
(
(
η
+
2
μ
)
/
ρ
)
1
/
2 (dilatational wave). The stress tensor is denoted by σ and n is the normal vector, from domain 1 to domain 2. Equation (17a) states that the displacements are continuous across the interface, whereas (17b) states that the tractions are continuous across the interface if there are no surface forces on it.One main advantage of the discussed coupling algorithm is that other iterative processes can be carried out in the same iterative loop needed for the coupling. Thus, consideration of coupled nonlinear models, as for example, may not demand a superior amount of computational effort, which is very beneficial.In the present application, a nonlinear model is considered and elastoplastic analyses are carried out (for details about elastoplastic analyses, one is referred to [33–35, 62–64]). Moreover, two discretization approaches are employed here, one taking into account FEM-BEM coupling procedures, and another considering BEM-BEM coupled techniques (for more details about these coupled models, one is referred to [37, 43]). In this context, a nonlinear infinity domain is analyzed here, in which a circular cavity is loaded. The region expected to develop plastic strains is discretized by the finite element method, in the case of the FEM-BEM coupled analysis, or by the domain boundary element method (D-BEM), in the case of the BEM-BEM coupled analysis. The remainder of the infinity domain is discretized by the time-domain boundary element method (TD-BEM). A sketch of the model is depicted in Figure 11, as well as the adopted discretizations.(a) Sketch of the circular cavity; (b) FEM-BEM discretization; (c) BEM-BEM discretization.
(a)
(b)
(c)The FEM-BEM discretization is depicted in Figure11(b). In this case, 1944 linear triangular finite elements and 80 linear boundary elements are employed in the coupled analysis. The BEM-BEM discretization is depicted in Figure 11(c). In this case, 46 linear boundary elements are employed in the BEM-BEM coupled analysis (20 linear boundary elements for the TD-BEM and 26 linear boundary elements for the D-BEM), as well as 270 linear triangular cells (D-BEM formulation). In the BEM-BEM coupled analysis, the double symmetry of the problem is taken into account. An interesting feature of the boundary element formulation is that symmetric bodies under symmetric loads can be analysed without discretization of the symmetry axes. This can be accomplished by an automatic condensation process, which integrates over reflected elements and performs the assemblage of the final matrices in reduced size [64]. The time discretization adopted is given by Δ
t
=
0.04
s for the FEM and Δ
t
=
0.20
s for the D-BEM and the TD-BEM.The physical properties of the model areμ
=
2.652
·
1
0
8
N
/
m
2, η
=
2.274
·
1
0
8
N
/
m
2, and ρ
=
1.804
·
1
0
3
kg
/
m
3. A perfectly plastic material obeying the Mohr-Coulomb yield criterion is assumed, where c
o
=
4.8263
·
1
0
6
N
/
m
2 (cohesion) and ϕ
=
3
0
° (internal friction angle). The geometry of the problem is defined by R
=
3.048
m (the radius of the TD-BEM circular mesh is given by 5R).In Figure12, the displacement time history at point A is depicted, considering linear and nonlinear analyses. As one can notice, good agreement is observed between the FEM-BEM and BEM-BEM results. It is important to highlight that, for the FEM-BEM analyses, a difference of 5 times between the FEM and BEM time steps is considered, illustrating the effectiveness of the time interpolation/extrapolation procedures adopted in the analyses.Figure 12
Radial displacement considering coupled TD-BEM/D-BEM and coupled TD-BEM/FEM analyses: linear and nonlinear results at point A.The number of iterations per time step and the optimal relaxation parameters, evaluated at each iterative step, are depicted in Figure13, taking into account the FEM-BEM coupled analyses. As one may observe, basically the same computational effort (i.e., number of iterative steps) is necessary for both linear and nonlinear analyses, highlighting the efficiency of the proposed methodology for complex phenomena modeling. It is also important to remark the low number of iterative steps necessary for convergence, with a maximum of 7 iterations being necessary, within a time step, taking into account the entire linear and nonlinear analyses. For the focused configurations, the optimal relaxation parameters are intricately distributed within the interval (
0.75
;
1.00
), as depicted in Figure 13(b).TD-BEM/FEM analyses: (a) number of iterations per time step considering optimal relaxation parameters; (b) optimal relaxation parameters for each iterative step.
(a)
(b)In Table2, the total number of iterations is presented, considering analyses with optimal relaxation parameters and with some constant preselected λ values. As one may observe, an inappropriate selection for the relaxation parameter can considerably increase the associated computational effort. Thus, the optimization technique is extremely important in order to provide a robust and efficient iterative coupling formulation. In Figure 14, the computed σ
x
y stresses are depicted, considering the BEM-BEM elastoplastic analysis. An advantage of the D-BEM is that it employs nodal stress equations [37], allowing computing continuous stress fields, in counterpart to the FEM, which computes stresses based on displacement derivatives, obtaining discontinuous stress fields at element interfaces.Table 2
Total number of iterations (considering all time steps) for the dynamic model.
Relaxation parameter
Elastic analysis
Elastoplastic analysis
1.00
3730
3740
0.90
3392
3443
0.80
3973
3993
0.70
4772
4777
Optimal
3287
3346Figure 14
Spatial and temporal evolution ofσ
x
y considering elastoplastic analysis (scale factor: 6.8947
·
10
6
N
/
m
2); (a) t
=
4 s; (b) t
=
8 s; (c) t
=
12 s; (d) t
=
16 s.


### 4.4. Coupled Acoustic-Mechanical Waves

In this case, different wave equations, as indicated in Sections4.2 and 4.3 (see (14) and (16)), describe different domains of the global model. The interface conditions for the acoustic-dynamic coupling (3) can then be written as (in this subsection, frequency-domain analyses are focused)(18a)
(
n
(
x
)
·
σ
(
x
+
,
ω
)
n
(
x
)
-
p
(
x
-
,
ω
)
)
=
0
,
(18b)
(
-
ω
2
n
(
x
)
·
u
(
x
+
,
ω
)
-
υ
(
x
-
)
2
q
(
x
-
,
ω
)
)
=
0
,where u is the displacement vector and σ is the stress tensor of the dynamic model (domain 1). p is the hydrodynamic pressure and q is the hydrodynamic flux of the acoustic model (domain 2). n is the normal vector, from domain 1 to domain 2. The acoustic wave propagation velocity is denoted by υ. Equation (18a) states that the normal components of the dynamic tractions are equal to the acoustic pressures and (18b) relates the normal components of the dynamic accelerations to the acoustic fluxes.In the present application, a model in which a concrete wall of 10.0 m high is coupled to a fluid waveguide, filled with water, is analyzed. A sketch of the model is depicted in Figure15(a). For this case, a pressure source is positioned in the waveguide, at (−10.0; 0.5), illuminating the system. The concrete structure corresponds to a wall with variable cross-section, exhibiting thicknesses of 4.0 m at its basis and of 2.0 m at its top.(a) Sketch of the coupled acoustic-dynamic model; (b) MFS collocation points (o) and virtual sources (x); (c) FEM mesh when 20 nodes are used along the solid-fluid interface; (d) node distribution for the meshless methods when 20 nodes are used along the solid-fluid interface.
(a)
(b)
(c)
(d)To simulate this coupled system, several approaches are employed. For the fluid medium, the MFS is used in all cases, allowing the use of the Green’s function for a waveguide (see [61] for details concerning this function). This Green’s function is written as a summation of modes, and its convergence is very difficult when the source and the receiver are positioned along the same vertical line, thus posing severe difficulties for its use together with a BEM formulation. The structure is modelled using three different methods, namely, the FEM, a local collocation method (CM) (see [65] for details about this procedure), and a meshless local Petrov-Galerking technique (MLPG) (see [65, 66] for details about this procedure). Representations of the node distribution for each one of the methods can be found in Figures 15(b)–15(d).Since no analytical solution can be found in the literature for the present case, a numerical solution making use of a full BEM model ensuring the correct coupling between the solid and the fluid is used. For this case, the rigid bottom of the waveguide is accounted for using an image-source Green’s function, while the free surface is fully discretized up to a distance of 60.0 m from the concrete wall; after this, an anechoic termination is considered, imposing adequate Robin boundary conditions. For the solid, full-space Green’s functions are adopted, and 60 nodal points are used along the solid-fluid interface, ensuring that accurate results can be obtained. A total of 796 boundary elements are used to build this model. Details on the mathematical formulation of this technique can be found in the works of Tadeu and Godinho  [7].Figure16(a) illustrates the reference response in terms of real and imaginary components of the acoustic pressure at a receiver located at (−1.0; 3.0), for frequencies between 1 Hz and 150 Hz. This position is chosen so that the effect of the vibration of the concrete wall and thus of the solid-fluid coupling can be evident in the responses. As can be seen in the figure, two peaks with significant amplitude can be observed, corresponding to vibration modes of the wall coupled to the fluid; after these peaks, the response exhibits a smoother form. Figures 16(b) and 16(c) illustrate the absolute difference, to the reference solution, calculated for the three different approaches, namely, the MFS-FEM, the MFS-CM, and the MFS-MLPG. For the fluid, 10 nodes are positioned along the interface, while for the solid, results are presented for 20 (Figure 16(b)) and 40 nodes (Figure 16(c)). When 20 nodes are used, the responses provided by the three approaches are very similar, with the two meshless methods exhibiting a lower error level at the lower frequencies, and with a worse behaviour of the CM being observable in the higher frequencies. Observing the figure, it is apparent that the MLPG is providing a more accurate response throughout the analysed frequency range, exhibiting a lower error than the FEM even at high frequencies. When more nodes are used (Figure 16(c)), the error levels provided by all methods improve, although the MLPG still exhibits a better overall behaviour than the remaining methods.(a) Reference pressure result; absolute error considering (b) 20 and (c) 40 boundary nodes in the solid, along the solid-fluid interface.
(a)
(b)
(c)It is important to notice that, although different error levels are observed for each method, the iterative coupling algorithm always quickly converges even for frequencies in the vicinity of the response peaks referred to before. Figure17 illustrates the number of iterations required for convergence, for the three approaches, considering 40 nodes along the interface to model the solid. Clearly, all three approaches exhibit very similar curves, requiring similar numbers of iterations for the iterative process to reach convergence at each frequency. It is also very clear that, for this case, the number of iterations is always small, slightly exceeding 20 iterations only at a few specific frequencies. For the remaining frequencies, only about 10 to 15 iterations are necessary to attain convergence. In the same figure, the number of iterations required, when a fixed relaxation parameter is used (i.e., λ
=
0.5), is also depicted. Comparison between the curves calculated with optimal and fixed parameters reveals a striking difference, with the optimal parameter always leading to significantly less iterations and ensuring convergence for all frequencies. When the relaxation parameter is fixed, convergence cannot be reached at two sets of frequencies, associated with specific dynamic behaviours of the system. These results once again illustrate the importance of using a well-chosen relaxation parameter to ensure that effective analyses are obtained.Number of iterations required for convergence: (a) MFS-FEM; (b) MFS-CM; (c) MFS-MLPG.
(a)
(b)
(c)The set of plots shown in Figure18 illustrates the deformation of the structure at frequency 125 Hz. In the plotted results, a grey patch is used to identify the reference response, while marks are used to depict the (amplified) deformed shape of the structure when analysed by the three iteratively coupled approaches. The left column reveals the response for 20 nodes positioned in the solid along the interface, whereas the right column shows the equivalent result computed for 40 nodes. It can be observed that, for all approaches, the response improves significantly when more nodes are used, indicating that the convergent behaviour of the methods can be observed. The computed responses reveal very similar shape and displacement amplitudes when compared to the reference solution, with the response provided by the MFS-CM approach being somewhat worse than the remaining two. In fact, the MFS-MLPG and the MFS-FEM exhibit very similar behaviours, with lower discrepancies being registered for the meshless method. Finally, Figure 19 illustrates the variation of the complex relaxation parameter throughout the iterative process, showing its real and imaginary components, together with its absolute value. Those plots reveal a very similar evolution of the parameter for all combinations of methods; again, this indicates that the discussed iterative procedure is quite independent of the discretization methods involved in the analyses.Real ((a) and (b)) and imaginary ((c) and (d)) parts of the deformation of the solid structure (amplified) when the excitation frequency is 125 Hz. Results are shown for 20 ((a) and (c)) and 40 ((b) and (d)) boundary nodes in the solid, along the fluid-solid interface, when using the FEM (×), CM (∘), and MLPG (•).
(a)
(b)
(c)
(d)Variation of the complex relaxation parameter throughout the iterative process: (a) MFS-FEM; (b) MFS-CM; (c) MFS-MLPG.
(a)
(b)
(c)


## 4.1. Electromagnetic Waves

In electromagnetic models, vectorial wave equations describe the electric and the magnetic field evolution [56, 57]. In this case, (1a) can be rewritten as (in this subsection, time-domain analyses are focused on):(12a)
∇
×
(
μ
(
x
)
-
1
∇
×
E
(
x
,
t
)
)
+
ε
(
x
)
E
¨
(
x
,
t
)
=
-
J
˙
(
x
,
t
)
,
(12b)
∇
×
(
ε
(
x
)
-
1
∇
×
H
(
x
,
t
)
)
+
μ
(
x
)
H
¨
(
x
,
t
)
=
∇
×
(
ε
(
x
)
-
1
J
(
x
,
t
)
)
,and (3) can be rewritten as(13a)
n
(
x
)
×
(
E
(
x
+
,
t
)
-
E
(
x
-
,
t
)
)
=
0
,
(13b)
(
D
(
x
+
,
t
)
-
D
(
x
-
,
t
)
)
·
n
(
x
)
=
ρ
(
x
,
t
)
,
(13c)
(
B
(
x
+
,
t
)
-
B
(
x
-
,
t
)
)
·
n
(
x
)
=
0
,
(13d)
n
(
x
)
×
(
H
(
x
+
,
t
)
-
H
(
x
-
,
t
)
)
=
J
(
x
,
t
)
,where E and H are the electric and magnetic field intensity vectors, respectively; D and B represent the electric and magnetic flux densities, respectively; and J and ρ stand for the electric current and electric charge density, respectively. The parameters ε and μ denote, respectively, the permittivity and permeability of the medium and its wave propagation velocity is specified as c
=
(
ε
μ
)
-
1
/
2. n is the normal vector, from domain 1 to domain 2. Equations (13a) and (13b) state that the tangential component of E is continuous across the interface and that the normal component of D has a step of surface charge on the interface surface, respectively. Equations (13c) and (13d) state that the normal component of B is continuous across the interface and that the tangential component of H is continuous across the interface if there is no surface current present, respectively.In the present application, the electromagnetic fields surrounding infinitely long wires are studied [41]. Two cases of analysis are focused here, namely, (a) case 1, where one wire is considered; (b) case 2, where two wires are employed. For both cases, the wires are carrying time-dependent currents (i.e., I
(
t
)
=
t or I
(
t
)
=
t
2) and they are located along the adopted z-axis. A sketch of the model is depicted in Figure 3.Sketch of the electromagnetic models and adopted FEM/BEM spatial discretizations: (a) case 1, one wire; (b) case 2, two wires.
(a)
(b)The spatial and temporal evolution of the electric field intensity vector is analyzed here taking into account a finite element method (FEM)—boundary element method (BEM) coupled formulation. In this context, the FEM is applied to model the region close to the wires, whereas the BEM simulates the remaining infinity domain. As it is well known, the BEM employs fundamental solutions which fulfill the radiation condition. Thus, this formulation is very suitable to perform infinite domain analysis, once reflected waves from infinity are avoided [58].The adopted spatial discretization is also described in Figure3. In this case, 2344 linear triangular finite elements and 80 linear boundary elements are employed in the analyses (see references [57, 58] for more details regarding the FEM and the BEM applied to electromagnetic analyses). The radius of the FEM-BEM interface is defined by R
=
1
m and matching nodes are considered at the interface. For temporal discretization, the selected time step is given by Δ
t
=
5
·
1
0
-
11
s for both domains. The physical properties of the medium (air) are μ
=
1.2566
·
1
0
-
6
H
/
m and ε
=
8.8544
·
1
0
-
12
F
/
m.Figure4 shows the modulus of the electric field intensity obtained at points A and B (see Figure 3) considering the iterative coupling methodology. Analytical time histories [58] are also depicted in Figure 4, highlighting the good accuracy of the numerical results. In Figure 5, charts are displayed, indicating the percentage of occurrence of different relaxation parameter values (evaluated according to expression (10)), in each analysis. As can be observed, for all considered cases, optimal relaxation parameters are mostly in the interval 0.7
≤
λ
≤
0.8. In fact, an optimal relaxation parameter selection is extremely case dependent. It is function of the physical properties of the model, geometric aspects, adopted spatial and temporal discretizations, and so forth. Equation (10) provides a simple expression to evaluate this complex parameter.Time history results for the electric field intensity at points A and B consideringI
(
t
)
=
t and (a) case 1 and (b) case 2; I
(
t
)
=
t
2 and (c) case 1 and (d) case 2.
(a)
(b)
(c)
(d)Percentage of occurrence of different relaxation parameter values during the analysis, consideringI
(
t
)
=
t and (a) case 1 and (b) case 2; I
(
t
)
=
t
2 and (c) case 1 and (d) case 2.
(a)
(b)
(c)
(d)In order to illustrate the effectiveness of the methodology when considering different time discretizations for different domains, Figure6 depicts results that are computed considering Δ
t
=
2.5
·
1
0
-
11
s for the FEM and Δ
t
=
2.0
·
1
0
-
10
s for the BEM (i.e., a difference of 8 times between the time steps). For simplicity, results are presented considering just the first case of analysis, that is, case 1 and I
(
t
)
=
t. As one can observe in Figure 6(a), good results are still obtained taking into account the iterative formulation, in spite of the existing time disconnections at the interface. In Figure 6(b), the evolution of the relaxation parameter is depicted, taking into account this last configuration. As one can observe, in this case, optimal relaxation parameter values are between 0.7 and 1.0 and mostly concentrate on the interval (
0.9,1.0
). In fact, it is expected that these values get closer to 1.0 when smaller time steps are considered. In the present analysis, an average number of 4.92 iterations per time step is obtained (taking into account 800 FEM time steps), which is a relatively low number, illustrating the good performance of the technique (it must be remarked that a tight tolerance criterion was adopted for the convergence of the iterative analysis).Results considering different time steps for each domain: (a) electric field intensity at points A and B; (b) optimal relaxation parameters for each iterative step.
(a)
(b)


## 4.2. Acoustic Waves

In acoustic models, a scalar wave equation describes the acoustic pressure field evolution [1]. In this case, (1b) can be rewritten as (in this subsection, frequency-domain analyses are focused)
(14)
∇
·
(
κ
(
x
)
∇
p
(
x
,
ω
)
)
+
ω
2
ρ
(
x
)
p
(
x
,
ω
)
=
γ
(
x
,
ω
)
and (3) can be rewritten as(15a)
(
p
(
x
+
,
ω
)
-
p
(
x
-
,
ω
)
)
=
0
,
(15b)
(
q
(
x
+
,
ω
)
-
q
(
x
-
,
ω
)
)
=
g
(
x
,
ω
)
,where p is the hydrodynamic pressure and γ and g stand for domain and surface sources, respectively. The parameters ρ and κ denote, respectively, the mass density and compressibility of the medium and its wave propagation velocity is specified as υ
=
(
κ
/
ρ
)
1
/
2. The hydrodynamic fluxes on the interfaces are represented by q, and they are defined by q
=
κ
∇
p
·
n, where n is the normal vector, from domain 1 to domain 2. Equation (15a) states that the pressure is continuous across the interface, whereas (15b) states that the flux is continuous across the interface if there is no surface source.The advantages of using iterative coupling procedures are revealed when more complex configurations are analyzed. In this subsection, the case of a heterogeneous domain, composed of a homogeneous fluid incorporating multiple circular inclusions with different properties, is analyzed. For this purpose, consider the host medium to allow the propagation of sound with a velocity of 1500 m/s, and this medium is excited by a line source located atx
s
=
-
5
.
0 m and y
s
=
0.0 m. Within this fluid, consider the presence of 8 circular inclusions; all of them are with unit radius and filled with a different fluid, allowing sound waves to travel at 3000 m/s, as depicted in Figure 7.(a) Sketch for the heterogeneous medium with multiple subregions; (b) boundary and domain point distribution, considering the spatial discretization of an inclusion and adjacent fluid.
(a)
(b)The above-described system has been analyzed taking into account the proposed iterative coupling procedure making use of the Kansa’s method (KM) to model all the inclusions and of the method of fundamental solutions (MFS) to model the host fluid (see references [12, 59–61] for more details regarding the KM and the MFS applied to acoustic analyses). One should note that, since the real source is positioned at the outer region, the iterative process is initialized with the analysis of the MFS model, considering prescribed Neumann boundary conditions at the common interfaces. Once the boundary pressures for the outer region are computed, these values are transferred to the closed regions by imposing Dirichlet boundary conditions, incorporating information about the influence of each inclusion on the remaining heterogeneities. Then, each KM subregion is analysed independently and the internal boundary values (normal fluxes) are evaluated autonomously for each inclusion. The iterative procedure then goes further, including the calculation of the relaxation parameter at each iterative step, as well as the correction of boundary variables, until convergence is achieved.To model the system, each MFS boundary is discretized by 55 points. 331 KM domain points are equally distributed within each inclusion, and 66 KM boundary points (around 31 points per wavelength) are used (see Figure7(b)). The complexity of the model hinders the definition of a closed form solution; thus, the results are checked against a numerical model which performs the direct (i.e., noniterative) coupling between both methods. In that model, 66 boundary points are used in the MFS to define the boundary of each inclusion, and 66 and 331 KM boundary and domain points, respectively, are adopted for the discretization of each inclusion (analogously to the iterative coupling procedure). Figure 8 compares the responses computed by the iterative and the direct coupling methodologies. Results are depicted along the boundary of the 8th inclusion, for excitation frequencies of 400 Hz (Figure 8(a)) and 1000 Hz (Figure 8(b)). As can be observed in the figure, there is a perfect match between both approaches, with the iterative procedure clearly converging to the correct solution.Hydrodynamic pressures along the common interface of the 8th inclusion for (a)ω
=
400 Hz and (b) ω
=
1000 Hz (—, real-iterative; - - - -, Imag-iterative; ○, Real-direct; □, Imag-direct).
(a)
(b)It is important, at this point, to highlight the differences in the computational times of the direct and of the iterative coupling approaches. For the present model configuration, the direct coupling approach had to deal simultaneously with 528 boundary points and a total of 2648 internal points (i.e., considering a coupled matrix of dimension 3704), which is implied in 373.89 s of CPU time in a Matlab implementation (being this CPU time independent of the frequency in focus). For the iterative coupling approach, using 55 boundary points for the MFS and 66 boundary points for the KM, it was possible to obtain analogous results considering 12.06 s of CPU time for the frequency of 200 Hz and 32.32 s for the frequency of 1000 Hz (i.e., 3.23% and 8.64% of the computational cost of the direct coupling methodology, resp.). Even if the same number of boundary points is used in the iterative coupling approach for the MFS (i.e., 66 points), the final CPU time would just increase up to 35.71 s (9.55% of the computational cost of the direct coupling methodology). These results are summarized in Table1, where the number of iterations and the CPU time are presented for the first scenarios (i.e., 55 boundary points for the MFS) and for frequencies between 50 Hz and 1000 Hz. The values described in Table 1 further confirm that the difference in calculation times between the iterative and the direct coupling approach is striking and reveal an excellent gain in performance favouring the iterative coupling technique. It is important to understand that this gain is strongly related to the possibility of dealing with smaller-sized matrices when using the iterative coupling procedure. Moreover, it is possible to invert (or triangularize, etc.) the relevant matrices only at the first iterative step and then proceed with the calculations using the inverted matrices (or forward/back substituting, etc.). As a consequence, after the first iteration, only matrix vector multiplication operations are required, and very high savings in what concerns computational time are achieved. In fact, considering that the number of operations required for matrix inversion can be assumed to be of the order of N
3 (N being the matrix size), a simple calculation allows concluding that, for the current model, the relative cost of inverting the eight KM matrices (each one being a square matrix with 397 × 397 entries) and the MFS matrix (with 440 × 440 entries) is less than 2% of the cost of inverting a larger 3704 × 3704 matrix, as required for the direct coupling strategy. Similar conclusions can be obtained considering other solver procedures, such as matrix triangularizations, demonstrating that a considerably less expensive methodology is obtained if the different subdomains are analysed separately (even considering an eventual high number of iterative steps in the iterative analysis).Table 1
Total number of iterations and relative CPU time (iterative/direct coupling) for the acoustic model.
Frequency (Hz)
Iterations
Relative CPU time (%)
50
14
2.09
100
18
2.43
150
29
3.19
200
31
3.23
250
26
2.78
300
40
3.86
350
31
3.13
400
32
3.20
450
30
3.02
500
49
4.49
550
48
4.40
600
45
4.18
650
68
5.94
700
34
3.28
750
39
3.66
800
110
9.20
850
98
8.28
900
78
6.75
950
104
8.83
1000
100
8.64Analyzing the difference in computational times between the two analyzed frequencies (i.e., 200 Hz and 1000 Hz) also reveals a significant difference between them. This difference is related to the number of iterations required for convergence, which was higher when the excitation frequency of 1000 Hz was considered. The plot in Figure9 indicates the number of iterations required for convergence along a range of frequencies between 10 Hz and 1000 Hz, using a constant number of boundary (55 for the MFS and 66 for the KM) and internal points (331 for the KM). As expected, the number of iterations increases with the frequency. It is interesting to note that the maximum necessary number of iterations occurred for a frequency of 990 Hz, requiring 170 iterations and a CPU time of 54.80 s to converge, which is less than 15% of the CPU time required by the direct coupling for the same frequency.Figure 9
Total number of iterations considering different frequencies and 55 collocation points for the MFS and 66 boundary points for the KM, at each circular inclusion.In Figure10, the wavefield produced within and around the inclusions is illustrated for excitation frequencies of 600 Hz and 1000 Hz. As expected, as the frequency increases, the multiple inclusions generate progressively more complex wave fields, with the interaction between them becoming very significant for the higher frequency. Observation of these results also reveals a strong shadow effect produced by the inclusions, with much lower amplitudes being registered in the region behind the inclusions placed further away from the source. This effect is even more pronounced for the higher frequency. Interestingly, for both frequencies, the space between the two lines of inclusions works as a guiding path, along which the sound energy travels with less attenuation.3D plots of the sound field for frequencies of (a) 600 Hz and (b) 1000 Hz.
(a)
(b)


## 4.3. Mechanical Waves

In dynamic models, a vectorial wave equation describes the displacement field evolution [1]. In this case, considering linear behaviour, (1a) can be rewritten as (in this subsection, time-domain analyses are focused)
(16)
∇
×
(
μ
(
x
)
∇
×
u
(
x
,
t
)
)
-
∇
(
(
η
(
x
)
+
2
μ
(
x
)
)
∇
·
u
(
x
,
t
)
)
+
ρ
(
x
)
u
¨
(
x
,
t
)
=
f
(
x
,
t
)
,
and (3) can be rewritten as:(17a)
(
u
(
x
+
,
t
)
-
u
(
x
-
,
t
)
)
=
0
,
(17b)
(
σ
(
x
+
,
t
)
-
σ
(
x
-
,
t
)
)
n
(
x
)
=
τ
(
x
,
t
)
,where u is the displacement vector and f and τ stand for domain and surface forces, respectively. The terms η and μ denote the so-called Lamé parameters, and ρ is the mass density of the medium. In this case, the wave propagation velocities are specified as c
s
=
(
μ
/
ρ
)
1
/
2 (shear wave) and c
d
=
(
(
η
+
2
μ
)
/
ρ
)
1
/
2 (dilatational wave). The stress tensor is denoted by σ and n is the normal vector, from domain 1 to domain 2. Equation (17a) states that the displacements are continuous across the interface, whereas (17b) states that the tractions are continuous across the interface if there are no surface forces on it.One main advantage of the discussed coupling algorithm is that other iterative processes can be carried out in the same iterative loop needed for the coupling. Thus, consideration of coupled nonlinear models, as for example, may not demand a superior amount of computational effort, which is very beneficial.In the present application, a nonlinear model is considered and elastoplastic analyses are carried out (for details about elastoplastic analyses, one is referred to [33–35, 62–64]). Moreover, two discretization approaches are employed here, one taking into account FEM-BEM coupling procedures, and another considering BEM-BEM coupled techniques (for more details about these coupled models, one is referred to [37, 43]). In this context, a nonlinear infinity domain is analyzed here, in which a circular cavity is loaded. The region expected to develop plastic strains is discretized by the finite element method, in the case of the FEM-BEM coupled analysis, or by the domain boundary element method (D-BEM), in the case of the BEM-BEM coupled analysis. The remainder of the infinity domain is discretized by the time-domain boundary element method (TD-BEM). A sketch of the model is depicted in Figure 11, as well as the adopted discretizations.(a) Sketch of the circular cavity; (b) FEM-BEM discretization; (c) BEM-BEM discretization.
(a)
(b)
(c)The FEM-BEM discretization is depicted in Figure11(b). In this case, 1944 linear triangular finite elements and 80 linear boundary elements are employed in the coupled analysis. The BEM-BEM discretization is depicted in Figure 11(c). In this case, 46 linear boundary elements are employed in the BEM-BEM coupled analysis (20 linear boundary elements for the TD-BEM and 26 linear boundary elements for the D-BEM), as well as 270 linear triangular cells (D-BEM formulation). In the BEM-BEM coupled analysis, the double symmetry of the problem is taken into account. An interesting feature of the boundary element formulation is that symmetric bodies under symmetric loads can be analysed without discretization of the symmetry axes. This can be accomplished by an automatic condensation process, which integrates over reflected elements and performs the assemblage of the final matrices in reduced size [64]. The time discretization adopted is given by Δ
t
=
0.04
s for the FEM and Δ
t
=
0.20
s for the D-BEM and the TD-BEM.The physical properties of the model areμ
=
2.652
·
1
0
8
N
/
m
2, η
=
2.274
·
1
0
8
N
/
m
2, and ρ
=
1.804
·
1
0
3
kg
/
m
3. A perfectly plastic material obeying the Mohr-Coulomb yield criterion is assumed, where c
o
=
4.8263
·
1
0
6
N
/
m
2 (cohesion) and ϕ
=
3
0
° (internal friction angle). The geometry of the problem is defined by R
=
3.048
m (the radius of the TD-BEM circular mesh is given by 5R).In Figure12, the displacement time history at point A is depicted, considering linear and nonlinear analyses. As one can notice, good agreement is observed between the FEM-BEM and BEM-BEM results. It is important to highlight that, for the FEM-BEM analyses, a difference of 5 times between the FEM and BEM time steps is considered, illustrating the effectiveness of the time interpolation/extrapolation procedures adopted in the analyses.Figure 12
Radial displacement considering coupled TD-BEM/D-BEM and coupled TD-BEM/FEM analyses: linear and nonlinear results at point A.The number of iterations per time step and the optimal relaxation parameters, evaluated at each iterative step, are depicted in Figure13, taking into account the FEM-BEM coupled analyses. As one may observe, basically the same computational effort (i.e., number of iterative steps) is necessary for both linear and nonlinear analyses, highlighting the efficiency of the proposed methodology for complex phenomena modeling. It is also important to remark the low number of iterative steps necessary for convergence, with a maximum of 7 iterations being necessary, within a time step, taking into account the entire linear and nonlinear analyses. For the focused configurations, the optimal relaxation parameters are intricately distributed within the interval (
0.75
;
1.00
), as depicted in Figure 13(b).TD-BEM/FEM analyses: (a) number of iterations per time step considering optimal relaxation parameters; (b) optimal relaxation parameters for each iterative step.
(a)
(b)In Table2, the total number of iterations is presented, considering analyses with optimal relaxation parameters and with some constant preselected λ values. As one may observe, an inappropriate selection for the relaxation parameter can considerably increase the associated computational effort. Thus, the optimization technique is extremely important in order to provide a robust and efficient iterative coupling formulation. In Figure 14, the computed σ
x
y stresses are depicted, considering the BEM-BEM elastoplastic analysis. An advantage of the D-BEM is that it employs nodal stress equations [37], allowing computing continuous stress fields, in counterpart to the FEM, which computes stresses based on displacement derivatives, obtaining discontinuous stress fields at element interfaces.Table 2
Total number of iterations (considering all time steps) for the dynamic model.
Relaxation parameter
Elastic analysis
Elastoplastic analysis
1.00
3730
3740
0.90
3392
3443
0.80
3973
3993
0.70
4772
4777
Optimal
3287
3346Figure 14
Spatial and temporal evolution ofσ
x
y considering elastoplastic analysis (scale factor: 6.8947
·
10
6
N
/
m
2); (a) t
=
4 s; (b) t
=
8 s; (c) t
=
12 s; (d) t
=
16 s.


## 4.4. Coupled Acoustic-Mechanical Waves

In this case, different wave equations, as indicated in Sections4.2 and 4.3 (see (14) and (16)), describe different domains of the global model. The interface conditions for the acoustic-dynamic coupling (3) can then be written as (in this subsection, frequency-domain analyses are focused)(18a)
(
n
(
x
)
·
σ
(
x
+
,
ω
)
n
(
x
)
-
p
(
x
-
,
ω
)
)
=
0
,
(18b)
(
-
ω
2
n
(
x
)
·
u
(
x
+
,
ω
)
-
υ
(
x
-
)
2
q
(
x
-
,
ω
)
)
=
0
,where u is the displacement vector and σ is the stress tensor of the dynamic model (domain 1). p is the hydrodynamic pressure and q is the hydrodynamic flux of the acoustic model (domain 2). n is the normal vector, from domain 1 to domain 2. The acoustic wave propagation velocity is denoted by υ. Equation (18a) states that the normal components of the dynamic tractions are equal to the acoustic pressures and (18b) relates the normal components of the dynamic accelerations to the acoustic fluxes.In the present application, a model in which a concrete wall of 10.0 m high is coupled to a fluid waveguide, filled with water, is analyzed. A sketch of the model is depicted in Figure15(a). For this case, a pressure source is positioned in the waveguide, at (−10.0; 0.5), illuminating the system. The concrete structure corresponds to a wall with variable cross-section, exhibiting thicknesses of 4.0 m at its basis and of 2.0 m at its top.(a) Sketch of the coupled acoustic-dynamic model; (b) MFS collocation points (o) and virtual sources (x); (c) FEM mesh when 20 nodes are used along the solid-fluid interface; (d) node distribution for the meshless methods when 20 nodes are used along the solid-fluid interface.
(a)
(b)
(c)
(d)To simulate this coupled system, several approaches are employed. For the fluid medium, the MFS is used in all cases, allowing the use of the Green’s function for a waveguide (see [61] for details concerning this function). This Green’s function is written as a summation of modes, and its convergence is very difficult when the source and the receiver are positioned along the same vertical line, thus posing severe difficulties for its use together with a BEM formulation. The structure is modelled using three different methods, namely, the FEM, a local collocation method (CM) (see [65] for details about this procedure), and a meshless local Petrov-Galerking technique (MLPG) (see [65, 66] for details about this procedure). Representations of the node distribution for each one of the methods can be found in Figures 15(b)–15(d).Since no analytical solution can be found in the literature for the present case, a numerical solution making use of a full BEM model ensuring the correct coupling between the solid and the fluid is used. For this case, the rigid bottom of the waveguide is accounted for using an image-source Green’s function, while the free surface is fully discretized up to a distance of 60.0 m from the concrete wall; after this, an anechoic termination is considered, imposing adequate Robin boundary conditions. For the solid, full-space Green’s functions are adopted, and 60 nodal points are used along the solid-fluid interface, ensuring that accurate results can be obtained. A total of 796 boundary elements are used to build this model. Details on the mathematical formulation of this technique can be found in the works of Tadeu and Godinho  [7].Figure16(a) illustrates the reference response in terms of real and imaginary components of the acoustic pressure at a receiver located at (−1.0; 3.0), for frequencies between 1 Hz and 150 Hz. This position is chosen so that the effect of the vibration of the concrete wall and thus of the solid-fluid coupling can be evident in the responses. As can be seen in the figure, two peaks with significant amplitude can be observed, corresponding to vibration modes of the wall coupled to the fluid; after these peaks, the response exhibits a smoother form. Figures 16(b) and 16(c) illustrate the absolute difference, to the reference solution, calculated for the three different approaches, namely, the MFS-FEM, the MFS-CM, and the MFS-MLPG. For the fluid, 10 nodes are positioned along the interface, while for the solid, results are presented for 20 (Figure 16(b)) and 40 nodes (Figure 16(c)). When 20 nodes are used, the responses provided by the three approaches are very similar, with the two meshless methods exhibiting a lower error level at the lower frequencies, and with a worse behaviour of the CM being observable in the higher frequencies. Observing the figure, it is apparent that the MLPG is providing a more accurate response throughout the analysed frequency range, exhibiting a lower error than the FEM even at high frequencies. When more nodes are used (Figure 16(c)), the error levels provided by all methods improve, although the MLPG still exhibits a better overall behaviour than the remaining methods.(a) Reference pressure result; absolute error considering (b) 20 and (c) 40 boundary nodes in the solid, along the solid-fluid interface.
(a)
(b)
(c)It is important to notice that, although different error levels are observed for each method, the iterative coupling algorithm always quickly converges even for frequencies in the vicinity of the response peaks referred to before. Figure17 illustrates the number of iterations required for convergence, for the three approaches, considering 40 nodes along the interface to model the solid. Clearly, all three approaches exhibit very similar curves, requiring similar numbers of iterations for the iterative process to reach convergence at each frequency. It is also very clear that, for this case, the number of iterations is always small, slightly exceeding 20 iterations only at a few specific frequencies. For the remaining frequencies, only about 10 to 15 iterations are necessary to attain convergence. In the same figure, the number of iterations required, when a fixed relaxation parameter is used (i.e., λ
=
0.5), is also depicted. Comparison between the curves calculated with optimal and fixed parameters reveals a striking difference, with the optimal parameter always leading to significantly less iterations and ensuring convergence for all frequencies. When the relaxation parameter is fixed, convergence cannot be reached at two sets of frequencies, associated with specific dynamic behaviours of the system. These results once again illustrate the importance of using a well-chosen relaxation parameter to ensure that effective analyses are obtained.Number of iterations required for convergence: (a) MFS-FEM; (b) MFS-CM; (c) MFS-MLPG.
(a)
(b)
(c)The set of plots shown in Figure18 illustrates the deformation of the structure at frequency 125 Hz. In the plotted results, a grey patch is used to identify the reference response, while marks are used to depict the (amplified) deformed shape of the structure when analysed by the three iteratively coupled approaches. The left column reveals the response for 20 nodes positioned in the solid along the interface, whereas the right column shows the equivalent result computed for 40 nodes. It can be observed that, for all approaches, the response improves significantly when more nodes are used, indicating that the convergent behaviour of the methods can be observed. The computed responses reveal very similar shape and displacement amplitudes when compared to the reference solution, with the response provided by the MFS-CM approach being somewhat worse than the remaining two. In fact, the MFS-MLPG and the MFS-FEM exhibit very similar behaviours, with lower discrepancies being registered for the meshless method. Finally, Figure 19 illustrates the variation of the complex relaxation parameter throughout the iterative process, showing its real and imaginary components, together with its absolute value. Those plots reveal a very similar evolution of the parameter for all combinations of methods; again, this indicates that the discussed iterative procedure is quite independent of the discretization methods involved in the analyses.Real ((a) and (b)) and imaginary ((c) and (d)) parts of the deformation of the solid structure (amplified) when the excitation frequency is 125 Hz. Results are shown for 20 ((a) and (c)) and 40 ((b) and (d)) boundary nodes in the solid, along the fluid-solid interface, when using the FEM (×), CM (∘), and MLPG (•).
(a)
(b)
(c)
(d)Variation of the complex relaxation parameter throughout the iterative process: (a) MFS-FEM; (b) MFS-CM; (c) MFS-MLPG.
(a)
(b)
(c)


## 5. Conclusions

This paper presents an overview of the application of iterative coupling strategies to the analysis of wave propagation problems. Different methods were considered, including mesh-based and meshless methods, ranging from the more classic BEM and FEM to the less usual MLPG, collocation methods, or MFS. Several examples of the iterative coupling technique were presented in Section4, including the application of the scheme to electromagnetic, acoustic, elastic/elastoplastic, and acoustic-elastic interaction problems. The generality and flexibility of the iterative scheme allowed an efficient analysis of these problems, either using time or frequency domain models. The use of an optimized relaxation parameter (which is the basis of this scheme) proved to be quite important, clearly accelerating (or in some cases ensuring) convergence; for all tested cases, this parameter was shown to unpredictably vary throughout the iterative process, and thus its appropriate recalculation at each iterative step becomes important. The illustrated analyses clearly indicated that the strategy can be effectively used for different methods and that the performance of the iterative technique is quite insensitive to the discretization methods employed in the analyses.It should be highlighted that the coupling technique presented here is based on previous experience and works of the authors and, although the paper is focused in wave propagation problems, the iterative strategy presently discussed can be regarded as a quite generic framework to perform the coupling between different methods in many types of applications.


## Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.



---

*Source: 426283-2014-01-14.xml*