# Adaptive Analysis of Acoustic-Elastodynamic Interacting Models Considering Frequency Domain MFS-FEM Coupled Formulations

Authors: D. Soares; L. Godinho
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This work discusses adaptive iterative coupling strategies for the frequency domain analysis of interacting acoustic-elastodynamic models. The method of fundamental solutions (MFS) is used to analyze acoustic fluids, whereas the finite element method (FEM) is employed to discretize elastodynamic solids. Flexible and optimized iterative MFS-FEM coupling procedures are considered, allowing independent discretizations to be adopted for each subdomain. In this context, it is easy to implement adaptive refinements and enable enhanced analyses. Two adaptive coupling approaches are discussed, based on multiple and single iterative algorithms. Numerical results are presented to illustrate the performance of the proposed techniques.


---

## Body

## 1. Introduction

The analysis of coupled acoustic-elastodynamic systems is a complex and demanding task, which involves different physical phenomena and governing equations. It can be more efficiently handled by making use of different numerical methods for each part of the problem. A number of the research works that have been published on this topic have looked at BEM-FEM coupling algorithms [1–8]. Indeed, the BEM (boundary element method) is quite suitable to handle infinite or semi-infinite media, while the FEM (finite element method) has been widely studied and developed to handle structural systems; a combination of both methods thus seems like a natural choice.Many of the above references account for the full-coupling between the BEM and the FEM directly by establishing a fully coupled system matrix [1–5]; however, this option can be quite inefficient. This is because the symmetric and sparse/banded character of the FEM matrix is significantly affected by the presence of the BEM components, and optimized solvers, usually used by the FEM, cannot be employed anymore. Additionally, the different properties of fluid and solid media may lead to badly conditioned matrices, which affect the accuracy of the methodology. Finally, standard direct coupling methodologies require compatible discretizations (with matching nodes along the interfaces) and this greatly affects the flexibility and versatility of the technique. Although alternatives can be used to allow incompatible discretizations (such as interpolation along the interface), this usually increases the complexities when setting up the problem.In order to overcome these difficulties, iterative coupling procedures have been presented, considering both time and frequency domain acoustic-elastodynamic interaction analyses, and taking into account not only BEM-FEM coupling, but also considering other possibilities [6–10]. Studies have shown that iterative coupling approaches allow BEM and FEM subdomains to be analyzed separately, leading to smaller and better-conditioned systems of equations (different solvers, suitable for each subdomain, may be employed). Furthermore, a small number of iterations are usually required for the algorithm to converge and the matrices related to the smaller governing systems of equations do not need to be treated (inverted, triangularized, etc.) at each iterative step, thus providing an efficient methodology. Different methods have been tested for frequency domain analysis. For the elastodynamic subdomain, the FEM [8, 10] and the newer meshless methods, including the Kansa’s method [9, 10] and the meshless local Petrov-Galerkin method [10], have been employed with success, while for the acoustic fluid, both the BEM [8, 10] and the method of fundamental solutions (MFS) [9, 10] have been used in iterative coupling algorithms. Previous works [9–12] have also shown that the MFS can be particularly well suited for the analysis of the acoustic subdomain, even surpassing the BEM in terms of efficiency and accuracy. Thus, iteratively coupling the MFS with the FEM can be a robust and efficient alternative strategy to handle acoustic-elastodynamic problems in the frequency domain. It is worth noting that although the MFS is quite well suited for this type of analysis, it requires previously defining the position of a set of virtual sources, located outside the domain of analysis, which are used to construct the approximation of the acoustic pressure field. The definition of such positions is still an open problem, which has been widely discussed in the literature, and some methods have been proposed for their optimal determination, such as in [13]; a number of alternative methods have also been proposed to overcome this difficulty, such as those described in [14–16].Differently from earlier studies, the present paper reports on the analysis of acoustic-elastodynamic systems using an efficient adaptive iterative coupling between the MFS and the FEM. In fact, adaptive remeshing can be employed in association with the analyses of different physical problems. For the case of elastoplasticity, remeshing has been applied in the context of the BEM-FEM coupled analysis by Elleithy et al. [17, 18], who updated the subdomains modelled by the FEM and the BEM according to the evolution of the plastic zones of the model. Regarding dynamic problems, although contributions can be found in the literature on elastodynamics (such as [19–21]) and acoustics (such as [22, 23]), not many studies report on adaptive methods for coupled acoustic-elastodynamic interaction problems. Of the few there are, we would mention García et al. [24], who applied an adaptive FEM algorithm to solve solid-fluid interaction problems governed by the Navier-Stokes equation and Demkowicz and Oden [25, 26], who applied an hp-adaptive coupled BEM-FEM model to study problems of elastic scattering, assuming the fluid to be governed by the Helmholtz equation.In these works, the BEM-FEM coupling is performed in a direct form, requiring the construction of a fully coupled matrix. Here, a very different approach is developed: the regions modelled by the MFS do not change during the analysis, and thus the fully populated matrices of the MFS are computed only once, which increases efficiency; as for the FEM, a coarse discretization is used to start the algorithm, and this is adaptively enriched as the solution evolves. In this context, and since refining only takes place within the FEM meshes, nonmatching nodes must be allowed at the MFS-FEM interfaces; otherwise changing the MFS point distributions would require the MFS matrices to be recomputed, with a consequent loss of efficiency. The use of an iterative coupling algorithm is thus highly desirable because of its superior performance and because it is very flexible.In this work, two different iterative coupling strategies are proposed and discussed, both offering mesh adaptivity. The first is a more standard approach in which the coupling between the two domains is computed iteratively, with the FEM mesh then being analyzed and refined according to the final coupled solution. The problem is then reanalyzed using the new mesh. The second approach, on the other hand, incorporates the mesh adaptivity within the coupling iterative process, thus producing a much more efficient procedure with significant computational savings while still providing the same accuracy level. For both cases, the use of different methods to address the acoustic and elastodynamic problems (MFS and FEM, respectively), which are handled separately, is a significant advantage of the iteratively coupled approaches. Indeed, the MFS allows directly accounting for the infinite character of the problem, while the FEM is quite adequate to model a finite solid subdomain, which may contain complex details and spatially variable properties. Thus, their joint use in iterative coupling strategies, without the requirement of matching discretizations in the contact interfaces, leads to a powerful numerical tool, which incorporates what can be called the “best of two worlds”.The paper is organized as follows: first, we present the governing equations of the acoustic and elastodynamic models and briefly discuss the basic aspects of the MFS and FEM. We then describe the proposed adaptive iterative MFS-FEM coupling algorithms and look at some numerical applications, illustrating the performance and potentialities of the proposed techniques.


## 2. Numerical Modelling

This section sets out the basic governing equations of the model and the main aspects of the numerical techniques focused upon. Here, acoustic fluids are handled by the MFS and elastodynamic solids by the FEM.


### 2.1. Acoustic Medium Discretization

The Helmholtz wave equation can be written as(1)px,ω,ii+κ2px,ω+sx,ω=0where p stands for the acoustic pressure distribution, s represents body source terms, κ=ω/c stands for the wavenumber, c represents the wave propagation velocity, and ω and x represent the frequency and the spatial domain of analysis, respectively. Subscript commas indicate partial space derivatives (index notation is adopted). The boundary conditions of the model are given by(2a)px,ω=p-x,ωforx∈Γp(2b)vx,ω=v-x,ωforx∈Γvwhere the prescribed values are indicated by overbars and v=-1/(iωρ)(∂p/∂n) represents the normal velocity of the acoustic fluid along the boundary, whose unit outward normal vector is represented by n (ρ stands for the mass density). The boundary of the model is denoted by Γ (Γ=Γp∪Γv) and the domain by Ω. Equation (2a) stands for essential (or Dirichlet) boundary conditions and (2b) stands for natural (or Neumann) boundary conditions.In the MFS, the solutionp is approximated by a linear combination of fundamental solutions centred on ns virtual sources (xs), placed outside the domain of interest to avoid singularities in the response:(3)px,ω=∑n=1nsGx,xns,ωαnω+ςx,ωIn (3), G stands for the Green’s function of the model, ς is related to domain source terms, and α represents coefficients to be determined. By applying approximation (3) to the model, a system of algebraic equations can be obtained, as indicated below:(4)Hωαω=βωwhere the entries of matrix H and vector β are given by(5a)Hmnω=Gmnωforxm∈Γp(5b)Hmnω=-1iωρ∂Gmnω/∂nmforxm∈Γv(6a)βmω=-ςmω+p-mωforxm∈Γp(6b)βmω=1iωρ∂ςmω/∂nm+v-mωforxm∈Γvand for 2D analysis the Green’s function expression can be given for an infinite medium by Gmn(ω) = G(xm,xns,ω) = -(i/4)H0(2)(κrmn), where H0(2) stands for the second type Hankel’s function of order 0, and the term rmn represents the distance between the collocation (xm) and the virtual source (xns) points.Once the system of (4) is solved (i.e., vector α is computed), the approximate solution at any point of interest can be obtained using definition (3). For more details considering the numerical modelling by the method of fundamental solutions, see [11, 27].


### 2.2. Solid Medium Discretization

The elastic wave equation is given by(7)cd2-cs2ujx,ω,ji+cs2uix,ω,jj+ω2uix,ω+bix,ω=0where ui and bi stand for displacement and body force distribution components, respectively. The notation for space derivatives employed in (1) is once again adopted. In (7), cd is the dilatational wave velocity and cs is the shear wave velocity; they are given by cd2=(γ+2μ)/ρ and cs2=μ/ρ, where ρ is the mass density and γ and μ are the Lamé constants. Equation (7) can be obtained from the combination of the following basic mechanical equations:(8a)σijx,ω,j+ρω2uix,ω+ρbix,ω=0(8b)σijx,ω=γδijεkkx,ω+2μεijx,ω(8c)εijx,ω=12uix,ω,j+ujx,ω,iwhere σij and εij are stress and strain tensor components, respectively, and δij is the Kronecker delta (δij = 1, for i=j and δij = 0, for i≠j). Equation (8a) is the momentum equilibrium equation; (8b) represents the constitutive law of the elastic model; and (8c) stands for kinematic relations. The boundary conditions of the model are given by(9a)uix,ω=u-ix,ωforx∈Γu(9b)τix,ω=τ-ix,ωforx∈Γτwhere, once again, the prescribed values are indicated by overbars and τi=σijnj denotes the traction vector along the boundary.In the FEM, the solutionui is approximated by a local (i.e., at element level) interpolation: (10)uix,ω=∑n=1neNnxuinωIn (10), N stands for local interpolating functions, ne describes the number of nodes in the element, and ui represents nodal values to be determined. Applying approximations (10), a system of algebraic equations is obtained once the FEM is implemented, as indicated below:(11)-ω2M+Kuω=fωwhere matrices M and K and vector f are given by(12a)M=⋃e∫ΩeNTρNdΩ(12b)K=⋃e∫ΩeBTDBdΩ(13)fω=⋃e∫ΩeNTbx,ωdΩ+∫ΓτeNTτ-x,ωdΓandN,B, andD stand for the interpolation, the strain, and the constitutive matrices of element e, respectively.Once the system of (11) is solved (i.e., vector u is computed), the approximate solution at any point of interest can be obtained using definition (10). For more details considering the numerical modelling by the finite element method, see [28–30].


## 2.1. Acoustic Medium Discretization

The Helmholtz wave equation can be written as(1)px,ω,ii+κ2px,ω+sx,ω=0where p stands for the acoustic pressure distribution, s represents body source terms, κ=ω/c stands for the wavenumber, c represents the wave propagation velocity, and ω and x represent the frequency and the spatial domain of analysis, respectively. Subscript commas indicate partial space derivatives (index notation is adopted). The boundary conditions of the model are given by(2a)px,ω=p-x,ωforx∈Γp(2b)vx,ω=v-x,ωforx∈Γvwhere the prescribed values are indicated by overbars and v=-1/(iωρ)(∂p/∂n) represents the normal velocity of the acoustic fluid along the boundary, whose unit outward normal vector is represented by n (ρ stands for the mass density). The boundary of the model is denoted by Γ (Γ=Γp∪Γv) and the domain by Ω. Equation (2a) stands for essential (or Dirichlet) boundary conditions and (2b) stands for natural (or Neumann) boundary conditions.In the MFS, the solutionp is approximated by a linear combination of fundamental solutions centred on ns virtual sources (xs), placed outside the domain of interest to avoid singularities in the response:(3)px,ω=∑n=1nsGx,xns,ωαnω+ςx,ωIn (3), G stands for the Green’s function of the model, ς is related to domain source terms, and α represents coefficients to be determined. By applying approximation (3) to the model, a system of algebraic equations can be obtained, as indicated below:(4)Hωαω=βωwhere the entries of matrix H and vector β are given by(5a)Hmnω=Gmnωforxm∈Γp(5b)Hmnω=-1iωρ∂Gmnω/∂nmforxm∈Γv(6a)βmω=-ςmω+p-mωforxm∈Γp(6b)βmω=1iωρ∂ςmω/∂nm+v-mωforxm∈Γvand for 2D analysis the Green’s function expression can be given for an infinite medium by Gmn(ω) = G(xm,xns,ω) = -(i/4)H0(2)(κrmn), where H0(2) stands for the second type Hankel’s function of order 0, and the term rmn represents the distance between the collocation (xm) and the virtual source (xns) points.Once the system of (4) is solved (i.e., vector α is computed), the approximate solution at any point of interest can be obtained using definition (3). For more details considering the numerical modelling by the method of fundamental solutions, see [11, 27].


## 2.2. Solid Medium Discretization

The elastic wave equation is given by(7)cd2-cs2ujx,ω,ji+cs2uix,ω,jj+ω2uix,ω+bix,ω=0where ui and bi stand for displacement and body force distribution components, respectively. The notation for space derivatives employed in (1) is once again adopted. In (7), cd is the dilatational wave velocity and cs is the shear wave velocity; they are given by cd2=(γ+2μ)/ρ and cs2=μ/ρ, where ρ is the mass density and γ and μ are the Lamé constants. Equation (7) can be obtained from the combination of the following basic mechanical equations:(8a)σijx,ω,j+ρω2uix,ω+ρbix,ω=0(8b)σijx,ω=γδijεkkx,ω+2μεijx,ω(8c)εijx,ω=12uix,ω,j+ujx,ω,iwhere σij and εij are stress and strain tensor components, respectively, and δij is the Kronecker delta (δij = 1, for i=j and δij = 0, for i≠j). Equation (8a) is the momentum equilibrium equation; (8b) represents the constitutive law of the elastic model; and (8c) stands for kinematic relations. The boundary conditions of the model are given by(9a)uix,ω=u-ix,ωforx∈Γu(9b)τix,ω=τ-ix,ωforx∈Γτwhere, once again, the prescribed values are indicated by overbars and τi=σijnj denotes the traction vector along the boundary.In the FEM, the solutionui is approximated by a local (i.e., at element level) interpolation: (10)uix,ω=∑n=1neNnxuinωIn (10), N stands for local interpolating functions, ne describes the number of nodes in the element, and ui represents nodal values to be determined. Applying approximations (10), a system of algebraic equations is obtained once the FEM is implemented, as indicated below:(11)-ω2M+Kuω=fωwhere matrices M and K and vector f are given by(12a)M=⋃e∫ΩeNTρNdΩ(12b)K=⋃e∫ΩeBTDBdΩ(13)fω=⋃e∫ΩeNTbx,ωdΩ+∫ΓτeNTτ-x,ωdΓandN,B, andD stand for the interpolation, the strain, and the constitutive matrices of element e, respectively.Once the system of (11) is solved (i.e., vector u is computed), the approximate solution at any point of interest can be obtained using definition (10). For more details considering the numerical modelling by the finite element method, see [28–30].


## 3. Coupling Formulation

For the coupled analysis in question, the following continuity and equilibrium equations must hold at the interfaces between the acoustic fluid and solid subdomains:(14a)iωuNx,ω+vx,ω=0(14b)τNx,ω+px,ω=0where uN and τN stand for normal (normal to the common interface) displacements and tractions, respectively.The coupling between the acoustic fluid (MFS) and the elastic solid (FEM) subdomains of the model is enabled by implementing an iterative procedure that performs a successive renewal of the relevant variables at the acoustic-elastodynamic interface. The proposed approach is based on the imposition of prescribed velocities and tractions at the fluid and solid common interfaces, respectively, following relations (14a) and (14b). Since the subdomains are analyzed separately, independent discretizations can easily be considered at each subdomain without requiring matching nodes on common interfaces, thereby enhancing the flexibility of the technique. This is especially important when remeshing is considered and adaptive refinement is implemented, as it is the case here. To ensure and/or to speed up convergence, a relaxation parameter λ is introduced in the iterative coupling algorithm. The effectiveness of the iterative process depends to a great extent on the selection of this relaxation parameter, since an inappropriate choice for λ can significantly increase the number of iterations in the analysis or, even worse, make convergence unfeasible.In the subsections that follow, first, the basic steps of the iterative coupling procedure are described. Afterwards, we give an expression for an optimal relaxation parameter and then discuss the introduction of the adaptive refinement.


### 3.1. Iterative Coupling

In thekth iterative step of the MFS-FEM acoustic-elastodynamic coupling, the MFS subdomain is analyzed and vector αk is computed, as described in Section 2.1 (the superscript indicates the iterative step of the analysis). In this case, the MFS analysis takes into account prescribed velocities at the collocation points on the common interfaces, which are given by the previous iterative step (in the first iterative step, null prescribed velocities are assumed). Once αk is computed, the fluid pressures pk at the FEM nodes of the common interfaces are evaluated, using (3). These pressures are then employed to compute the prescribed tractions at the FEM interface, following (14b) (i.e., τi=-nip):(15)τ-ikx,ω=-nix∑n=1nεNnxpnkωHere, an interpolation approach, analogous to that described by (10), is used to describe the pressure distribution along the FEM interface, as indicated by the r.h.s. of (15), where nε describes the number of nodes at the interface of the element.Once the prescribed values at the common interfaces of the FEM are known, the FEM subdomains can be analyzed so that the displacement vectoruk can be computed as described in Section 2.2. Once uk is computed, the velocities vk+λ at the MFS collocation points of the common interface can be evaluated, following approximation (10) and (14a) (i.e., v=-iωujnj):(16)v-k+λx,ω=-iω∑n=1nεNnxujknωnjxAs previously explained, we have used relaxation parameters to ensure and/or to speed up the convergence of the iterative process. Thus, the prescribed velocities that are employed by the MFS subdomains in the next iterative step are computed as follows:(17)v-mk+1ω=λv-mk+λω+1-λv-mkωwhere λ stands for the relaxation parameter. Once the prescribed values at the common interfaces of the MFS are known, the algorithm goes on to the next iterative step, repeating all the procedures described above until convergence is achieved.


### 3.2. Optimal Relaxation Parameter

In order to evaluate an optimal relaxation parameter, the following square error functional is minimized here:(18)fλ=v-k+1ω-v-kω2where v- stands for the MFS prescribed values at the solid-fluid interfaces.Taking into account the relaxation of the prescribed values for the (k+1) and (k) iterations, (19a) and (19b) may be written, based on the definition in (17):(19a)v-k+1ω=λv-k+λω+1-λv-kω(19b)v-kω=λv-k+λ-1ω+1-λv-k-1ωSubstituting (19a) and (19b) into (18) yields(20)fλ=λWk+λω+1-λWkω2=λ2Wk+λω2+2λ1-λWk+λω,Wkω+1-λ2Wkω2where the inner product definition is employed (e.g., (W,W)=W2) and new variables, as defined in (21), are considered.(21)Wk+λω=v-k+λω-v-k+λ-1ωTo find the optimalλ that minimizes the functional f(λ), (20) is differentiated with respect to λ and the result is set to zero, as described below:(22)λWk+λω2+1-2λWk+λω,Wkω+λ-1Wkω2=0Rearranging the terms in (22) yields(23)λ=Wkω,Wkω-Wk+λωWkω-Wk+λω2which is an easy to implement expression that provides an optimal value for the relaxation parameter λ, at each iterative step. It should be noted that other alternatives for the calculation of the relaxation parameter can be found in the literature, such as in [31–33]; however, this expression leads to lower computational costs than some of these alternatives (see, for instance, [33]).Additionally, one should keep in mind that the computed relaxation parameter is a complex number since the problem is formulated in the frequency domain. This complex number computation could be ranged (e.g., by imposingλ≤1), but we have observed that faster convergence is usually achieved in the iterative process if a nonrestricted relaxation parameter selection, provided by (23), is considered. Moreover, although we found that the iterative process is relatively insensitive to the value of the relaxation parameter used for the first step, a real value of λ=0.5 is assumed in all the cases discussed here.


### 3.3. Adaptive Discretization

This work uses triangular finite elements since discretizations considering this type of element are easier to adaptively refine. Adaptive discretizations are widely used nowadays and a typical loop for an adaptive FEM through local refinement basically involves 4 steps, namely: (i) solve; (ii) estimate; (iii) mark; (iv) refine/coarsen. Thus, the FEM subdomains are analyzed and the solution in the current triangular mesh obtained (first step). The error is then estimated using the computed solution (second step) and it is used to mark (third step) a set of triangles that are to be refined, which is done while keeping the triangulation shape regularity and conformity (fourth step). The adaptive procedure implemented here is based on the package provided by Chen and Zhang [34].We implemented two iterative approaches with respect to the adaptive coupling formulation. The first approach used a multiple iterative algorithm. In this case, the entire coupling iterative loop is carried out within each iterative step of the adaptive iterative loop. This approach is expected to be very computationally demanding because engaged iterative procedures are employed. To avoid this computational cost, we propose a second approach, which is expected to be more efficient. In this alternative formulation, a unified single iterative loop is considered and the coupling and adaptive analyses are carried out together in the same iterative step. As illustrated in the next section, this seems to be a good approach because the number of iterative steps required by the iterative coupling procedure is not significantly increased by handling the two iterating procedures together in the same iterative loop.Figure1 shows sketches for the multiple and the single iterative algorithms, taking into account the application being studied.Figure 1
Sketches for the iterative analysis: (a) multiple and (b) single iterative algorithms.
(a) (b)


## 3.1. Iterative Coupling

In thekth iterative step of the MFS-FEM acoustic-elastodynamic coupling, the MFS subdomain is analyzed and vector αk is computed, as described in Section 2.1 (the superscript indicates the iterative step of the analysis). In this case, the MFS analysis takes into account prescribed velocities at the collocation points on the common interfaces, which are given by the previous iterative step (in the first iterative step, null prescribed velocities are assumed). Once αk is computed, the fluid pressures pk at the FEM nodes of the common interfaces are evaluated, using (3). These pressures are then employed to compute the prescribed tractions at the FEM interface, following (14b) (i.e., τi=-nip):(15)τ-ikx,ω=-nix∑n=1nεNnxpnkωHere, an interpolation approach, analogous to that described by (10), is used to describe the pressure distribution along the FEM interface, as indicated by the r.h.s. of (15), where nε describes the number of nodes at the interface of the element.Once the prescribed values at the common interfaces of the FEM are known, the FEM subdomains can be analyzed so that the displacement vectoruk can be computed as described in Section 2.2. Once uk is computed, the velocities vk+λ at the MFS collocation points of the common interface can be evaluated, following approximation (10) and (14a) (i.e., v=-iωujnj):(16)v-k+λx,ω=-iω∑n=1nεNnxujknωnjxAs previously explained, we have used relaxation parameters to ensure and/or to speed up the convergence of the iterative process. Thus, the prescribed velocities that are employed by the MFS subdomains in the next iterative step are computed as follows:(17)v-mk+1ω=λv-mk+λω+1-λv-mkωwhere λ stands for the relaxation parameter. Once the prescribed values at the common interfaces of the MFS are known, the algorithm goes on to the next iterative step, repeating all the procedures described above until convergence is achieved.


## 3.2. Optimal Relaxation Parameter

In order to evaluate an optimal relaxation parameter, the following square error functional is minimized here:(18)fλ=v-k+1ω-v-kω2where v- stands for the MFS prescribed values at the solid-fluid interfaces.Taking into account the relaxation of the prescribed values for the (k+1) and (k) iterations, (19a) and (19b) may be written, based on the definition in (17):(19a)v-k+1ω=λv-k+λω+1-λv-kω(19b)v-kω=λv-k+λ-1ω+1-λv-k-1ωSubstituting (19a) and (19b) into (18) yields(20)fλ=λWk+λω+1-λWkω2=λ2Wk+λω2+2λ1-λWk+λω,Wkω+1-λ2Wkω2where the inner product definition is employed (e.g., (W,W)=W2) and new variables, as defined in (21), are considered.(21)Wk+λω=v-k+λω-v-k+λ-1ωTo find the optimalλ that minimizes the functional f(λ), (20) is differentiated with respect to λ and the result is set to zero, as described below:(22)λWk+λω2+1-2λWk+λω,Wkω+λ-1Wkω2=0Rearranging the terms in (22) yields(23)λ=Wkω,Wkω-Wk+λωWkω-Wk+λω2which is an easy to implement expression that provides an optimal value for the relaxation parameter λ, at each iterative step. It should be noted that other alternatives for the calculation of the relaxation parameter can be found in the literature, such as in [31–33]; however, this expression leads to lower computational costs than some of these alternatives (see, for instance, [33]).Additionally, one should keep in mind that the computed relaxation parameter is a complex number since the problem is formulated in the frequency domain. This complex number computation could be ranged (e.g., by imposingλ≤1), but we have observed that faster convergence is usually achieved in the iterative process if a nonrestricted relaxation parameter selection, provided by (23), is considered. Moreover, although we found that the iterative process is relatively insensitive to the value of the relaxation parameter used for the first step, a real value of λ=0.5 is assumed in all the cases discussed here.


## 3.3. Adaptive Discretization

This work uses triangular finite elements since discretizations considering this type of element are easier to adaptively refine. Adaptive discretizations are widely used nowadays and a typical loop for an adaptive FEM through local refinement basically involves 4 steps, namely: (i) solve; (ii) estimate; (iii) mark; (iv) refine/coarsen. Thus, the FEM subdomains are analyzed and the solution in the current triangular mesh obtained (first step). The error is then estimated using the computed solution (second step) and it is used to mark (third step) a set of triangles that are to be refined, which is done while keeping the triangulation shape regularity and conformity (fourth step). The adaptive procedure implemented here is based on the package provided by Chen and Zhang [34].We implemented two iterative approaches with respect to the adaptive coupling formulation. The first approach used a multiple iterative algorithm. In this case, the entire coupling iterative loop is carried out within each iterative step of the adaptive iterative loop. This approach is expected to be very computationally demanding because engaged iterative procedures are employed. To avoid this computational cost, we propose a second approach, which is expected to be more efficient. In this alternative formulation, a unified single iterative loop is considered and the coupling and adaptive analyses are carried out together in the same iterative step. As illustrated in the next section, this seems to be a good approach because the number of iterative steps required by the iterative coupling procedure is not significantly increased by handling the two iterating procedures together in the same iterative loop.Figure1 shows sketches for the multiple and the single iterative algorithms, taking into account the application being studied.Figure 1
Sketches for the iterative analysis: (a) multiple and (b) single iterative algorithms.
(a) (b)


## 4. Numerical Analyses

We assessed the behaviour of the proposed adaptive coupling procedures by analyzing two test problems. The computed responses are compared to reference solutions, evaluated either analytically or numerically. The selected application examples correspond to (i) a circular solid structure, embedded within a fluid medium, excited by a load positioned either in the fluid or in the solid subdomain; and (ii) a concrete dam, with internal galleries and localized cracking, coupled to a fluid reservoir.


### 4.1. Circular Inclusion Embedded in a Fluid Medium

This first test case consists of a circular inclusion defined by a unit radius and made of one of the elastic materials listed in Table1. The fluid hosting the solid inclusion is assumed to be water, with a density of 1000kg/m3 and allowing acoustic waves to travel at 1500m/s. This system is illuminated by either an acoustic source, located in the fluid at x=-2.0 m and y=0 m, or by a vertical dynamic load, located in the solid at x=-0.5 m and y=0 m. A schematic representation of the model is displayed in Figure 2.Table 1
Material properties of the solid subdomain.
Solid S1 Solid S2 Young Modulus (GPa) 30 200  Mass Density (kg/m3) 2400 7000  Poisson Ratio 0.2 0.3Figure 2
Geometry of the first application example, considering a solid circular inclusion embedded in an infinite acoustic fluid medium.We analyzed this model with the proposed MFS-FEM strategy, making use of 40 sources and collocation points for the MFS (Figure3(a)), and an initial (coarse) mesh for the FEM, with 102 nodes and 170 elements (Figure 3(b)); refinement is performed until the mesh has at least 500 nodes. In addition, for comparison purposes a uniform fine FEM mesh with 559 nodes and 1052 elements was also considered (Figure 3(c)).Figure 3
Numerical setup for the analysis of the first application example: (a) sources and collocation points used for the MFS; (b) initial FEM mesh used for the adaptive analysis; (c) uniform FEM mesh.
(a) (b) (c)Figure4 gives a first set of results for an excitation frequency of 200 Hz and considering “source 1”, which is located in the fluid. Figures 4(a) and 4(b) show that both iterative adaptive approaches lead to very similar results and in both cases the computed pressures along the line of receivers LR (see Figure 2) perfectly match the analytical solution (see Figure 4(a)). Given the final FEM meshes obtained for each adaptive iterative strategy (Figure 4(b)), we can see that they are quite similar, with refinement occurring mostly in the same regions. Figure 4(c) describes the dynamic behaviour of the optimal relaxation parameter, depicting its computed values during the iterative analysis (single iterative algorithm). In Figure 4(d), the evolution of the relative error (calculated by ascertaining the continuity conditions at the interface) is displayed, taking into account both iterative adaptive approaches (a tight tolerance error of 10−7 is adopted in this work). We can see that, in this specific case, the two proposed iterative approaches are equally efficient requiring quite similar numbers of iterations to achieve a solution. Observing the convergence curve for the so-called multiple iterative approach, it can be seen that the required tolerance of 10−7 was reached after iteration 9; at this point, the refinement process starts and a mesh refinement is performed whenever the error is below the required tolerance. Thus, in the last 5 iterations, refinements are performed until the specified number of nodes is achieved, and only after that the algorithm stops. Figure 4(d) also presents results for the uniform mesh depicted in Figure 3(c) and a nonadaptive analysis (fixed mesh). As we can see, neither the single nor the multiple iterative approaches increase significantly the number of iterative steps in the analysis.Figure 4
Results for a source in the fluid and a frequency of 200 Hz (solid S1): (a) pressures along a line of receivers; (b) final refined meshes; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)Figure5 gives analogous results for a higher frequency of 800 Hz. In this case, some differences are observed, with the single iterative strategy reaching convergence faster, with 25 iterations, while 42 iterative steps are required for the multiple iterative approach (Figure 5(d)). It is interesting to notice that the so-called multiple iterative process starts to converge quite well, but when the first refinement of the mesh is performed, this fast convergence seems to be lost and the process needs to find the convergence path again. Once again and as expected, maintaining a constant fine FEM mesh leads to a faster convergence (17 iterations); however, this is achieved at the cost of using larger matrices right from the beginning of the analysis, and of not refining the mesh at the key regions of the problem. Comparing the results obtained for the solid inclusion, in terms of absolute displacements and final meshes (Figures 5(a) and 5(b)) very similar patterns and amplitudes are observed, indicating that both iterative approaches converge to the same solution.Figure 5
Results for a source in the fluid and a frequency of 800 Hz (solid S1): absolute displacement map in the solid for the (a) single iterative approach; and (b) multiple iterative approach; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)Although the results obtained for the frequencies of 200 Hz and 800 Hz are presented in detail in Figures4 and 5, illustrating the performance of the proposed techniques, results have also been computed for a full range of frequencies, varying from 10 Hz to 1000 Hz. In this context, Figure 6(a) describes the acoustic pressure calculated at the receiver R (see Figure 2), taking into account the single iterative approach and analytical answers. As we can see, there is an excellent match between the two responses for all the analyzed frequencies. Figure 6(b) shows the absolute errors of the computed pressures for both iterative approaches. The error curves displayed in the figure reveal that comparable error levels are registered, although with a visible advantage for the single iterative procedure. The total number of iterations for each frequency is given in Figure 6(c). It can be seen that in most cases only a small number of iterative steps are required by the iterative approaches (less than 20), with higher numbers only registered at higher frequencies (a maximum number of 55 iterations are observed in Figure 6(c), for the single iterative approach). In the case of the multiple iterative approach, when higher frequencies are considered this approach seems to require a somewhat higher number of iterative steps, revealing more difficulties to achieve convergence. The results given in Figure 6 indicate that the proposed single iterative strategy can be quite effective and amounts to a very efficient and accurate technique.Figure 6
Results for a source in the fluid for a complete frequency range from 10 Hz to 1000 Hz (solid S1): (a) pressure at a receiver located at x=3 m and y=1 m; (b) absolute errors; (c) total number of iterations for each frequency.
(a) (b) (c)Tests were also performed assuming different properties for the solid material, and Figure7 presents the computed results for a frequency of 2000 Hz, for solid S2. Similar features are observed for this second material and the good performance of the single iterative approach is once again highlighted. As observed before, when the multiple iterative approach is used, the first mesh refinement is performed after convergence of the solution is reached within the specified tolerance for the initial mesh; it is clear that the refinement of the mesh originates significant changes in the response, deviating from the initially estimated response and leading to an increase of the relative error before convergence is once again achieved. In fact, in this case the advantage of the single iterative strategy is even more noticeable since the total number of iterations is closer to that of the nonadaptive analysis and much lower than that of the multiple iterative approach (almost 100). In this case, the total computational time taken by the single iterative analysis is around 50% of the nonadaptive analysis.Figure 7
Results for a source in the fluid for a frequency of 2000 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses; (c) final refined meshes.
(a) (b) (c)A supplementary set of simulations is also performed for “source 2”, which is positioned within the solid inclusion (see Figure2). Figure 8 presents the results calculated for an excitation frequency of 800 Hz (solid S1). As can be seen, once again very similar responses are obtained at the line of receivers LR, covering all the discussed procedures (Figure 8(a)). As in previous cases, Figure 8(d) shows that the single iterative procedure seems to allow a faster convergence to the solution, with fewer iterations.Figure 8
Results for a source in the solid and a frequency of 800 Hz (solid S1): (a) pressures along a line of receivers; (b) final refined meshes; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)It is interesting to note that, as expected for this configuration, the adaptive analyses give final meshes with intense refinement around the loaded point (Figure8(b)). In fact, as further illustrated in Figure 9 where absolute displacement values are displayed, very sharp variations of the response are obtained in this region, with a very steep gradient being registered in the displacement field. From Figure 9 it can be seen that both adaptive strategies can accurately capture the sharp displacement peak generated at the loaded point, while this peak is significantly smoothed when the nonadaptive analysis is considered, thus highlighting the importance of adaptive approaches.Figure 9
Absolute displacement maps for a source in solid S1 (800 Hz): (a) single iterative approach; (b) multiple iterative approach; (c) uniform mesh.
(a) (b) (c)Finally, the results computed for 800 Hz and 2000 Hz are described in Figures10 and 11, respectively, for solid S2. Once again, as the figures show, the accuracy of both iterative approaches is similar, but the single iterative strategy is considerably more efficient than the multiple iterative approach. It is interesting to note that for the higher frequency (see Figure 11(b)), the use of the single iterative process leads to a reduction of the number of iterations, even when compared with the first stage of the multiple iterative algorithm (i.e., before the first refinement). Indeed, for that situation, the initial coarse mesh makes it difficult for the multiple iterative algorithm to reach a solution without refining the mesh, thus delaying the whole process (analogous behaviour can be observed in Figure 7(b)). This advantage is quite significant and leads to additional saving when the single iterative procedure is used.Figure 10
Results for a source in the solid and a frequency of 800 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses.
(a) (b)Figure 11
Results for a source in the solid and a frequency of 2000 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses.
(a) (b)


### 4.2. Concrete Dam-Reservoir System

This second numerical application concerns a more elaborate and realistic engineering scenario in which the acoustic-elastodynamic interaction occurs between a water reservoir and a dam structure. In this case, the geometry of the dam is based on the well-known Fontana gravity dam, described in many engineering works (see, for instance, Ingraffea [35], Newell and Wagner [36]), which has challenging modelling features. This concrete dam is depicted in Figure 12(a) and has two internal inspection galleries, and cracking is reported in the scientific literature at two opposite corners of the lower gallery, due to alkali-aggregate reaction. The presence of such cracks introduces highly localized tension concentrations (especially at the crack tips), which need to be captured accurately when using a numerical model.Figure 12
Dam-reservoir system: (a) sketch of the model; (b) discretization (circles and crosses stand for the MFS virtual sources and collocation points, respectively; dashed lines stand for the initial FEM mesh).
(a) (b)To set up this coupled model, the concrete dam is modelled using the FEM and the water reservoir, with a water depth of 20 m, is modelled using the MFS with 10 collocation points and virtual sources. In the case of the MFS, Green’s functions that account for the presence of a rigid floor and a free surface based on the summation of normal modes (see [37] for details of this function), are used. For the FEM, an initial coarse mesh (108 nodes) is defined to initiate the single iterative adaptive coupling analysis (Figure 12(b)). The solid material of the dam (concrete) has a Young's modulus of 22 GPa, a Poisson's ratio of 0.17, and a mass density of 2500 kg/m3; the material damping ratio is assumed to be 1.5%, and it is introduced in the form of a complex Young's modulus given by Ec=E(1+2iξ). For the adaptive process, successive refinements in the mesh are performed (one at each iterative step) until at least 800 nodes are obtained (given this criterion, on average 9 refinement steps are performed during the analyses in question).First, a harmonic pressure load is taken to be located in the fluid at x=-200 m and y=10 m, assuming different excitation frequencies. Figure13 illustrates the convergence and the final mesh generated for the model, for excitation frequencies of 25 Hz, 50 Hz and 100 Hz. Similar convergence patterns are registered for all three frequencies, with a more oscillatory behaviour being observed at the beginning of the iterative process (when the adaptive refinement is performed), followed by progressively decreasing error, until convergence is achieved. As expected, a larger number of iterations are required for higher frequencies since a more complex displacement pattern is then generated. The final refined meshes illustrated in Figure 13 exhibit slight differences between the three frequencies. Indeed, since the adaptive refinement deals with the specific displacement field computed at each frequency and tries to refine the mesh at critical locations, it is natural that differences occur between frequencies. However, at certain specific places, such as the crack tips and the corners of the galleries, the mesh is always refined, indicating that sharper variations of the field occur.Figure 13
Convergence of the single iterative approach and final refined mesh (load in the fluid): (a) 25 Hz; (b) 50 Hz; (c) 100 Hz.
(a) (b) (c)To assess the accuracy of the obtained results, Figure14 illustrates the (exaggerated) deformed shape (real and imaginary part) of the dam at a frequency of 50 Hz, computed using (i) the initial coarse mesh, without refinement; (ii) the adaptive refinement procedure; (iii) a fine mesh with 3024 uniformly distributed nodes. As this figure shows, results are related to configurations (ii) and (iii) match quite closely, while significant deviations are seen for configuration (i). Thus, the proposed strategy leads to accurate results and properly reproduces the displacement field of the dam.Figure 14
Deformed shapes for a frequency of 50Hz (scale factor of 4∙109): (a) real results; (b) imaginary results.
(a) (b)Calculations are also performed when the load is applied to the structure. In this case, a vertical load acting at x=3.25 m and y=37 m (middle of the crest) is considered. Figure15 illustrates the convergence curves of the iterative analyses and the final refined meshes that are obtained, for frequencies of 25 Hz, 50 Hz and 100 Hz. Again, the convergence of the single iterative approach is somewhat affected by the adaptive refinement at the beginning of the analyses before quickly converging to the solution, for the three frequencies. As expected, for the present load configuration the obtained refinement patterns are different from those depicted in Figure 13, with intense refinement occurring around the application point of the vertical load, particularly for the two lower frequencies. As described in the previous subsection, a very steep change in the displacement field occurs around this point, requiring the concentration of a large amount of reduced-size elements for proper reproduction. The algorithm also leads to intense refinement around the crack tips, where similar steep variations also occur.Figure 15
Convergence of the single iterative approach and final refined mesh (load in the solid): (a) 25 Hz; (b) 50 Hz; (c) 100 Hz.
(a) (b) (c)Finally, Figure16 shows the results computed at a receiver placed at x=-20 m and y=5 m, for a range of frequencies varying from 2 Hz to 100 Hz. In Figure 16(a), the pressures computed using the single iterative adaptive refinement and using a constant fine mesh are displayed. Very good agreement is observed between these solutions, although a visible difference is registered at 56 Hz. Since there is a strong peak in the pressure response at this frequency, it is possible that a specific dynamic behaviour is occurring, such as a natural mode of the structure; in this case, the results can become quite dependent on the mesh employed. In Figure 16(b) it is possible to see that the number of iterations necessary for the adaptive analyses is always higher than that required for the nonadaptive analyses, since the adaptive remeshing somewhat delays the initial convergence of the iterative process. However, this difference is attenuated for higher frequencies. The relative CPU times of the analyses are displayed in Figure 16(c) (i.e., the CPU time of the single iterative adaptive analysis divided by the CPU time of the uniform fine mesh analysis is displayed, as a percentage), clearly indicating the superior performance of the single iterative adaptive procedure. There is a significant gain at all frequencies for the proposed adaptive approach, with the related CPU time even reaching just 30% of that provided by the uniform mesh analysis. It can be seen that the proposed single iterative adaptive coupling algorithm is quite efficient and accurate and very useful in solving practical and complex engineering problems.Figure 16
Results throughout a frequency range from 2 Hz to 100 Hz (load in the solid): (a) pressure at x=-20 m and y=5 m; (b) number of iterations in the analyses; (c) percentage of CPU time (single iterative adaptive analysis/uniform fine mesh analysis).
(a) (b) (c)


## 4.1. Circular Inclusion Embedded in a Fluid Medium

This first test case consists of a circular inclusion defined by a unit radius and made of one of the elastic materials listed in Table1. The fluid hosting the solid inclusion is assumed to be water, with a density of 1000kg/m3 and allowing acoustic waves to travel at 1500m/s. This system is illuminated by either an acoustic source, located in the fluid at x=-2.0 m and y=0 m, or by a vertical dynamic load, located in the solid at x=-0.5 m and y=0 m. A schematic representation of the model is displayed in Figure 2.Table 1
Material properties of the solid subdomain.
Solid S1 Solid S2 Young Modulus (GPa) 30 200  Mass Density (kg/m3) 2400 7000  Poisson Ratio 0.2 0.3Figure 2
Geometry of the first application example, considering a solid circular inclusion embedded in an infinite acoustic fluid medium.We analyzed this model with the proposed MFS-FEM strategy, making use of 40 sources and collocation points for the MFS (Figure3(a)), and an initial (coarse) mesh for the FEM, with 102 nodes and 170 elements (Figure 3(b)); refinement is performed until the mesh has at least 500 nodes. In addition, for comparison purposes a uniform fine FEM mesh with 559 nodes and 1052 elements was also considered (Figure 3(c)).Figure 3
Numerical setup for the analysis of the first application example: (a) sources and collocation points used for the MFS; (b) initial FEM mesh used for the adaptive analysis; (c) uniform FEM mesh.
(a) (b) (c)Figure4 gives a first set of results for an excitation frequency of 200 Hz and considering “source 1”, which is located in the fluid. Figures 4(a) and 4(b) show that both iterative adaptive approaches lead to very similar results and in both cases the computed pressures along the line of receivers LR (see Figure 2) perfectly match the analytical solution (see Figure 4(a)). Given the final FEM meshes obtained for each adaptive iterative strategy (Figure 4(b)), we can see that they are quite similar, with refinement occurring mostly in the same regions. Figure 4(c) describes the dynamic behaviour of the optimal relaxation parameter, depicting its computed values during the iterative analysis (single iterative algorithm). In Figure 4(d), the evolution of the relative error (calculated by ascertaining the continuity conditions at the interface) is displayed, taking into account both iterative adaptive approaches (a tight tolerance error of 10−7 is adopted in this work). We can see that, in this specific case, the two proposed iterative approaches are equally efficient requiring quite similar numbers of iterations to achieve a solution. Observing the convergence curve for the so-called multiple iterative approach, it can be seen that the required tolerance of 10−7 was reached after iteration 9; at this point, the refinement process starts and a mesh refinement is performed whenever the error is below the required tolerance. Thus, in the last 5 iterations, refinements are performed until the specified number of nodes is achieved, and only after that the algorithm stops. Figure 4(d) also presents results for the uniform mesh depicted in Figure 3(c) and a nonadaptive analysis (fixed mesh). As we can see, neither the single nor the multiple iterative approaches increase significantly the number of iterative steps in the analysis.Figure 4
Results for a source in the fluid and a frequency of 200 Hz (solid S1): (a) pressures along a line of receivers; (b) final refined meshes; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)Figure5 gives analogous results for a higher frequency of 800 Hz. In this case, some differences are observed, with the single iterative strategy reaching convergence faster, with 25 iterations, while 42 iterative steps are required for the multiple iterative approach (Figure 5(d)). It is interesting to notice that the so-called multiple iterative process starts to converge quite well, but when the first refinement of the mesh is performed, this fast convergence seems to be lost and the process needs to find the convergence path again. Once again and as expected, maintaining a constant fine FEM mesh leads to a faster convergence (17 iterations); however, this is achieved at the cost of using larger matrices right from the beginning of the analysis, and of not refining the mesh at the key regions of the problem. Comparing the results obtained for the solid inclusion, in terms of absolute displacements and final meshes (Figures 5(a) and 5(b)) very similar patterns and amplitudes are observed, indicating that both iterative approaches converge to the same solution.Figure 5
Results for a source in the fluid and a frequency of 800 Hz (solid S1): absolute displacement map in the solid for the (a) single iterative approach; and (b) multiple iterative approach; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)Although the results obtained for the frequencies of 200 Hz and 800 Hz are presented in detail in Figures4 and 5, illustrating the performance of the proposed techniques, results have also been computed for a full range of frequencies, varying from 10 Hz to 1000 Hz. In this context, Figure 6(a) describes the acoustic pressure calculated at the receiver R (see Figure 2), taking into account the single iterative approach and analytical answers. As we can see, there is an excellent match between the two responses for all the analyzed frequencies. Figure 6(b) shows the absolute errors of the computed pressures for both iterative approaches. The error curves displayed in the figure reveal that comparable error levels are registered, although with a visible advantage for the single iterative procedure. The total number of iterations for each frequency is given in Figure 6(c). It can be seen that in most cases only a small number of iterative steps are required by the iterative approaches (less than 20), with higher numbers only registered at higher frequencies (a maximum number of 55 iterations are observed in Figure 6(c), for the single iterative approach). In the case of the multiple iterative approach, when higher frequencies are considered this approach seems to require a somewhat higher number of iterative steps, revealing more difficulties to achieve convergence. The results given in Figure 6 indicate that the proposed single iterative strategy can be quite effective and amounts to a very efficient and accurate technique.Figure 6
Results for a source in the fluid for a complete frequency range from 10 Hz to 1000 Hz (solid S1): (a) pressure at a receiver located at x=3 m and y=1 m; (b) absolute errors; (c) total number of iterations for each frequency.
(a) (b) (c)Tests were also performed assuming different properties for the solid material, and Figure7 presents the computed results for a frequency of 2000 Hz, for solid S2. Similar features are observed for this second material and the good performance of the single iterative approach is once again highlighted. As observed before, when the multiple iterative approach is used, the first mesh refinement is performed after convergence of the solution is reached within the specified tolerance for the initial mesh; it is clear that the refinement of the mesh originates significant changes in the response, deviating from the initially estimated response and leading to an increase of the relative error before convergence is once again achieved. In fact, in this case the advantage of the single iterative strategy is even more noticeable since the total number of iterations is closer to that of the nonadaptive analysis and much lower than that of the multiple iterative approach (almost 100). In this case, the total computational time taken by the single iterative analysis is around 50% of the nonadaptive analysis.Figure 7
Results for a source in the fluid for a frequency of 2000 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses; (c) final refined meshes.
(a) (b) (c)A supplementary set of simulations is also performed for “source 2”, which is positioned within the solid inclusion (see Figure2). Figure 8 presents the results calculated for an excitation frequency of 800 Hz (solid S1). As can be seen, once again very similar responses are obtained at the line of receivers LR, covering all the discussed procedures (Figure 8(a)). As in previous cases, Figure 8(d) shows that the single iterative procedure seems to allow a faster convergence to the solution, with fewer iterations.Figure 8
Results for a source in the solid and a frequency of 800 Hz (solid S1): (a) pressures along a line of receivers; (b) final refined meshes; (c) evolution of the relaxation parameter (single iterative algorithm); (d) convergence of the iterative analyses.
(a) (b) (c) (d)It is interesting to note that, as expected for this configuration, the adaptive analyses give final meshes with intense refinement around the loaded point (Figure8(b)). In fact, as further illustrated in Figure 9 where absolute displacement values are displayed, very sharp variations of the response are obtained in this region, with a very steep gradient being registered in the displacement field. From Figure 9 it can be seen that both adaptive strategies can accurately capture the sharp displacement peak generated at the loaded point, while this peak is significantly smoothed when the nonadaptive analysis is considered, thus highlighting the importance of adaptive approaches.Figure 9
Absolute displacement maps for a source in solid S1 (800 Hz): (a) single iterative approach; (b) multiple iterative approach; (c) uniform mesh.
(a) (b) (c)Finally, the results computed for 800 Hz and 2000 Hz are described in Figures10 and 11, respectively, for solid S2. Once again, as the figures show, the accuracy of both iterative approaches is similar, but the single iterative strategy is considerably more efficient than the multiple iterative approach. It is interesting to note that for the higher frequency (see Figure 11(b)), the use of the single iterative process leads to a reduction of the number of iterations, even when compared with the first stage of the multiple iterative algorithm (i.e., before the first refinement). Indeed, for that situation, the initial coarse mesh makes it difficult for the multiple iterative algorithm to reach a solution without refining the mesh, thus delaying the whole process (analogous behaviour can be observed in Figure 7(b)). This advantage is quite significant and leads to additional saving when the single iterative procedure is used.Figure 10
Results for a source in the solid and a frequency of 800 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses.
(a) (b)Figure 11
Results for a source in the solid and a frequency of 2000 Hz (solid S2): (a) pressures along a line of receivers; (b) convergence of the iterative analyses.
(a) (b)


## 4.2. Concrete Dam-Reservoir System

This second numerical application concerns a more elaborate and realistic engineering scenario in which the acoustic-elastodynamic interaction occurs between a water reservoir and a dam structure. In this case, the geometry of the dam is based on the well-known Fontana gravity dam, described in many engineering works (see, for instance, Ingraffea [35], Newell and Wagner [36]), which has challenging modelling features. This concrete dam is depicted in Figure 12(a) and has two internal inspection galleries, and cracking is reported in the scientific literature at two opposite corners of the lower gallery, due to alkali-aggregate reaction. The presence of such cracks introduces highly localized tension concentrations (especially at the crack tips), which need to be captured accurately when using a numerical model.Figure 12
Dam-reservoir system: (a) sketch of the model; (b) discretization (circles and crosses stand for the MFS virtual sources and collocation points, respectively; dashed lines stand for the initial FEM mesh).
(a) (b)To set up this coupled model, the concrete dam is modelled using the FEM and the water reservoir, with a water depth of 20 m, is modelled using the MFS with 10 collocation points and virtual sources. In the case of the MFS, Green’s functions that account for the presence of a rigid floor and a free surface based on the summation of normal modes (see [37] for details of this function), are used. For the FEM, an initial coarse mesh (108 nodes) is defined to initiate the single iterative adaptive coupling analysis (Figure 12(b)). The solid material of the dam (concrete) has a Young's modulus of 22 GPa, a Poisson's ratio of 0.17, and a mass density of 2500 kg/m3; the material damping ratio is assumed to be 1.5%, and it is introduced in the form of a complex Young's modulus given by Ec=E(1+2iξ). For the adaptive process, successive refinements in the mesh are performed (one at each iterative step) until at least 800 nodes are obtained (given this criterion, on average 9 refinement steps are performed during the analyses in question).First, a harmonic pressure load is taken to be located in the fluid at x=-200 m and y=10 m, assuming different excitation frequencies. Figure13 illustrates the convergence and the final mesh generated for the model, for excitation frequencies of 25 Hz, 50 Hz and 100 Hz. Similar convergence patterns are registered for all three frequencies, with a more oscillatory behaviour being observed at the beginning of the iterative process (when the adaptive refinement is performed), followed by progressively decreasing error, until convergence is achieved. As expected, a larger number of iterations are required for higher frequencies since a more complex displacement pattern is then generated. The final refined meshes illustrated in Figure 13 exhibit slight differences between the three frequencies. Indeed, since the adaptive refinement deals with the specific displacement field computed at each frequency and tries to refine the mesh at critical locations, it is natural that differences occur between frequencies. However, at certain specific places, such as the crack tips and the corners of the galleries, the mesh is always refined, indicating that sharper variations of the field occur.Figure 13
Convergence of the single iterative approach and final refined mesh (load in the fluid): (a) 25 Hz; (b) 50 Hz; (c) 100 Hz.
(a) (b) (c)To assess the accuracy of the obtained results, Figure14 illustrates the (exaggerated) deformed shape (real and imaginary part) of the dam at a frequency of 50 Hz, computed using (i) the initial coarse mesh, without refinement; (ii) the adaptive refinement procedure; (iii) a fine mesh with 3024 uniformly distributed nodes. As this figure shows, results are related to configurations (ii) and (iii) match quite closely, while significant deviations are seen for configuration (i). Thus, the proposed strategy leads to accurate results and properly reproduces the displacement field of the dam.Figure 14
Deformed shapes for a frequency of 50Hz (scale factor of 4∙109): (a) real results; (b) imaginary results.
(a) (b)Calculations are also performed when the load is applied to the structure. In this case, a vertical load acting at x=3.25 m and y=37 m (middle of the crest) is considered. Figure15 illustrates the convergence curves of the iterative analyses and the final refined meshes that are obtained, for frequencies of 25 Hz, 50 Hz and 100 Hz. Again, the convergence of the single iterative approach is somewhat affected by the adaptive refinement at the beginning of the analyses before quickly converging to the solution, for the three frequencies. As expected, for the present load configuration the obtained refinement patterns are different from those depicted in Figure 13, with intense refinement occurring around the application point of the vertical load, particularly for the two lower frequencies. As described in the previous subsection, a very steep change in the displacement field occurs around this point, requiring the concentration of a large amount of reduced-size elements for proper reproduction. The algorithm also leads to intense refinement around the crack tips, where similar steep variations also occur.Figure 15
Convergence of the single iterative approach and final refined mesh (load in the solid): (a) 25 Hz; (b) 50 Hz; (c) 100 Hz.
(a) (b) (c)Finally, Figure16 shows the results computed at a receiver placed at x=-20 m and y=5 m, for a range of frequencies varying from 2 Hz to 100 Hz. In Figure 16(a), the pressures computed using the single iterative adaptive refinement and using a constant fine mesh are displayed. Very good agreement is observed between these solutions, although a visible difference is registered at 56 Hz. Since there is a strong peak in the pressure response at this frequency, it is possible that a specific dynamic behaviour is occurring, such as a natural mode of the structure; in this case, the results can become quite dependent on the mesh employed. In Figure 16(b) it is possible to see that the number of iterations necessary for the adaptive analyses is always higher than that required for the nonadaptive analyses, since the adaptive remeshing somewhat delays the initial convergence of the iterative process. However, this difference is attenuated for higher frequencies. The relative CPU times of the analyses are displayed in Figure 16(c) (i.e., the CPU time of the single iterative adaptive analysis divided by the CPU time of the uniform fine mesh analysis is displayed, as a percentage), clearly indicating the superior performance of the single iterative adaptive procedure. There is a significant gain at all frequencies for the proposed adaptive approach, with the related CPU time even reaching just 30% of that provided by the uniform mesh analysis. It can be seen that the proposed single iterative adaptive coupling algorithm is quite efficient and accurate and very useful in solving practical and complex engineering problems.Figure 16
Results throughout a frequency range from 2 Hz to 100 Hz (load in the solid): (a) pressure at x=-20 m and y=5 m; (b) number of iterations in the analyses; (c) percentage of CPU time (single iterative adaptive analysis/uniform fine mesh analysis).
(a) (b) (c)


## 5. Conclusions

This work has described a new strategy for addressing acoustic-elastodynamic interaction problems in the frequency domain, using optimized iterative coupling between two distinct numerical methods and incorporating adaptive refinement within the iterative process. The use of iterative coupling procedures is known to be quite efficient in such problems because it avoids complex direct coupling matrices and makes it simpler to exploit the individual advantages of each of the coupled numerical methods. The MFS (for the acoustic fluid) and the FEM (for the solid) are used in the proposed method, the first being a boundary meshless method, and the latter being a traditional and well established domain discretization method. The adaptive approach is applied to the FEM to make it possible to refine its mesh based on the progressive solutions that are computed during the iterative coupling analysis. The whole strategy is simple to implement and can be used to link the existing codes developed for each method without having to make significant changes. Indeed, the coupling procedure only makes use of the interface results provided by each individual method, while the adaptive refinement is applied to the FEM (solid) and does not interfere with the MFS (acoustic fluid) part of the model. The adaptive refinement could equally be applied to the fluid part of the model without interfering with the elastodynamic subdomain. The examples explored in this paper demonstrate that the proposed strategy is accurate and efficient and has a great potential for practical applications, in that it provides reliable results even when more complex structures are analyzed. It was noted that although the adaptive refinement caused a small delay to the convergence process, the gain of concentrating the new elements where they are indeed required is a significant advantage, which is even more important when working with dynamic (frequency dependent) problems.



---

*Source: 4834521-2019-01-06.xml*