A reduced domain strategy for local mesh movement application in unstructured grids

Applied Numerical Mathematics 61 (2011) 1001–1016

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A reduced domain strategy for local mesh movement application in unstructured grids Masoud Darbandi ∗ , Nematollah Fouladi Center of Excellence in Aerospace Systems, Department of Aerospace Engineering, Sharif University of Technology, Azadi Ave., P.O. Box 11365-8639, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 5 June 2010 Received in revised form 26 March 2011 Accepted 15 April 2011 Available online 11 May 2011 Keywords: Unstructured grid Local grid movement or deformation Moving boundary Automatic grid partitioning Grid data structure

a b s t r a c t Automatic control of mesh movement is mandatory in many fluid flow and fluid–solid interaction problems. This paper presents a new strategy, called reduced domain strategy (RDS), which enhances the efficiency of node connectivity-based mesh movement methods and moves the unstructured grid locally and effectively. The strategy dramatically reduces the grid computations by dividing the unstructured grid into two active and inactive zones. After any local boundary movement, the grid movement is performed only within the active zone. To enhance the efficiency of our strategy, we also develop an automatic mesh partitioning scheme. This scheme benefits from a new quasi-structured mesh data ordering, which determines the boundary of active zone in the original unstructured grid very easily. Indeed, the new partitioning scheme eliminates the need for sequential reordering of the original unstructured grid data in different mesh movement applications. We choose the spring analogy method and apply our new strategy to perform local mesh movements in two boundary movement problems including a multi-element airfoil with moving slat or deforming main body section. We show that the RDS is robust and cost effective. It can be readily employed in different node connectivity-based mesh movement methods. Indeed, the RDS provides a flexible local grid deformation tool for moving grid applications. © 2011 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction There are many fluid flow simulations, such as the aerodynamics shape optimization [1,13], fluid–structure interaction [24,43], store separation [15,19], moving objects [31,44], and free surface flows [12,28], in which the grid movement is mandatory. In such fluid flow simulations, the computational domain should be suitably moved to account for the deformation occurred at its boundaries. Generally, the motion or deformation of some parts of the computational domain, mostly the inner boundary, is known in grid movement problems and that the rest of grid should be moved or deformed suitably to accommodate the imposed displacements. However, the point is that the grid movement should be accomplished in a manner to guarantee the grid quality while minimizing the grid computations. These issues have been the concerns of many past researchers in mesh movement and grid deformation studies, e.g. [9,18,21,40,49]. Literature shows that there are two major unstructured mesh movement approaches used in fluid flow simulations widely. The first approach benefits from the connectivity of nodes in the computational domain to update the grid. The PDE-based mesh update methods and the spring analogy-based methods are categorized in this approach. These methods solve a system of equations as large as the number of grid nodes in the computational domain to update the grid. The PDEbased methods usually use some kinds of Laplacian (the second-order operator) smoothing procedure to update the grid [4, 27]. Furthermore, the use of biharmonic operator (the fourth-order generalization of the Laplacian operator) has improved

*

Corresponding author. Tel.: +9821 6616 4644; fax: +9821 6602 2731. E-mail address: [email protected] (M. Darbandi).

0168-9274/$30.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2011.04.005

1002

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

the mesh deformation quality through controlling the normal mesh spacing along the boundary [20,42]. In the spring analogy-based methods, the basic concept is to replace the connectivity of nodes with a network of springs or elastic solids inserted between each two neighboring grid nodes [3]. The basic method was suffering from the lack of collapse control mechanisms. Farhat et al. [11,14] considered the torsional spring context to avoid the grid crossover problem. Murayama et al. [30] further extended the method to generalize it for special wall treatments. Bottasso et al. [8] developed a 3D ballvertex method by introducing additional linear springs to resist the motion of a mesh vertex towards the corresponding region-opposed face. In case of large deformations, some researchers considered the grid as an elastic solid [2,23,32,48]. This consideration improved the robustness of mesh movement algorithm with more guaranteed mesh validity, however, in the expense of a much higher computational cost. Contrary to the node connectivity-based methods, there are mesh movement methods, which do not use the connectivity of nodes. For example, in the point-by-point method, the displacement of each point is determined based on its relative position with respect to the domain boundary [34]. However, this method has been applied to the boundary nodes of multiblock mesh movement cases. In the radial basis functions (RBFs)-based mesh movement method, the interpolation function of boundary displacements is constructed to extrapolate those displacements into the entire grid domain [10,22,35]. The RBF-based methods are established based on decaying global interpolation, which would normally be more robust than the basic spring analogy method. However, they are expensive for large mesh applications due to large data manipulation requirements. It is because the required matrix for mesh motion is strictly N surface points × N volume points. Rendall and Allen [36–38] improved the basic method through a point selection process to incorporate the reduced number of surface points for defining the surface motion. They indicated that the smart data reduction schemes would choose a significantly reduced data set based on minimizing an error function. Consequently, this would reduce the cost of mesh motion effectively in terms of both time and memory requirements. The method is currently underway to be improved suitably and it may be regarded as an important mesh movement strategy in near future. Comparing with RBF-based methods, the connectivity-based methods have been applied widely in past researches. Therefore, we focus on the node connectivity-based approaches in this paper. During moving or deforming a body enclosed with an unstructured grid, the grid deformation is usually started from the nodes located at the moving or deforming section. Therefore, the equations of motion should be applied iteratively to force the boundary movements into the unstructured grid distributed in the domain. However, due to the unstructured nature of the connectivity matrix, the equations of motion must be applied to all nodes located in the domain during each incremental boundary movement. We call this procedure the complete domain strategy (CDS). Definitively, the CDS is quite expensive for large mesh applications. Therefore, there have been numerous efforts to reduce the inefficiencies of the CDS movement solutions, e.g. Refs. [25,26,41,46,47,50]. Our literature review shows that many grid movement users choose the CDS to perform their mesh movements. Therefore, it is highly valuable to found new local mesh movement strategies to decrease the CDS inefficiencies. In this paper, we present a local mesh movement strategy, which moves the grid only in a small region close to the moving boundary. We call our strategy the reduced domain strategy (RDS). The RDS also benefits from our past innovated schemes, which are described shortly. Evidently, it is very advantageous to have easy access to the grid data near the deformed boundary in any local unstructured grid movement. This access helps to move the grid very locally without manipulating all mesh data. As is known, the grid data in a local region cannot be readily separated from the rest of unstructured grid domain. Normally, this task needs implementing suitable search processes in the entire unstructured grid data. Of course, the mesh ordering can provide the possibility of direct access (without searching the entire mesh data) to the target mesh regions [7]. Indeed, the mesh ordering process normally needs considerable additional computer memory storage. Meanwhile, it does not have flexibility enough to work effectively in alternative regional deformation tasks. In most of these cases, the computational inefficiencies increase due to a heavy sequential reordering requirements. To overcome these inefficiencies, Fouladi et al. [16] developed a new ordering-based flexible moving grid strategy. This strategy divides the deformation task into two steps. At the first step, it modifies the unstructured mesh data structure through renumbering its substructures in an ordering manner without requiring additional memory storage. At the second step, a suitable flexible deformation algorithm is applied to manage the grid movement into the computational domain. This strategy alleviates the need for heavy sequential grid reordering tasks. This procedure needs considerably less computational time if it is compared with past ordering strategies. This method was first developed in the context of an edge-based unstructured connectivity matrix [16]. However, it was suffering from the complexity in both code development and its implementation. Also, the total time for mesh ordering was relatively considerable and expensive. To improve the preliminary strategy, we later developed a more advanced ordering and renumbering strategy [17]. The improved strategy was basically founded in the context of a cell-based unstructured connectivity matrix. We called it a quasi-structured data structure (QSDS) strategy. In the QSDS, the un-structured data structure (USDS) of an original unstructured grid is broken into several bands of data around the main object. The QSDS is briefly described in Section 3.1. Following Refs. [16,17], we further extend a cost-effective local grid deformation strategy in this work. Using the QSDS, we first develop a new unstructured mesh partitioning strategy. This avoids unnecessary grid computations via dividing the computational domain into two active and inactive zones. The former zone includes the moving boundary and its adjacent moving grid, which is affected by the primitive boundary movement. Second, we develop a cost-effective algorithm to expand the active zone during the mesh movement or deformation procedure. In this algorithm, we start from the cells adjacent to the deformed boundary and expand the active zone gradually into the main domain until achieving negligible node

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1003

displacement magnitudes at the front of active zone boundary. As was mentioned before, we refer to our local moving grid strategy as the reduced domain strategy (RDS). Once more, the RDS is prone to be applied to any node connectivity-based mesh movement methods, e.g., the basic spring analogy, the torsional spring analogy, and the linear elasticity methods. In this paper, we simply choose the spring analogy method as a basic mesh deformation method to introduce the RDS and to elaborate its advantages. It should be re-emphasized that the choice of spring analogy does not mean that the RDS is restricted to this mesh movement method at all. Indeed, the RDS can be readily applied to different engineering problems engaged with unstructured grid movement, e.g., the aerodynamic shape optimization [5,29,33,39]. However, to reduce the size of paper, the RDS is examined for treating only two common mesh movement or mesh deformation cases. It is a multielement airfoil with moving slat and/or deforming main body section. Comparing with the CDS, the RDS provides rapid mesh movements with less memory requirements. 2. Mesh movement using the complete domain strategy (CDS) To illustrate the advantages of RDS, we need choosing a basic mesh movement method. As was mentioned in the introduction section, we simply choose the basic spring analogy method to carry out our study in both CDS and RDS. Therefore, we briefly describe the spring analogy method here. As is known, the spring analogy methods can be classified into two main vertex and segment types [6]. The vertex springs are always under tension unless the spring length is zero. They are mostly used for mesh smoothing purposes. However, the segment springs have zero tension before any deformation. The displacements occurred at the boundary points would create spring forces in the region around those boundary points. According to Hook’s Law, the force of neighboring nodes F i to node i can be calculated from

F i =

m 

ki j (δ j − δi )

(1)

j =1

where j counts the neighboring nodes of the ith node and δ j and δi present the displacement vectors for nodes j and i, respectively. As a result of one movement, the neighboring internal nodes should be displaced until reaching an equilibrium status in the entire domain. Given a set of nodal displacements, it eventually results in a set of linear algebraic equations as follows: m 

ki j δin+1 =

j =1

m 

ki j δin

(2)

j =1

where ki j represents the spring stiffness coefficient between node i and j. It is defined as



ki j = 1/ ( X i − X j )2 + (Y i − Y j )2

1/ p

(3)

where p = 2 [45]. It means that the spring stiffness coefficient is usually taken proportional to the inverse of segment length. Eq. (2) needs to be solved iteratively until reaching the equilibrium condition at all nodes. Then, the new position of each node is updated using

 new = X old + δn X i i i

(4)

The segment model, which is utilized in this paper, is suitably applied to move the interior nodes iteratively. The movement equations (like those of spring-based methods or other analogues smoothing or movement methods) govern the motion after each infinitesimal boundary movement. Then, the changes forced on the boundary nodes are gradually influenced into the solution domain. Therefore, the equations of motion need to be solved iteratively in the entire solution domain to propagate the movements enforced at the boundary into the entire domain and consequently to the nodes far from the primitive displaced boundary. To launch our study, we choose the multi-element airfoil as our test case. Fig. 1 shows two different views of the initial unstructured grid around the chosen multi-element airfoil. This grid consists of 8921 nodes and 17 361 triangles. The slat moves from position 1 to 2, see Fig. 2. The displacement length is d. The original distance between the slat and main airfoil was D. Although there is no limit in the magnitude of d, we typically chose d = (0.5, 1.5, and 2.5) D to illustrate small, moderate, and large deformations in grid. In fact, we only want to illustrate the basis of our local mesh movement strategy in the present work. To move the grid consistently with the slat motion, we need solving Eq. (2) for the entire grid nodes repeatedly. Fig. 3 shows the propagation of grid displacement front in the domain at the end of four different movement iteration cycles of 100, 150, 200, and 300. Figs. 3(a) and 3(b) respectively present the propagations of 5% and 10% of the slat displacement magnitudes in the computational domain. For example, the solid line in Fig. 3(a) splits the grid into two interior and exterior zones. The grid displacement is larger than 5% of the main slat displacement in the interior zone at the end of 100th movement iteration cycle. Similarly, the grid displacement is less than 5% of the main slat displacement in the exterior zone at the end of 100th movement iteration cycle. The figure also shows that the displacements originate from the displaced slat and propagate gradually into the domain during iterative solution of Eq. (2). A careful inspection in Fig. 3 shows that a small region in the vicinity of slat can mainly be affected by the imposed local

1004

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 1. Initial unstructured grid distributed around a multi-element airfoil: (a) the entire domain; (b) a closer view.

Fig. 2. The displacement of a multi-element airfoil slat.

Fig. 3. Grid displacement propagation in the domain at four different mesh movement cycles considering two displacement magnitudes of (a) 5 percentages and (b) 10 percentages of the main slat movement.

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1005

Fig. 4. A schematic view of element layers and node lines in the QSDS.

boundary displacement. In the CDS, the equations of motion need to be applied at all nodes in the domain in each iteration cycle. In other words, the movement equations are solved for both the affected and the non-affected grid nodes in the entire domain. This consequently results in computational inefficiency and a very time consuming mesh movement procedure. The issue seriously affects the solution of real moving boundary problems, where a large number of local grid movements is inevitable. In the next sections, we develop a new local ordering-based mesh movement strategy, which reduces these inefficiencies. 3. Introducing a cost-effective mesh partitioning framework To reduce the CDS inefficiencies, we try to eliminate unnecessary computations for stationary nodes in each iteration cycle by dividing the grid domain into two active and inactive zones. In this regard, we first convert the original unstructured mesh data structure (USDS) into a quasi-structured data structure (QSDS). Then, we construct the active zone in the QSDS. The details are provided in the following subsections. 3.1. The quasi-structured data structure (QSDS) scheme As is described in Ref. [17], the QSDS reorders the elements and nodes of an original USDS suitably into arrays of element layers and node lines bounded the main object in the domain. In this regard, the original unstructured connectivity matrix is broken into several parts. Each part represents the data for a specific layer wrapped around the main object. Fig. 4 shows two element layers and four node lines around the object. As is seen, each element layer is bounded by two interior and exterior node lines. After constructing the element layers, we renumber all the element and node indices in each layer and its two bounding lines in one chosen direction, clockwise or counterclockwise. This eventually yields a semi-structured data pattern for the constructed element layers and node lines. The elements in every element layer and the nodes at every node line are renumbered from a minimum index value to a maximum one. Also, the vertices of every element in each element layer are renumbered in a manner, in which the first and third columns of the new generated connectivity matrix in QSDS normally address the nodes at its interior and exterior node lines, respectively. This accelerates the search for grid substructures (i.e., elements, edges, and nodes) effectively. We benefit from the QSDS features and establish a cost-effective mesh partitioning procedure. We extend a cell-based data structure for an unstructured grid in the QSDS context without requiring additional memory to store the data. Similar to the original data structure, the new extended data structure only needs storing the node coordinates data, P I = ( X I , Y I ), and an element connectivity matrix knowledge, E J = { N J K , K =1,2,3 }. Indices I and J denote node and element numbers, respectively, and K counts the element vertices. To be consistent with the unstructured grid data, the dimensions of data structure matrices are the same for the original and improved frames in this method. However, two auxiliary vectors of L P and L E are defined and used to address each node index and each element index, respectively, appeared at the end of each line and each layer. Of course, they provide valuable information, which permits a quick access to the grid substructures in the extended grid data structure. In other words, the method simplifies any address to the grid substructures because it restricts the grid search to a local confined region of neighboring layers and lines, instead of the entire mesh data [17]. As is indicated in this reference, the data ordering is achieved without needing any additional memory storage capacity. 3.2. Constructing an active zone using the QSDS To avoid the CDS computational inefficiencies, we develop a new framework in which the movement equations are solved only for the nodes affected by the boundary displacements, i.e., nodes located in the active zone.

1006

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 5. A schematic view of an element sub-layer (ESL) and its two node sub-lines (NSLs).

3.2.1. The active zone framework To find the active zone, we need partitioning the grid automatically, of course, with the least computational cost. Normally, the initial active zone is located at a portion of the first element layer attached to the displaced boundary. This primitive active zone is gradually expanded into its neighboring element layers during the mesh movement computations. We use element sub-layer to refer to each new layer constructed around the active zone. Of course, the active zone is expanded due to the occurred boundary displacements. Fig. 5 shows a schematic view of an element sub-layer (ESL) constructed during one iteration cycle. In this case, the displaced boundary has occurred at the main object. As is seen, the ESL is located between two interior and exterior node sub-lines (NSLs). The ESLs and NSLs are extracted from the QSDS, which was constructed using the original USDS. However, contrary to the element layers and node lines in QSDS, which mostly encompass the main object, see Fig. 4, the ESLs and NSLs do not necessarily encompass the main object. Their starting and ending elements or nodes are mostly coincided at the main object, if the displacements occur at the main object surface. By default, the first ESL connected to the displaced boundary would be the initial active zone. In other words, the initial active zone will be a part of the first element layer wrapped around the main object. Therefore, the solution of grid movement equations is limited to this ESL at the first iteration cycle. In the next iteration cycles, the next ESLs are created step-by-step. Following this step-by-step influence, the active zone expands gradually into the main domain and this allows a suitable propagation of the initial displacement into the unstructured grid distributed in the domain. Obviously, the first NSL located at the main object surface is subject to the boundary movement or deformation, see Fig. 6. The figure shows a schematic view of the constructed ESLs and NSLs. There are four NSLs and three ESLs in this figure. As is seen, each ESL is terminated to two bounding elements in each element layer of the QSDS. They are labeled by two numbers of s and l, i.e. ±sl, in which the minus and plus signs refers to the left and right bounding elements, respectively. Also, the parameters s and l refer to the sub-layer and QSDS layer numbers, respectively. Such labels can be similarly attributed to the bounding nodes of each NSL of the original QSDS node lines. To distinguish them, we label them as ±sl N and ±sl E for NSLs and ESLs, respectively. 3.2.2. The active zone expansion in the QSDS Up to here, the ESLs and NSLs are defined in the QSDS using their bounding elements and bounding nodes, respectively. To expand the active zone, it needs to find the next bounding elements, bounding nodes and updating the bounding front data list. Having an NSL, its neighbor ESL is readily constructed because every element in the ESL has at least one vertex at that NSL. To create a new ESL, the existing active zone can be expanded in two main directions: 1) the lateral direction along the affected QSDS layers, and 2) the forward direction across neighboring unaffected QSDS layer, see Fig. 6. The flowchart in Fig. 7(a) shows the steps taken to expand the active zone in the left and right lateral directions. This flowchart indicates that the marching starts from layer l = 1 and ends up to l = s − 1. Normally, the marching is started from the bounding elements at the current front data list and is continued along every QSDS layer to search for the new bounding elements. The marching in every layer is continued until reaching the last elements, which have at least one vertex at the current active zone boundary. These elements are selected as the new bounding elements for the lth layer. Finally, the bounding front data list is updated. To expand the active zone in the forward direction, we must add a portion of the next neighboring element layer, i.e., the lth QSDS element layer, to the current active zone. Consequently, it needs finding the bounding elements at layer l = s in both left and right ends of the active zone front. To do this, we apply a consecutive estimation procedure. Fig. 8(a) shows the distribution of triangular elements located in a short length of an element layer. It indicates that there would normally be a simple relationship between the number of nodes and elements in this portion of the chosen element layer. In other words, if there are 2 × n nodes in the shaded area, there will be 2 × (n − 1) elements there. In our searching, we start from an arbitrary element of this layer and search for the bounding elements. However, the number of nodes from the starting

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1007

Fig. 6. A schematic view of four node sub-lines and three element sub-layers in the active zone and the expansion of active zone in lateral and forward directions.

Fig. 7. Flowcharts showing the expansion of active zone in the lateral (a), and forward (b) directions.

node (I 1 ) to a specified node ( I m ) in the interior node line can be calculated from I m − I 1 , where I m (= ±sl N ) indicates the current bounding nodes at the interior node line, see Fig. 8(b). Therefore, we can estimate the target bounding elements by supposing a regular distribution of triangles within a search length in this layer. As is seen in Fig. 8(b), the search length refers to a portion of chosen QSDS element layer located between the starting point I 1 and the target node I m . However, due to a quasi-structured distribution of triangles located in the layers, we must modify this estimation several times to reach element J . This minimizes the search length. We can formulate this using

C = min| I m − N J ,1 |

(5)

1008

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 8. (a) A regular arrangement of elements within a part of an element layer and (b) a schematic presentation of search length in element layer l and the target element there.

where C is an integer number. Its value is usually less than 3. We evaluate the condition only for the first vertex of the estimated element, i.e. N J ,1 . As mentioned before, the first vertex of every element in a given layer of QSDS is located at the interior node line. In this work, we use more useful relation to provide consecutive estimations of the target bounding elements as follows:



J sr =

J 1 + 2 × ( I m − I 1 ),

r=1

J sr −1 − 2 × ( N J sr −1 1 − I m ), r > 1

(6)

where r denotes the iteration number, J sr is the rth estimation of one bounding element, and J 1 indicates the index of starting element at layer l. The flowchart in Fig. 7(b) shows the expansion of active zone in the forward direction. Starting from an arbitrary element at the lth element layer, we can consecutively estimate the new bounding elements using Eq. (6) and compute the search length using Eq. (5) until reaching a minimum search length. This is achieved when the value of C does not change during two consecutive estimations, i.e. C = 0. During the consecutive estimations, the search length in the lth layer gradually reduces and the target bounding elements are eventually detected. In practice, our experience shows that the searching process reaches near the bounding elements very quickly after few successive estimations provided by Eq. (6). This is specifically shown in Section 5. In the next step, we need extending a local neighboring condition, which searches just in a small part of the proposed layer (i.e., the lth layer) to find the bounding elements. Finally, these bounding elements are added to the front data list. It is interesting to note that we do not need storing all the elements data located in the active zone at all. We only need storing the bounding elements indices of the last ESL. In any event, if one needs having access to the data of all created ESLs in specific circumstances, one can also store the bounding elements of all ESLs. Evidently, the required data storage is almost negligible if it is compared with the entire grid data storage requirement. 4. Mesh movement using the reduced domain strategy (RDS) After constructing the QSDS from the original USDS, we use our automatic partitioning scheme to determine the active zone expansion during a local mesh movement procedure. To implement this strategy more robustly, we present the extended local mesh movement algorithm in this section. According to the flowchart presented in Fig. 9, a CFD solver, e.g., a fluid–structure interaction (FSI) solver, may need performing a local mesh motion due to a local boundary movement. If so, a boundary displacement vector (BDV) should be constructed to define the locations of those solid boundary nodes subject to that boundary movement. The rest of our local mesh movement or deformation algorithm is continued as follows: 1) 2) 3) 4) 5) 6) 7)

a new ESL is constructed to expand the active zone, the movement equations (MEqs) are applied to the active zone, if the stop criterion is not satisfied, go to step 4 else go to step 6, compute absolute nodal displacement (ANDIS) at the front of active zone, if the propagation criterion is satisfied, go to step 1 else go to step 2, the global coordinates (GCs) need to be updated, and the CFD solver is re-called.

In step 1, the new boundary of the active zone is the same as the exterior NSL of the new constructed ESL over the preceding active zone. In step 2, the equilibrium equations, Eq. (2), are applied to nodes located in the new expanded

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1009

Fig. 9. Local mesh movement or deformation algorithm.

active zone. In other words, the movement or deformation equations are solved only for the active zone in each iteration cycle. This procedure may need to be repeated several times. The stop criterion for local mesh movement can be defined as exceeding a predefined maximum iteration cycle (MIC) magnitude in solving movement equations. The propagation criterion could be defined as a function of the relative displacement of nodes located at the boundary of the current active zone. This is indicated as the maximum nodal displacement (MNDIS). If the average summation of absolute nodal displacements (ANDIS) exceeds a given value at the boundary of the current active zone, a new ESL must be added to the preceding active zone. MNDIS can be defined based on the percentages of maximum nodal displacement occurred at the moving boundary. Back to Fig. 2, one may choose MNDIS = 5(%d) or 10(%d). We will study the effect of MNDIS value on the radius of active zone expansion in the solution domain shortly. 5. Results To show the performance of our local mesh movement strategy, we investigate two mesh movement problems. The first problem is the displacement of a multi-element airfoil slat, see Figs. 1 and 2. Considering a local slat motion, the mesh movement should be performed locally and confined to a small grid region around the slat. The grid includes 8921 nodes and 17 361 elements. According to the flowchart given in Fig. 9, we define a boundary displacement vector BDV to start our local mesh movement procedure. We choose MIC = 500 and MNDIS = 10. As was described before, the MNDIS = 10 means that the maximum nodal displacement at active zone boundary must be limited to 10 percentages of the main displacement performed by the slat. Fig. 10 shows the final resulting active zone. Comparing with the whole solution domain displayed in Fig. 1, the active zone is limited to a small region in front of the main airfoil. Evidently, choosing larger MNDIS values would result in smaller active zone sizes. Fig. 11 illustrates the resulting active zones choosing MNDIS = 10, 15, and 20. Fig. 12 shows more details of the node displacement distribution in the final active zone for MNDIS = 10. The displacements are non-dimensionalized with respect to the 0.01d. As is seen, the nodes in the active zone perform displacements from 10 at the front of active zone to 100 at the slat face. Fig. 13 presents the sub-layer creation number versus the movement iteration cycle number. This figure indicates that the final active zone includes 28 ESLs. Our calculation shows that the final active zone consists of 4394 elements and 2319 nodes. However, we only need 2 × 28 = 56 integer data to address the active zone in this problem. If one wishes to store all the active zone data, one will need storing 4394 + 2319 = 6713 integer data. Therefore, we only need 56/6713 ≈ 0.83% of the entire active zone data to record the active zone data. Additionally, this storage requirement is only 56/17 361 ≈ 0.32% of the entire storage required to record all grid element data. On the other hand, since the mesh displacement computation is only fulfilled for the active zone, the displacement vector, i.e. δ in Eq. (1), needs to be stored only for nodes inside the active zone instead of all grid nodes. Therefore, this leads to 1 − (2319)/(8921) = 74% less memory storage requirement. Comparing with the CDS, we can conclude that the use of RDS would result in 74% less memory requirement. Evidently, from computational time perspective, the RDS saving is overwhelming if compared with that of CDS.

1010

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 10. The final achieved active zone using MNDIS = 10, test 1.

Fig. 11. The final achieved active zones using MNDIS = 10, 15, and 20, test 1.

Fig. 12. Contours of nodal displacement inside the final active zone, MNDIS = 10, test 1.

Fig. 13. The history of element sub-layer creation (or active zone expansion) versus the movement iteration cycle, test 1.

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1011

Fig. 14. The resulted search length for bounding elements (l = s), test 1: (a) at the first estimation; (b) at the second estimation, test 1.

Fig. 15. Computational time ratio history versus movement iteration cycle number, test 1.

Fig. 16. The occurrence of cell vertex angles in the active zone of RDS and a similar region in CDS, test 1.

As was described in Section 3, to expand the active zone, the new ESLs and NSLs must be constructed iteratively and added to the active zone during the mesh movement procedure. To construct a new ESL, the movement algorithm should search for the new bounding elements. To illustrate the robustness of the current search in finding the new bounding elements, we monitor the estimations calculated by the search algorithm. Fig. 14 presents the results of using Eq. (6) to estimate the unknown left and right bounding elements (when l = s) for the created ESLs during the active zone expansion. In these figures, J m denotes the exact index of the bounding element. Therefore, the quantities J s1 − J m and J s2 − J m denote the search lengths performed in the first and second estimations, respectively. As the search length becomes narrower, the search process approaches the bounding elements. This figure shows that the search lengths converge to their minimum values very rapidly and mostly after the second estimation for every bounding elements’ set. To provide a better comparison on RDS time saving, we define computational time ratio as CTR = t RDS /t CDS , where t RDS and t CDS denote the time spent by RDS and CDS, respectively. Fig. 15 illustrates the computational time ratio for each movement iteration cycle. As was mentioned before, the grid movement equations must be solved for the entire grid domain in the CDS. This figure indicates that the RDS performs considerably better time performance than the CDS. Our calculations showed that the RDS is at least 5 times faster than the CDS, in the current mesh movement procedure considering MNDIS = 10 and MIC = 500. It should be mentioned that the active zone is gradually constructed during the boundary movement and mesh motion. In other words, we have considered the cost of local movement plus the overhead cost of computing the active zone in this figure. Fig. 16 shows the occurrence of element vertex angles appeared in the mesh moved by RDS and compares it with that of the CDS. As is shown, both RDS and CDS perform similar quality. We can conclude that the RDS eliminates the unnecessary computations without losing the grid quality achievable by employing the CDS. In other words, the RDS reduces the time for grid movement considerably while the quality of the deformed grid is the same as that of using the CDS. Fig. 17 illustrates the performance of RDS choosing various MNDIS values. Fig. 17(a) indicates that if we choose a smaller value for MNDIS the final achieved ESL number becomes larger. It means that the active zone expands more during the movement procedure if we select a smaller MNDIS value. This was also shown in Fig. 12. Fig. 17(b) shows the computational time ratio versus the MNDIS value. According to this figure, the CTR decreases as the MNDIS increases. This is quite reasonable. The next step is to investigate the performance of RDS choosing different maximum iteration cycle (MIC) values. We try MIC from 50 to 500. We also choose a small value for MNDIS, i.e., 5, to keep the same mesh quality in the resulting RDS and CDS grids for the entire range of MIC. Fig. 18 shows that CTR increases as MIC increases. It is reasonable to have farther displacement propagations inside the domain for larger MIC values. Evidently, this results in larger active zone and eventually higher time durations. Back to Fig. 3, we can conclude that the MIC must be selected sufficiently large to avoid

1012

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 17. Effect of choosing different MNDIS values on (a) active zone size; (b) computational time ratio, test 1.

Fig. 18. Effect of maximum iteration number cycle on the computational time ratio, MNDIS = 5, test 1.

Fig. 19. The performance of RDS at different slat movement positions, test 1.

Fig. 20. The initial and deformed boundaries, test 2.

large grid displacement gradients right at the vicinity of the deformed boundary. A large MIC propagates the displacements suitably far from the original displaced boundary. We further investigate the performance of RDS while changing the slat displacement from d = 0.5D to 2.5D. To achieve a mesh quality comparable with that of using CDS, we choose MNDIS = 5 and MIC = 500. The resulting CTR is presented in Fig. 19. As is seen, the CTR is almost the same despite using a relatively wide range of slat movements. Also, the computational time spent by the CDS is at least five times larger than that of the RDS in all cases. As the second test case, we study mesh movement problem due to a boundary deformation at the main airfoil section of the preceding multi-element airfoil shape. Fig. 20 shows both the initial boundary and the deformed one occurring at the upper surface of the airfoil. Zhao et al. [51] used traveling wave of flexible wall to reduce its turbulent drag. Therefore, they had to move their grid appropriately and consistent with the flexible wall motion. However, we are not interested in the reasons behind this behavior at the present. As was invoked previously, the grid deformation mission can be achieved by RDS without any need to reorder the primitive unstructured grid data. Choosing MNDIS = 10 and MIC = 100, the final achieved active zone is presented in Fig. 21. The active zone is limited to a small region at the top of the deformed boundary. Fig. 22 depicts the ESLs and NSLs constructed during our local mesh movement procedure. Fig. 23 shows the ESLs created in iteration cycle progress. The figure indicates that there are 19 ESLs in the final achieved active zone. The final active zone consists of 2047 elements and 918 nodes. However, we only need 2 × 19 = 38 integer data to store the active zone in this problem. If one wishes to store all the active zone data, one will need storing 2047 + 918 = 2965 integer data. However, we only need 38/2965 ≈ 1.28% of the entire active zone data to record the active zone data. Additionally, this storage requirement is only 38/17 361 ≈ 0.22%

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 21. The final achieved active zone, test 2.

1013

Fig. 22. Element sub-layers and node sub-lines created during our local mesh movement procedure, test 2.

Fig. 23. The ESL creation (or active zone expansion) history versus the movement iteration cycle, test 2.

of the entire storage required to store all grid element data. As said before, the displacement vector, i.e. δ in Eq. (1), needs to be stored only for nodes inside the active zone instead of all grid nodes. Therefore, this leads to 1 − 918/8921 = 89.7% less memory storage requirement. Comparing with the CDS, we can conclude that the use of RDS would result in up to 89% less memory requirement. Fig. 24 presents the results of two successive estimations of the left and right bounding elements (when l = s) for all the created ESLs. Following the discussion provided on Fig. 14, we have achieved good results almost for all ESLs after their second estimations, see specially Fig. 24(b). Fig. 25 shows the CTR at each movement iteration cycle. According to this figure, the RDS is at least 10 times faster than the CDS in this test. Fig. 26 shows the resulting occurrence of element vertex angle achieved using RDS procedure. They are compared with those obtained using the CDS at the same region, i.e. active zone. As is seen, the vertex angle occurrences are the same for both strategies. Similar to test 1, the RDS reduces the computational time considerably while the quality of the deformed grid remains the same as that of the CDS. We further investigate the case in which two separate but simultaneous boundary deformations occur. There are two ways to treat these two simultaneous deformations. First, we consider both boundary deformations as one boundary deformation and treat it in the manner described in tests 1 and 2. Therefore, we have only one active zone here. Fig. 27(a) shows the resulting active zone for this case. This treatment is justifiable for the cases in which the two deformed segments are close enough to treat them as one segment. Second, it is possible to treat the two deformed boundary segments in two separate active zones. These zones are constructed simultaneously. Fig. 27(b) shows the resulting active zones for this case. As is seen, the two deformed boundary segments are sufficiently far from each other and that the resulting active zones do not cross each other. The RDS has no limit if one wishes to use any of the above two mesh movements. A suitable choice mostly returns to the situation in which the deformed boundary segments occur. Similar strategy can be taken to treat the mesh movements due to more separate local boundary deformations. At the end, the authors would like to remind that the performance reported in these two test cases are only for one step of an unsteady simulation in treating real moving boundary problems. In reality, the users need to repeat the grid movement computations and procedures many times during their simulation steps, which may take a long time to be

1014

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

Fig. 24. The resulting search lengths for bounding elements (l = s), test 2: (a) at the 1st estimation; (b) at the 2nd estimation.

Fig. 25. CTR versus movement iteration cycle, test 2.

Fig. 26. The occurrence of cell vertex angles in the active zone of RDS and similar region in CDS, test 2.

Fig. 27. Two separate but simultaneous deformations and treating them as: (a) a single active zone; or (b) two separate active zones.

fulfilled [1,5,12,13,15,19,24,28,29,31,33,39,43,44,51]. Therefore, the current achieved performances are much higher in real practical problems, where grid computations are very heavy and may need several hours of computations. However, it should be noted that it is very difficult to quantify the current saving considering the overall costs of CFD simulation. It is because the CFD simulations can be performed using a wide range of scheme choices. One scheme may improve the accuracy of solution slightly while the resulting cost of computation is huge. In such cases, the time for local grid movement is almost constant, while the time of CFD simulation increases dramatically. In other words, the local mesh movement saving cannot be readily compared with the overall cost of CFD simulation. In fact, the saving is very case dependent.

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

1015

6. Conclusion A new local mesh movement strategy has been presented to treat the moving boundary or mesh deformation problems very efficiently. The strategy manipulated the cell-based data structure of an original unstructured grid and divided the computational domain into two active and inactive zones. This was achieved via constructing a quasi-structured data structure from an original unstructured data structure. This new structure enabled us to establish a fast partitioning process with a low memory storage requirement. It was shown that the grid computation time would decrease due to confining the grid computations only to the active zone. Additionally, benefiting from the quasi-structured pattern provided for the original unstructured grid data, there was very low memory storage requirement to find the active zone or trace its gradual expansion. The robustness and efficiency of the reduced domain strategy (RDS) were evaluated by treating two common and typical mesh movement and deformation problems. The tests showed that the computational time would be reduced considerably using our strategy. The current study also showed that the RDS was at least five and ten times faster than the complete domain strategy (CDS) in test 1 and test 2, respectively. It was also shown that the RDS uses about 70% less computer memory storages than that of the CDS. The new developed strategy is highly flexible to perform different local mesh movement or deformation tasks and remarkably efficient to eliminate the need for sequential re-ordering of the entire grid data in various grid movement missions. References [1] K. Arens, P. Rentrop, S.O. Stoll, U. Wever, An adjoint approach to optimal design of turbine blades, Appl. Numer. Math. 53 (2) (2005) 93–105. [2] T. Baker, Mesh deformation and reconstruction for time evolving domains, AIAA Paper No. 2001-2535, in: 15th AIAA Computational Fluid Dynamics Conference, Anaheim, California, 2001. [3] J.T. Batina, Unsteady Euler airfoil solutions using unstructured dynamic meshes, AIAA Paper No. 1989-0115, in: 27th Aerospace Sciences Meeting and Exhibit, Reno, NV, 1989. [4] J.D. Bau, H. Luo, R. Loehner, E. Goldberg, A. Feldhun, Application of unstructured moving body methodology to the simulation of fuel tank separation from an F-16 fighter, AIAA Paper No. 1997-0166, in: 35th Aerospace Sciences Meeting and Exhibit, Reno, NV, 1997. [5] O. Baysal, M.E. Eleshaky, Aerodynamic design optimization using sensitivity analysis and computational fluid dynamics, AIAA Paper No. 91-0471, in: 29th Aerospace Sciences Meeting, Reno, NV, 1991. [6] F.J. Blom, Considerations on the spring analogy, Internat. J. Numer. Methods Fluids 32 (6) (2000) 647–668. [7] C.L. Bottasso, D. Detomi, A procedure for tetrahedral boundary layer mesh generation, Eng. Comput. 18 (1) (2002) 66–79. [8] C.L. Bottasso, D. Detomi, R. Serra, The ball-vertex method: a new simple spring analogy method for unstructured dynamic meshes, Comput. Methods Appl. Mech. Engrg. 194 (39–41) (2005) 4244–4264. [9] C.J. Budd, W. Huang, R.D. Russell, Adaptivity with moving grids, Acta Numer. 18 (2009) 111–241. [10] A. de Boer, M.S. van der Schoot, H. Bijl, Mesh deformation based on radial basis function interpolation, Comput. Struct. 85 (11–14) (2007) 784–795. [11] C. Degand, C. Farhat, A three-dimensional torsional spring analogy method for unstructured dynamic meshes, Comput. Struct. 80 (3–4) (2002) 305–316. [12] W. Dettmer, D. Peric, A computational framework for free surface fluid flows accounting for surface tension, Comput. Methods Appl. Mech. Engrg. 195 (23–24) (2006) 3038–3071. [13] H.N.V. Dutt, A.K. Sreekanth, Design of airfoils in incompressible viscous flows by numerical optimization, Comput. Methods Appl. Mech. Engrg. 23 (3) (1980) 355–368. [14] C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional spring for two-dimensional dynamic unstructured fluid meshes, Comput. Methods Appl. Mech. Engrg. 163 (1–4) (1998) 231–245. [15] L. Formaggia, J. Peraire, K. Morgan, Simulation of a store separation using the finite element method, Appl. Math. Model. 12 (2) (1988) 175–181. [16] N. Fouladi, M. Darbandi, G.E. Schneider, Developing a robust ordering-based unstructured moving grid strategy, AIAA Paper No. 2007-1460, in: 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 2007. [17] N. Fouladi, M. Darbandi, G.E. Schneider, A directional renumbering strategy for improving unstructured grid data structure, AIAA Paper No. 2010-1266, in: 48th AIAA Aerospace Sciences Meeting, Orlando, Florida, 2010. [18] M. Grajewski, M. Köster, S. Turek, Numerical analysis and implementational aspects of a new multilevel grid deformation method, Appl. Numer. Math. 60 (8) (2010) 767–781. [19] O. Hassan, E.J. Probert, K. Morgan, N.P. Weatherill, Unsteady flow simulation using unstructured meshes, Comput. Methods Appl. Mech. Engrg. 189 (4) (2000) 1247–1275. [20] B.T. Helenbrook, Mesh deformation using the biharmonic operator, Internat. J. Numer. Methods Engrg. 56 (7) (2003) 1007–1021. [21] W. Huang, R.D. Russell, Adaptive Moving Mesh Methods, Springer, New York, 2011. [22] S. Jakobsson, O. Amoignon, Mesh deformation using radial basis functions for gradient-based aerodynamic shape optimization, Comput. & Fluids 36 (6) (2007) 1119–1136. [23] A.A. Johnson, T.E. Tezduyar, Simulation of multiple spheres falling in a liquid-filled tube, Comput. Methods Appl. Mech. Engrg. 134 (3–4) (1996) 351– 373. [24] P. Le Tallec, J. Mouro, Fluid structure interaction with large structural displacements, Comput. Methods Appl. Mech. Engrg. 190 (24–25) (2001) 3039– 3067. [25] E. Lefrançois, A simple mesh deformation technique for fluid–structure interaction based on a submesh approach, Internat. J. Numer. Methods Engrg. 75 (9) (2008) 1085–1101. [26] X. Liu, N. Qin, H. Xia, Fast dynamic grid deformation based on Delaunay graph mapping, J. Comput. Phys. 211 (2) (2006) 405–423. [27] R. Loehner, C. Yang, Improved ALE mesh velocities for moving bodies, Comm. Numer. Methods Engrg. 12 (10) (1996) 599–608. [28] R. Löhner, C. Yang, E. Onate, S. Idelssohn, An unstructured grid-based, parallel free surface solver, Appl. Numer. Math. 31 (3) (1999) 271–293. [29] A.M. Morris, C.B. Allen, T.C.S. Rendall, Domain element method for aerodynamic shape optimization applied to modern transport wing, AIAA J. 47 (7) (2009) 1647–1659. [30] M. Murayama, K. Nakahashi, K. Matsushima, Unstructured dynamic mesh for large movement and deformation, AIAA Paper No. 2002-0122, in: 40th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 2002. [31] A. Naderi, M. Darbandi, M. Taeibi-Rahni, Developing a unified FVE-ALE approach to solve unsteady fluid flow with moving boundaries, Internat. J. Numer. Methods Fluids 63 (1) (2010) 40–68. [32] E.J. Nielsen, W.K. Anderson, Resent improvements in aerodynamic design optimization on unstructured meshes, AIAA J. 40 (6) (2002) 1155–1163.

1016

M. Darbandi, N. Fouladi / Applied Numerical Mathematics 61 (2011) 1001–1016

[33] A. Oyama, M.S. Liou, S. Obayashi, Transonic axial-flow blade shape optimization using evolutionary algorithms and three-dimensional Navier–Stokes solver, AIAA Paper No. 2002-5642, in: 9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Atlanta, Georgia, 2002. [34] M.A. Potsdam, G.P. Guruswamy, A parallel multiblock mesh movement scheme for complex aeroelastic applications, AIAA Paper No. 2001-0716, in: 39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, 2001. [35] T.C.S. Rendall, C.B. Allen, Unified fluid–structure interpolation and mesh motion using radial basis functions, Internat. J. Numer. Methods Engrg. 74 (10) (2008) 1519–1559. [36] T.C.S. Rendall, C.B. Allen, Efficient mesh motion using radial basis functions with data reduction algorithms, J. Comput. Phys. 228 (17) (2009) 6231– 6249. [37] T.C.S. Rendall, C.B. Allen, Reduced surface point selection options for efficient mesh deformation using radial basis functions, J. Comput. Phys. 229 (8) (2010) 2810–2820. [38] T.C.S. Rendall, C.B. Allen, Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors, Internat. J. Numer. Methods Engrg. 81 (2010) 89–105. [39] D. Sasaki, S. Obayashi, K. Nakahashi, Navier–Stokes optimization of supersonic wings with four objectives using evolutionary algorithm, J. Aircr. 39 (4) (2002) 621–629. [40] T. Tang, Moving mesh methods for computational fluid dynamics, in: Contemp. Math., vol. 383, 2005, pp. 141–174. [41] H.M. Tsai, A.S.F. Wong, J. Cai, Y. Zhu, F. Liu, Unsteady flow calculations with a parallel multiblock moving mesh algorithm, AIAA J. 39 (6) (2001) 1021–1029. [42] L. Tysell, Grid deformation of 3D hybrid grids, in: Proceedings of the 8th International Conference on Numerical Grid Generation in Computational Field Simulations, Honolulu, Hawaii, 2002, pp. 265–274. [43] E.H. van Brummelen, K.G. van der Zee, R. de Borst, Space/time multigrid for a fluid–structure-interaction problem, Appl. Numer. Math. 58 (12) (2008) 1951–1971. [44] A.H. van Zuijlen, A. de Boer, H. Bijl, Higher-order time integration through smooth mesh deformation for 3D fluid–structure interaction simulations, J. Comput. Phys. 224 (1) (2007) 414–430. [45] V. Venkatakrishnan, D.J. Mavriplis, Implicit method for the computation of unsteady flows on unstructured grids, J. Comput. Phys. 127 (2) (1996) 380–397. [46] M. Wu, L. Wang, Y. Song, Preconditioned AOR iterative method for linear systems, Appl. Numer. Math. 57 (5–7) (2007) 672–685. [47] H.Q. Yang, Reduced order method for deforming unstructured grid in CFD, AIAA Paper No. 2010-4617, in: 40th Fluid Dynamics Conference and Exhibit, Chicago, Illinois, 2010. [48] Z. Yang, D.J. Mavriplis, Mesh deformation strategy optimized by the adjoint method on unstructured meshes, AIAA J. 45 (12) (2007) 2885–2896. [49] P.A. Zegeling, in: J.F. Thompson, B.K. Soni, N.P. Weaherill (Eds.), Handbook of Grid Generation, CRC Press, 1999, pp. 37-1–37-18. [50] D. Zeng, C.R. Ethier, A semi-torsional spring analogy model for updating unstructured meshes in 3D moving domains, Finite Elem. Anal. Des. 41 (11–12) (2005) 1118–1139. [51] H. Zhao, J.Z. Wu, J.S. Luo, Turbulent drag reduction by traveling wave of flexible wall, Fluid Dyn. Res. 34 (3) (2004) 175–198.