Generation of smooth grids with line control for scattering from multiple obstacles

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 79 (2009) 2506–2520

Generation of smooth grids with line control for scattering from multiple obstacles Vianey Villamizar a,∗ , Sebastian Acosta b b

a Department of Mathematics, Brigham Young University, Provo, UT 84602, USA Department of Mechanical Engineering and Department of Mathematics, Brigham Young University, Provo, UT 84602, USA

Available online 31 January 2009

Abstract A new approach to generate structured grids for two-dimensional multiply connected regions with several holes is proposed. The bounding curves may include corners or cusps. The new algorithm constitutes an extension of the Branch Cut Grid Line Control (BCGC) technique introduced byVillamizar et al. [V. Villamizar, O. Rojas, J. Mabey, Generation of curvilinear coordinates on multiply connected regions with boundary-singularities, J. Comput. Phys. 223 (2007) 571–588] to domains with a finite number of holes. Regions with multiple holes are reduced to several contiguous single hole subregions. Then, the BCGC algorithm is applied to each single hole subregion producing a smooth grid with line control. Finally, the subregions with their respective grids are joined and their interfaces are smoothed resulting a globally smooth grid. The advantages of the novel grids are revealed by employing them to numerically solve acoustic scattering problems in the presence of multiple complexly shaped obstacles. © 2009 IMACS. Published by Elsevier B.V. All rights reserved. PACS: 65N50; 74J20; 74S20 Keywords: Elliptic grid generation; Grid line control; Smoothing process; Multiple scattering; Radar cross-section

1. Introduction There are numerous important problems in physical applications that require the development of numerical methods over complex geometries with several holes. For example, viscous flow over arbitrarily shaped bodies [8], or, acoustic scattering from multiple complexly shaped obstacles [9] among many others. These problems naturally lead to grid generation over complex geometries. Our focus in this work is on the construction of smooth, non-self-overlapping grids with line control over geometrically complex domains and its application to wave scattering. Generation of structured grids has been the subject of intensive research for more than 30 years. One of the most popular approaches is based on the numerical solution of quasi-linear elliptic systems of partial differential equations. An important attribute of grids obtained from elliptic grid generators is their smoothness. As a result, slope discontinuities on the boundaries are not propagated into the field [5]. Dealing with complexly shaped regions constitutes a challenge for any grid generator. For elliptic grid generators, the common practice consists of dividing the physical region into subregions and generating a structured grid called a subgrid. Neighboring subregions form internal interfaces across which information must be transferred [10]. Smoothness across the interfaces needs a special treatment of the interface ∗

Corresponding author. E-mail addresses: [email protected] (V. Villamizar), [email protected] (S. Acosta).

0378-4754/$36.00 © 2009 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2009.01.006

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2507

nodes. It requires an additional data indexing procedure to link the blocks across the interfaces throughout the generation process. The smoothness of the grids happens to be a crucial property in numerical simulation. Accurate numerical modeling of physical phenomena requires smooth grids in general. For instance in [6], the authors reported a deficit in the computation of the pressure field produced by the acoustic scattering from multiple rigid circular cylinders. They attributed the deficit to the use of the standard multi-block grid generation strategy which, for some geometries, cannot avoid discontinuities of grid lines through the boundaries between adjacent blocks. In [12], a domain with a single hole bounded by an astroid-type curve was treated. A grid with slope discontinuity at certain locations was generated for this domain using the CAE package ANSYS. The transverse displacement of a membrane with this shape and fixed at its boundary was numerically modelled employing this grid. Unstable oscillations were observed after a moderate number of iterations. Grid line control is also a highly desired property in numerical simulation [10,1]. A possible approach to establish grid line control consists of defining appropriate control functions forcing the elliptic system. These functions are called grid control functions. In [12], a practical method to define the control functions for two-dimensional singlehole multiply connected regions was proposed. This method only depends upon the initial distribution of nodes on the branch cut (used to define a mapping from a rectangular computational domain to a physical domain) and the physical boundaries. A novel smoothing mechanism that depends on branch cut interchanges for successive iterations was also developed. In this work, the boundary-conforming coordinates generation algorithm with grid control (BCGC), developed in [12], is improved and extended. Here, the purpose is to introduce a new formulation of the BCGC that can handle twodimensional multiply connected regions with multiple holes while maintaining grid line control inside the region. To illustrate the construction of the algorithm, a prototypical multiply connected domain with two holes has been chosen (see Fig. 1). The treatment of this region consists of dividing it into two regions or blocks with a single hole each. Then, the BCGC procedure is applied to each individual block, the blocks are rejoined, and a smoothing mechanism is applied. More general multiply connected domains with several holes are treated in a similar manner by reducing them to several single hole regions. To measure the grid smoothness at each node, two new parameters “maximum deviation from smoothness (MDS)” and “average deviation from smoothness (ADS)” are defined in Section 3. In our prototypical domain depicted in Fig. 1, the holes are bounded by the three-leafed rose and the astroid closed curves which were originally defined in [12]. The astroid curve includes severe boundary singularities which are difficult to handle for most numerical grid generators. In what follows, we describe how to generate smooth grids with a desired mesh spacing for complex multiply connected regions with multiple holes even in the presence of singular bounding curves.

Fig. 1. Multiply connected domain with two complexly shaped holes.

2508

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2. BCGC extended to multiply connected regions with several holes The BCGC method was initially introduced for single hole multiply connected regions in [12]. It is based on the numerical solution of an elliptic boundary value problem. The method requires a definition of a transformation T from the computational domain D with rectangular coordinates (ξ, η) (1 ≤ ξ ≤ N1 and 1 ≤ η ≤ N2 ), to the multiply connected physical domain D with curvilinear coordinates (x(ξ, η), y(ξ, η)) (see Fig. 1). The transformation T is defined by numerically solving a Dirichlet boundary value problem governed by the familiar quasi-linear elliptic system of partial differential equations [11,12] given by αxξξ − 2βxξη + γxηη = −αψ(ξ, η) xξ − γφ(ξ, η) xη ,

(1)

αyξξ − 2βyξη + γyηη = −αψ(ξ, η) yξ − γφ(ξ, η) yη .

(2)

The symbols α, β and γ represent the scale metric factors of the transformation T. They are defined as α = xη2 + yη2 , β = xξ xη + yξ yη , γ = xξ2 + yξ2 . The control functions ψ and φ are defined from the initial distribution of the grid points at the physical boundaries and at the branch cut used to define the transformation T. In fact, the control function φ is initially defined along the branch cut used to define the transformation T (vertical rectangular boundaries of the computational domain) as φ(ξ, η) =

−(xη xηη + yη yηη ) , xη 2 + y η 2

when ξ = 1 and ξ = N1 .

(3)

Also, the control function ψ is initially defined along the horizontal rectangular boundaries of the computational domain (physical and fictitious infinite boundaries) as ψ(ξ, η) =

−(xξ xξξ + yξ yξξ ) , xξ 2 + y ξ 2

when η = 1 and η = N2 .

(4)

It is noticed that the control functions ψ and φ involve derivatives with respect to only one variable. For instance, the control function ψ only contains derivatives with respect to ξ, which is also the free variable along the horizontal straight segments η = 1 and η = N2 corresponding to the physical and fictitious infinite boundaries, respectively. Similarly, the control function φ only contains derivatives with respect to η, which is the only free variable along the vertical straight segments ξ = 1 and ξ = N1 corresponding to a branch cut C. Therefore, the derivative terms in the above definitions for ψ and φ can be easily approximated from the initial distribution of grid points on the physical boundaries, fictitious infinite boundary, and the branch cut, respectively. In fact, by using centered finite difference in the variable η the following approximations are obtained   1 (x1,j+1 − x1,j−1 )(x1,j+1 − 2x1,j + x1,j−1 ) (y1,j+1 − y1,j−1 )(y1,j+1 − 2y1,j + y1,j−1 ) φ1,j = − + , α1,j 2 2 φN1 ,j = φ1,j ,

j = 2, . . . N2 − 1.

(5)

If j = N2 then, second order backward finite-difference approximations of the derivative terms present in φ(1, N2 ) are employed. Analogously, values of ψ at the physical and fictitious infinite boundaries can be easily approximated from an initial distribution of points using centered differences in the variable ξ. Finally, the values of φ(ξ, η) and ψ(ξ, η) at the interior points of the computational domain D are defined by linear interpolation along the η-curves (η = constant) and the ξ-curves (ξ = constant), respectively. Numerical experiments on several multiply connected domains with a single hole in [12,13] revealed that the desired spacing among the ξ-curves and the η-curves of the final grid is governed by these control functions. Since the control functions are defined from the initial distribution of the grid points at the physical boundaries and at the branch cut, it is this initial distribution of points which ultimately governs the spacing among grid lines. Details about the definition of the initial distribution of points are offered below. Now that more holes are considered, the goal is to have grid line control in the vicinity of all of them. For multiple holes, it is necessary to divide the region into blocks. For the region shown in Fig. 1, the domain is decomposed into two blocks. They are denoted by Bk , where k = 1, 2. The two blocks are separated by two different η-curves (annular cut in Fig. 1), ACl and ACr , that share a common interface. These η-curves are called annular cuts. The importance of this procedure is that generation of boundary-conforming grids with line control for regions with multiple holes can be easily reduced to an application of the BCGC algorithm to regions with a single hole.

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2509

Each block is treated independently as described in [12]. To begin, a uniformly spaced rectangular grid with step sizes ξ = 1 and η = 1 is defined on the computational domain D . The discrete values for ξ and η are represented by ξi = iΔξ and ηj = jΔη, for i = 1, . . . Nk,1 and j = 1, . . . Nk,2 , respectively. The approximation of the curvilinear coordinates x and y at the grid point (ξi , ηj ) are denoted by xi,j and yi,j , respectively. A consideration of Fig. 1 helps to understand the notation and some definitions. The branch cuts BCk for each block Bk (k = 1, 2) are chosen such that they are coincident with an ξ-curve. The length of these branch cuts are lk and the number of grid points in each branch cut along the η-direction is Nk,2 , then a uniform step size hk = lk /(Nk,2 − 1) can be defined. The control parameters for the block Bk are given by fkl × nkl /fkr × nkr . It means that there is a cluster consisting of nkl grid points on BCk in the vicinity of the most inner end separated by a uniform distance hkl = fkl hk . Similarly, the second half of the above expression means that there is another cluster of nkr grid points on BCk in the vicinity of the most outer end separated by a uniform distance hkr = fkr hk . The scaling factors fkl and fkr produce a contraction of the distance between the grid points forming the clusters (nl and nr ) if they are less than one, or a stretching, if their values are greater than one. The distribution of the intermediate points along the branch cut of each block is determined by a fifth degree osculating polynomial interpolating the position of grid points on each end of the branch cut. The use of a fifth degree osculating polynomial guarantees a smooth cell size transition between the cluster cells located at the two ends and those cells immediately adjacent located in the intermediate region. The osculating polynomial can be easily obtained by applying divided-difference interpolation methods. It can be written as x(ξ1 , ηj ) = a0 + (ηj − ηnl ) hl + a3 (ηj − ηnl )3 + a4 (ηj − ηnl )3 (ηj − ηj−N2 +nr −1 ) +a5 (ηj − ηnl )3 (ηj − ηj−N2 +nr −1 )2 , y(ξ1 , ηj ) = 0,

(6)

j = nl + 1 . . . N2 − nr ,

where a0 = r0 + (nl − 1) hl , r0 = aux1 = r∞ − (nr − 1) hr .

(7)

 (x(1, 1)2 + y(1, 1)2 ),

Here, r0 is the distance from the origin to the first point along the branch cut. The symbol r∞ represent the distance from the origin to the farthest point along the branch cut. Also, aux1 − a0 , aux2 hr − 2 aux3 + hl aux4 − a3 aux4 = , a4 = , aux2 aux22 aux3 − hr aux5 − 2 aux4 + a3 aux5 = , a5 = . aux22 aux22 aux2 = N2 − (nr + nl − 1),

aux3 =

a3 =

aux3 − hl , aux22

From the initial node distributions over the branch cuts and the physical boundaries the control functions φ and ψ are defined. Then, an initial algebraic grid for the physical domain D that includes the nodes present in the initial distribution is constructed. The BCGC algorithm is applied to the blocks B1 and B2 shown in Fig. 2. These two blocks ˜ for the domain with their recent constructed grids are joined along the ACl and ACr η-curves resulting a new grid G ˜ D, as seen at the left portion of Fig. 3. Since the blocks are treated independently, the grid G is not necessarily smooth across the block interface. An additional smoothing mechanism that consists of transforming the region D with the ˜ into a simply connected block S1 is introduced. The block S1 is obtained by removing a portion of the region grid G D as indicated at the right of Fig. 3. It should be noticed that the non-smooth interface is contained in the interior of S1 . Therefore, an application of the BCGC algorithm to S1 corrects the non-smoothness along the interface between blocks B1 and B2 . The partial process just described, consisting of the application of the BCGC algorithm to the multiply connected blocks B1 , B2 , and to the simply connected one S1 , will be called a cycle. A new cycle is started if the total number of SOR iterations within the previous cycle exceeds 10. This constitutes the stop criteria of the algorithm. In Table 1, the

2510

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

Fig. 2. Block B1 BCGC grid (left) and block B2 BCGC grid (right)

BCGC grid generation convergence process for a region with two holes is recorded. A step by step description of the BCGC process applied to the physical domain D with two holes in its interior is the subject of the following section. The generalization to regions with more than two holes is straightforward. The prototypical region showed in Figs. 2–4 has branch cut lengths: l1 = 5 and l2 = 5. The grid is formed by N1 = 81 ξ-curves and N2 = 31 η-curves in total. The region is divided into two blocks B1 of size N1 × N2 = 81 × 16; and B2 of size N1 × N2 = 81 × 15. The control parameters are 1.5 × 4/0.5 × 4 for both blocks. The BCGC grids corresponding to the two blocks B1 and B2 are depicted in Fig. 2. As discussed above, the grid formed by joining the two blocks is non-smooth across the interface (see the left part of Fig. 3). The simply connected region used to smooth

Fig. 3. Non-smooth grid at interface (left). Smoothed simply connected block S1 (right).

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2511

Table 1 BCGC grid convergence with grid control parameters: 1.5 × 4/0.5 × 4/1.5 × 4/0.5 × 4. Tol = 1e−4. Grid size

81 × 81

101 × 101

121 × 121

161 × 161

201 × 201

B1 B2 S1

132 89 83

162 102 80

196 167 75

282 191 101

335 251 80

G0 − G1 G1 − G2 G2 − G3 G3 − G4 G4 − G5 G5 − G6 G6 − G7 G7 − G8

304 124 201 90 45 9 – –

344 114 261 89 47 13 3 –

438 98 381 71 35 12 6 –

574 174 683 160 80 43 23 9

666 138 908 107 57 34 12 8

Total iter

773

871

1041

1746

1930

the interface is also shown at the right in Fig. 3. The final grid obtained after iterative application of the above process is illustrated in Fig. 4. 2.1. The smoothing process We now, summarize the steps required to construct BCGC grids on complex regions with two holes. Extension to regions with more than two holes is straightforward. Step 1: Partition of the domain into blocks and initial distribution of nodes - Divide the domain into two blocks, B1 and B2 , with a single hole in their respective interiors, as shown in Fig. 2.This is accomplished by selecting two η-curves as annular cuts ACl and ACr that include part of the boundary of one of the holes. Then, select a branch cut BC1 and BC2 for each block, as illustrated in Fig. 1.

Fig. 4. Final BCGC grid (left). Final BCGC grid with local adjustments (right).

2512

Step 2:

Step 3:

Step 4:

Step 5:

Step 6:

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

- Specify nodes on the branch cuts according to the desired spacing among the grid lines in the final grid. This requires the introduction of stretching and/or contraction parameters fkl , fkr , and also the definition of an osculating polynomial for the node distribution, as explained in the previous section. - Create an initial mesh G0 for the physical domain D containing all the grid points initially prescribed on the branch cuts and on the physical boundaries. Control function definitions - Define the control functions φ and ψ from the initial distribution of points on the branch cuts and on the physical boundaries as described in the previous section. Application of the BCGC algorithm to each individual block - Apply the BCGC algorithm to block B1 with inner boundary given by the bounding curve of hole 1 (threeleafed rose) and outer boundary ACl . The outer boundary ACl contains the interface and the left part of the bounding curve of hole 2 (see left of Fig. 2). - Apply the BCGC algorithm to block B2 . Its inner boundary is ACr , which also contains the interface and the right part of the bounding curve of hole 2. The outer boundary of B2 consists of the most outer η-curve, η = N2 (see right of Fig. 2). The blocks with their respective grids are joined. The interface is smoothed - At this step the blocks B1 and B2 with their new grids are rejoined and the original domain D is recovered. In general, D is not smooth at the interface, as illustrated at the left of Fig. 3. Thus, a smoothing needs to take place. - A region Sh containing hole 2 (astroid) is cut out of D. It creates a new simply connected region S1 , as shown at the right of Fig. 3. The region S1 contains the common interface of the two annular cuts ACl and ACr in its interior. - The BCGC algorithm is applied to the simply connected region S1 . As a consequence, the grid is smoothed along the interface of blocks B1 and B2 . The original domain is rejoined - The two subregions Sh and S1 with their respective grids are rejoined. Therefore, the original region D is recovered with a new mesh G1 . This partial process is called a cycle. The cycle is repeated until a stop criteria is satisfied - The process (cycle) described in steps 3–5 is repeated. New meshes G2 , G3 , . . . GM are obtained. The main purpose is to generate a globally smooth grid. - The process stops when the number of iterations to compute a new grid Gi from the previous grid Gi−1 is 10 or less. The final grid is depicted at the left of Fig. 4.

The key idea of the generation technique is the use of the smoothness property of the solution for the quasi-linear elliptic system of partial differential equations defining the grid generation algorithm. When this property is properly coupled with successive interchangeable cuts (branch and annular cuts) on the physical domain, a new boundary conforming and globally smooth grid with line control is created. BCGC grids of various sizes were generated for the prototypical region shown in Fig. 1. The convergence process for each grid is reported in Table 1. The first intermediate grid G1 is obtained at the end of the first cycle. The total number of iterations to obtain G1 from the initial grid G0 are listed in the fifth row of the Table. It includes the iterations employed to smooth each individual block (B1 and B2 ), as well as the common interface S1 . The tolerance employed in this first cycle was 5 × 10−4 . The number of iterations to obtain the following intermediate grids are also reported. The last row contains the total number of iterations performed to obtain the final grid G7 or G8 from the initial grid G0 . The tolerance for the second cycle was decreased to 2 × 10−4 . For the subsequent cycles, a tolerance of 10−4 was used. This is why the number of iterations after the third cycle are more than those in the second cycle. It took less than a minute of CPU-time on a Dell Notebook D810, with a 2 GHz processor and 1 GB of RAM, to run a MATLAB code and generate each one of the grids reported in Table 1. 2.2. A posteriori local adjustment of grid cell sizes A final 81 × 31grid, with 1.5 × 4/0.5 × 4/1.5 × 4/0.5 × 4 as control parameters is shown at the left of Fig. 4. This is obtained after the application of the extended version of the BCGC algorithm, described in this work. It is noticeable

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2513

Fig. 5. Original control function φold (top). Modified control function φnew (bottom).

how the clustering at the most outer boundary and the stretching at the three-leafed rose boundary of the η-curves, imposed on the initial grid, is maintained throughout the domain D. However, the closest cells to the upper and lower cusps of the astroid are much larger than their neighbors. Indeed, the η-curves bounding these cells tend to separate more at this location than at other locations, as seen at the left of Fig. 4. The control function φ affecting the η-curve spacings is depicted at the top in Fig. 5. Our numerical experiments show that an increment of the value of φ at a given grid point produces a local displacement of the corresponding η-curves in the positive direction (outward from the origin). Similarly, a decrement of φ values brings the grid lines inward (negative direction). Moreover, the displacement of the η-curves is proportional to the magnitude of the change experienced by φ. The relationship between changes of ψ values and displacement of ξ-curves is analogous. The above observations helped us to define a posteriori local grid line control spacing mechanism for the η-curves in the final BCGC grid. In fact, by adding an appropriate quantity to the original control function φ at certain subregion H ⊂ D the η-curve spacings are modified there. More precisely, a new φ is defined by φnew (ξi , ηj ) = φold (ξi , ηj ) + j (ξi , ηj ),

for (ξi , ηj ) ∈ H.

(8)

If an outward local displacement of the η-curves is needed to improve the spacing, then a positive  is added. Otherwise, if an inward displacement of η-curves is sought, then a negative  is added. The local adjustment mechanism should be designed such that smoothness of the grid be preserved. First, the relatively large or relatively small grid cell that wants to be modified, together with its neighboring cells, are identified (target cells). This is followed by defining smooth functions along part of the η-curves bounding the target cells. The sign of these smooth functions depends on the direction in which the η-curve is to be displaced. We found convenient to use parabolas for this purpose. More precisely, a grid point (ξc , ηd ) belonging to the target cell and the η-curve

2514

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

(η = ηd ) containing it are chosen. A parabola with its domain along part of the grid line η = ηd is constructed. This parabola is only nonzero on the interval (ξb , ξu ) whose middle point is ξc . If the value at ξc = (ξu + ξb )/2 is Md , then the parabola is given by d (ξi ) = 2

Md Md ξb + ξu ), (ξi − ξb ) − 4 (ξi − ξb )(ξi − 2 ξu − ξ b 2 (ξu − ξb )

u ≤ i ≤ b.

(9)

If Md > 0 the parabola takes only positive values. If Md = 0, then d = 0. For Md < 0, the parabola takes only negative values. For a local modification of the spacing among the ηd -curve and some neighboring ηj -curves, parabolas j (ξ) are also constructed along a portion of the neighboring ηj -curves whose spacing will to be modified. Then, the j values are added to the original φ (which is denoted as φold ) for all the grid points in the subregion H, as indicated in (8). For smoothness of the new grid, a gradual change in the extreme value Mj of the parabolas j is also required. Once, the new control function φnew is defined, the extended BCGC algorithm is reapplied. It employs the final BCGC grid, which is obtained from using the original control functions, as initial grid. In Fig. 4, the smooth grid depicted at the right was obtained by constructing parabolas j with vertices at the point (ξc , ηj ) (j = d − 4 . . . d), respectively. Each parabola includes the grid points (ξi , ηj )(i = c − 3 . . . c + 3). They were constructed such that j > 0 over the selected η-curves to the left of the astroid. Also, j < 0 over the selected η-curves to the right of the astroid. The separation among η-curves immediately adjacent to the three-leafed rose boundary is also shortened by constructing parabolas j , over a portion of these neighboring η-curves, with Mj < 0. For the particular experiment illustrated in Fig. 4, values of Md ∈ (−0.1, 0.1) were enough to obtain the desired modification for φ. The resulting control function φnew is illustrated at the bottom of Fig. 5. A comparison of the two grids in Fig. 4 shows the effectiveness of the application of the modified control functions to the subregions presenting undesired separation of η-curves. Space control between neighboring ξ-curves can also be established by modifying the control function ψ. The mechanism is analogous to the one just described for the η-curves. 3. Grid quality The reliability of the grids used to compute numerical approximations of field variables depends on how well they satisfy certain properties such as smoothness, non-self-overlapping, orthogonality and appropriate cell areas among others. In [4], important parameters such as the maximum deviation from orthogonality (MDO) MDO = max|90◦ − θi,j |, 1 ≤ i ≤ N1 − 1, 2 ≤ j ≤ N2 − 1, and the average deviation from orthogonality (ADO)   N 1 −1N 2 −1  βi,j 1 ◦ , (|90 − θi,j |), θi,j = arccos ADO = (N2 − 1)(N1 − 2) (αi,j γi,j )1/2 i=1 j=2 were introduced. The angle θi,j is a discrete approximation for the local distortion angle between grid lines at the point (xi,j , yi,j ). Here, we introduce new parameters that measure the smoothness of the grids. This is done by considering approximations of the tangent vectors at every interior grid point xi,j = (xi,j , yi,j ) along the η and ξ coordinate curves . In particular, for a given η-curve (η=constant) two approximations of the tangent vector xξ (xi,j , yi,j ) at the grid point xi,j , using one-sided derivatives, are computed. More precisely, we approximate the tangent vector at xi,j in the ξ-direction by using the following backward and forward finite difference approximations, fwd

xξ

(xi,j ) = xi+1,j − xi,j

xbwd ξ (xi,j ) = xi,j − xi−1,j .

(10) fwd

The angle κi,j between the tangent vector approximations xbwd and xξ ξ

at the grid point xi,j satisfies

fwd

cos κi,j =

xbwd ξ (xi,j ) · xξ xbwd ξ (xi,j )

(xi,j ) fwd xξ (xi,j )

(11)

A good measure of the grid smoothness in the ξ-direction is given by the deviation of the angle κi,j from the value zero. Therefore, we define the maximum deviation from smoothness (MDS) as MDSξ = max|κi,j |,

1 ≤ i ≤ N1 − 1, 2 ≤ j ≤ N2 − 1,

(12)

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2515

Table 2 Generalized BCGC grid quality for a domain with two holes (Fig. 1). Grid size

81 × 81

101 × 101

121 × 121

161 × 161

201 × 201

Jmin Jmax Javg ADO ADSξ ADSη

6.45e−04 3.71e−01 1.92e−01 5.9◦ 5◦ 0.2◦

5.04e−04 2.88e−01 1.23e−01 5.6◦ 4◦ 0.2◦

4.10e−04 2.38e−01 8.53e−02 5.6◦ 3.3◦ 0.1◦

4.40e−05 1.30e−01 4.80e−02 5.3◦ 2.5◦ 0.1◦

2.80e−5 1.07e−01 3.07e−02 5.5◦ 2◦ 0.1◦

and the average deviation from smoothness (ADS) in the ξ-direction as ADSξ =

N 1 −1N 2 −1  1 |κi,j |. (N2 − 1)(N1 − 2) i=1 j=2

Analogously, the maximum deviation MDSη and the average deviation from smoothness ADSη in the η-direction can be defined from one-sided approximations of the tangent vector at xi,j in the η-direction. In Table 2, the mesh quality for the grids contained in the previous Table 1 is reported by using the above parameters. The Jacobian of the transformation J = xξ yη − xη yξ serves to determine non-self-overlapping [12] and also serves to approximate the area of the cells. Bounding curves describing the epicycloid and astroid experience abrupt changes of 360◦ at the cusps points. Neighboring grid lines also experience abrupt changes as they go near these singular points. As a consequence, the maximum deviation from orthogonality is relatively high for domains containing these type of holes. This is the case for the grids previously obtained and reported in Table 1. However, the BCGC grids obtained have an average deviation from orthogonality (ADO) below 6◦ . On the other hand, the average deviation from smoothness (ADS) is 5◦ or less in the ξ-direction while in the ηdirection is only 0.2◦ or less, for the grids in Table 1. It is greater in the ξ-direction because the bounding curves for the complexly shaped holes (astroid and three-leafed rose) are η-curves and they are strongly non-smooth at the corners or cusps. However, from Table 2, it is noticed that the smoothness in the ξ-direction greatly improves as the mesh is refined. 4. Scattering from multiple obstacles The BCGC grids obtained in the previous section are now used to obtain numerical approximations of the acoustic pressure field for multiple scattering problems. The initial boundary value problem in BCGC curvilinear coordinates canbe written as (psc )tt =

c2 (α(psc ))ξξ − 2β(psc )ξη + γ(psc )ηη + αψ (psc )ξ + γφ (psc )η , J2

psc (ξs , ηs , t) = −pinc (ξs , ηs , t) = −eik(x(ξs ,ηs ) cos δ+y(ξs ,ηs ) sin δ) e−iwt , psc (ξ, η, 0) = 0, (psc )t +

(psc )t (ξ, η, 0) = 0,

(ξ, η) ∈ D , t > 0,

(x(ξs , ηs ), y(ξs , ηs )) ∈ S

(ξ, η) ∈ D

c c psc ((xyη − yxη )(psc )ξ + (yxξ − xyξ )(psc )η ) + → 0, rJ 2r

(13) (14) (15)

when r → r∞ .

(16)

where D is the computational domain with coordinates ξ and η, and S represents all the closed curves bounding the multiple obstacles. The function pinc represent an incident pressure plane wave propagating along the xy-plane The variable δ corresponds to the angle of inclination of the incident wave. The parametersk, ω, and c = ω/k corresponds to the wave number, the frequency, and the wave speed respectively. The incident wave is scattered from several infinitely long cylindrical objects of arbitrary cross-section, whose length extends to infinity in the z-direction. Eq. (16) is the curvilinear version of an absorbing boundary condition of O(1/r 5/2 ) [2]. The IBVP (13)–(16) is numerically solved using the leap-frog second order, marching in time, finite-difference scheme in curvilinear coordinates. The following discrete equation, for the scattered pressure field at the interior points

2516

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

of the physical domain D, results 2 2 n (psc )n+1 i,j = 2(1 − αi,j ζi,j − γi,j ζi,j )(psc )i,j     ψi,j ψi,j φi,j 2 αi,j 1 + + ζi,j (psc )ni+1,j + αi,j 1 − (psc )ni−1,j + γi,j 1 + (psc )ni,j+1 2 2 2   φi,j βi,j

+ γi,j 1 − (psc )ni,j−1 − (psc )ni+1,j+1 − (psc )ni+1,j−1 − (psc )ni−1,j+1 + (psc )ni−1,j−1 2 2

− (psc )n−1 i,j .

(17)

The symbols αi,j , βi,j , γi,j and Ji,j correspond to the discrete approximations of the scale metric factors α, β and γ and the Jacobian J at the point (ξi , ηj ). The parameter ζi,j is such that ζi,j = c t/Ji,j . Eq. (17) is used to obtain the numerical approximations at all the grid points for 1 ≤ i ≤ N1 ,1 ≤ j ≤ N2 − 1, except for those points corresponding to obstacles’ boundaries. Special care should be exercised for the computation at the interior points located on the branch cuts and at the block interfaces (annular cuts). At these points, periodic conditions for the approximations of the derivatives should be used. Another discrete equation is employed at r = r∞ (j = N2 and 1 ≤ i ≤ N1 ). This is obtained by combining the finite difference approximations of the governing Eq. (13) and the absorbing boundary condition (16). In reference [13], a detailed derivation of this discrete equations is provided. At each time step, Eq. (17) and the discrete equation valid at r = r∞ are used until a harmonic steady state of the scattered pressure wave psc is obtained. The computation stops when the max point-wise difference between two successive time levels falls below a specified tolerance . Thus, the stop criteria that helps to determine the uniform convergence of the successive approximations of psc to a steady state N−1 N is given by Es = Max |(psc )i,j | − |(psc )i,j | < . i,j

The leap-frog finite difference method is explicit and should satisfy stability conditions. In reference [7], it is shown that a sufficient condition for stability is given by t ≤

x y c( x2

1/2

+ y2 )

.

(18)

This means that the time step size chosen should be smaller than the ratio of the area of the cells to the diagonal of the cells. Hence, the time step size is determined by those cells with the smallest ratio. As shown in [13], the BCGC grid generator algorithm plays a key role in modifying a given grid to satisfy the stability condition. In many instances, there is no need to reduce the time step t, or to increase the number of cells, to improve the accuracy of the numerical method. An adjustment of the spacing control parameters fl and/or fr is usually enough to guarantee stability. As a result, the computational cost is less than the cost incurred by employing commonly used structured grids. 5. Numerical results The finite difference method described in the previous section is now applied to multiple scattering problems. We consider a plane wave propagating in the xy-plane in the positive direction of the x-axis with frequency k = 4 and wave speed c = 1. This incident wave is scattering from a region containing three cylindrical obstacles with their crosssections bounded by the following closed curves: epicycloid, astroid, and four-petal rose. Their parametric equations are given in [13]. 5.1. Total pressure field First, we study a configuration having an epicycloid in the middle, the astroid located to the right and below at 45◦ angle, and the four-petal rose to the right and above at 45◦ angle too, as seen in Fig. 6. We call this region “nonaligned three obstacles” (NATO). Fig. 6 shows a close view of a 221 × 121 grid. This grid is generated by decomposing the domain in three single-hole blocks. They are constructed by extending the two-block procedure described in Section 2 to three blocks. According to the smoothness parameters defined in Section 3, this grid is highly smooth. In fact, the average deviation from

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2517

Fig. 6. Generalized BCGC grid for the nonaligned three obstacles. Table 3 Numerical simulations for NATO configuration. Tolerance = 5e−4, t = 2e−3. Grid size

ADO

ADSξ

ADSη

Jmin

ttheo

Avgerr

Iter

281 × 181 301 × 201 341 × 241 381 × 281

4.9◦

1.5◦

0.2◦

9.81e−06 5.97e−06 1.63e−06 9.13e−07

7.88e−03 6.87e−03 5.77e−03 5.15e−03

3.86e−04 3.85e−04 4.04e−04 3.21e−04

40,001 40,001 40,001 40,001

4.7◦ 4.3◦ 3.6◦

1.4◦ 1.2◦ 1.1◦

0.2◦ 0.2◦ 0.1◦

smoothness (ADS) in the two curvilinear coordinate directions are ADSξ = 1.8◦ and ADSη = 0.4◦ , respectively. For more refined grids the deviation form smoothness is even much smaller as seen in Table 3. The amplitude of the total pressure field for the NATO obstacle problem is illustrated in Fig. 7. The grid used for the numerical experiment has 301 × 201 grid points. A maximum amplitude of 3.32 is reached at the grid point (5.9,7),

Fig. 7. Amplitude of total pressure field for the NATO configuration.

2518

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

Fig. 8. Amplitude of total pressure field for the ATO configuration.

very close to the four-petal rose boundary. There is not a clear shadow zone. Indeed, the incident wave is deviated by the epicycloid surrounding obstacles, penetrating to the otherwise unperturbed shadow zone, creating a complex wave pattern. However, the amplitude there is still much smaller than on the rest of the domain. In Table 3, numerical results for the NATO obstacles on generalized BCGC grids of different sizes are recorded. It is important to note that ADS is below 1.5◦ in the ξ-direction and 0.25◦ in the η-direction. Due to the cusp points present in the epicycloid and astroid curves, the maximum deviation from orthogonality is very high about these points. However, the average deviation from orthogonality (ADO) is below 5◦ . The sixth column contains the theoretical upper bound of the time step t for stability. It decreases as the grid is refined. The theoretical minimum found for a 381 × 281 grid is t = 5.15 × 10−3 . We employed a common t = 2 × 10−3 for the numerical simulation on all grids. The seventh column contains values for the Avgerr , which is the average error of all the point-wise differences between two successive time iterations. The tolerance used by the stop criteria is 5 × 10−4 , but it is not activated until the wave has traveled back and forth twice over the entire domain to eliminate all transient effects. This occurs at about 40,000 iterations. A second numerical simulation was also performed. Its domain consists of the same set of obstacles analyzed for the NATO problem. However, they are aligned along the vertical y-axis, as seen in Fig. 8. We call this configuration “aligned three obstacles” (ATO). The other grid and field parameters are identical to the previous NATO configuration. But, the incoming wave incident angle is δ = 45◦ for this case. The maximum amplitude of 3.42 for the total pressure still occurs about the four-petal rose boundary, at the node (−0.94,9.2). An interesting pattern of waves is observed between the obstacles. The smoothness of the grid allows a smooth wave propagation throughout the entire domain. 5.2. Radar cross-section The radar cross-section can be easily approximated from the discrete approximation of the harmonic steady state t scattered pressure in the far field (psc )N i,N2 (i = 1, . . . N1 ). Following a procedure detailed in [13], as a first step, a set of coefficients cˆ m areobtained N1 −1 1  −imθ t cˆ m = (psc )N , i,N2 e N1 l=0

for m =

−N1 , . . . , (N1/2 ) − 1. 2

(19)

These finite series can be directly evaluated, or a FFT algorithm can be used to compute them. Then, the radar cross-section is approximated by the semi-analytical formula

Â0 (θ) =

2 −iπ/4 e πk

N1 /2−1 m=−N1 /2

bˆ m (−i)m eimθ ,

(20)

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

2519

Fig. 9. Radar cross-sections for the NATO (left) and ATO (right) configurations.

where the coefficients bˆ m are obtained from the coefficients cˆ m and Hankel function of order m, Hm(1) evaluated in the far field. Indeed, bˆ m = cˆ m /Hm(1) (kr∞ ). The formula (20) is extremely accurate as shown in [3]. The error in the approximation of A0 (θ) using (20) depends almost entirely upon the error made in the approximation of the coefficients bˆ m . As expected, the radar cross-sections for these two numerical experiments exhibit a complex pattern, as shown in Fig. 9. For the NATO configuration the maximum value of A0 (θ) is 7.73 at θ = 0◦ , which is the direction of propagation of the incident wave. Similarly, A0 (θ) = 7.33 at θ = 44.6◦ for the ATO configuration with an incident wave propagating at a 45◦ angle. 6. Concluding remarks A practical grid generation algorithm for two-dimensional multiply connected regions with multiple holes has been constructed. It can handle regions with bounding curves that includes corners and cusps that constitutes a challenge for any grid generator. There are two novel key features of the generalized algorithm. First, the smoothing process between adjacent blocks to obtain globally smooth grids. Secondly, the posteriori local adjustment of grid cell sizes to improve the grid quality. The generalized BCGC grids have been successfully used to compute the scattered pressure field over complex geometry configurations. In most experiments, stability conditions were obtained by adjusting grid control parameters without decreasing the time step. The generalized BCGC grids will have an impact not only in wave scattering, but in the numerical solutions of physical problems defined on multiply connected regions with several holes. In particular, computational fluid dynamics can benefit from these BCGC grids to numerically simulate fluid problems that may be impossible to simulate using current grid generation techniques. References [1] V. Akcelik, B. Jaramaz, O. Ghattas, Nearly orthogonal two-dimensional grid generation with aspect ratio control, J. Comp. Phys. 171 (2001) 805–821. [2] A. Bayliss, E. Turkel, Boundary conditions for the numerical solution of elliptic equations in exterior regions, SIAM J. Appl. Math. 42 (1982) 430–451. [3] O.P. Bruno, E.M. Hyde, Higher-order Fourier approximation in scattering by two-dimensional, inhomogeneous media, SIAM J. Numer. Anal. 42 (2005) 2298–2319. [4] L. Eca, 2D orthogonal grid generation with boundary point distribution control, J. Comp. Phys. 125 (1996) 440–453.

2520

V. Villamizar, S. Acosta / Mathematics and Computers in Simulation 79 (2009) 2506–2520

[5] A. Khamayseh, A. Kuprat, Surface grid generation systems, in: J.F. Thompson, B.K. Soni, N.P. Weatherill (Eds.), Handbook of Grid Generation, CRC Press, 1999, pp. 9.1–929.929. [6] E. Manoha, R. Guenanff, S. Redonnet, M. Terracol, Acoustic scattering from complex geometries, in: Proceedings 10th AIAA/CEAS Aeroacosutics Conference, AIAA, 2004, pp. 1543–1551. [7] G. Sewell, The Numerical Solution of Ordinary and Partial Differential Equations, Wiley, 2005. [8] S.E. Sherer, J.N. Scott, High-order compact finite-difference methods on general overset methods, J. Comput. Phys. 210 (2005) 459–496. [9] S. Sherer, M. Visbal, High-order overset-grid simulations of acoustic scattering from multiple cylinders, in: Proceedings of the Fourth Computational Aeroacoustics (CAA) Workshop on Benchmark Problems, NASA/CP-2004–212954, 2004. [10] J.F. Thompson, A general three-dimensional elliptic grid generation system on a composite block structure, Comput. Meth. Appl. Mech. Eng. 64 (1987) 377–411. [11] J.F. Thompson, B.K. Soni, N.G. Weatherill, Handbook of Grid Generation, CRC Press, 1999. [12] V. Villamizar, O. Rojas, J. Mabey, Generation of curvilinear coordinates on multiply connected regions with boundary singularities, J. Comput. Phys. 223 (2007) 571–588. [13] V. Villamizar, M. Weber, Boundary-conforming coordinates with grid line control for acoustic scattering from complexly shaped obstacles, Numer. Meth. Part. Differ. Equ. 23 (2007) 1445–1467.