- Email: [email protected]
The Multigrid Method
Chapter 5
The Multigrid Method 5.1 BASIC CONCEPTS The vast majority of computational fluid dynamics (CFD) techniques utilize some type of iterative or time-stepping approach that requires multiple sweeps through the flow field. The convergence of these techniques can be greatly enhanced by the use of a technique called multigrid (MG). The MG method was originally developed by Fedorenko [1] for elliptic equations and the full potential of it was put forward and systemized by Brandt [2]. Later, Jameson [3] applied it to the Euler equations by solving it on the fine grid and then progressively transfers these flow-field variables and residuals to a series of coarser grids. Subsequently, Jameson employs this method to solve the governing equations implicitly in time so as to avoid the explicit time step limitation [4]. In this book, first, the idea of the MG method is applied to both unsteady compressible and incompressible flow solvers, which perform early iterations on a fine grid and then progressively transfers these flow-field variables and residuals to a series of coarser grids. Thus the solutions on the coarser grids are driven by the fine-grid residuals. The advantage of this transfer and subsequent computation on the courser grids is to cancel the dominating lowfrequency error mode in the late stages of convergence and hence accelerate the convergence rate. The flow is then solved on the coarser grids, and the correction, the difference between the newly computed coarse-grid values and the transferred values, is then interpolated back to the fine grid. The purpose of this correction is to correct the solution on the fine grid, thus eliminating its low-frequency errors. The process is repeated over a sufficient number of times until satisfactory convergence on the fine grid is achieved. For the ease of implementation, the nonnested mesh method using independently generated nonnested coarse meshes is adopted [5]. There are two reasons for using the MG method. One is that the control volumes are larger in the coarse grid so that it permits larger time steps when an explicit MG time-stepping scheme is used to advance the solution. A coarse grid has roughly a fourth to one eighth of the number of points as the finer grid for unstructured grid, so the computational effort per time step is dramatically reduced on the coarse grids. The second reason is to cancel the dominant low-frequency error mode in stages of convergence because a low-frequency Computational Fluid-Structure Interaction. DOI: https://doi.org/10.1016/B978-0-12-814770-2.00005-2 © 2019 Elsevier Inc. All rights reserved.
63
64
Computational Fluid-Structure Interaction
error mode becomes a high frequency one once a fine-grid solution is transferred to a coarser grid. Moreover, the MG method is further developed and thus it can be equally applied to solve the equations of dynamic equilibrium, allowing convergence to be achieved with much less CPU time than the single-grid scheme. For a three-grid system, an MG cycle begins by solving the flow on a fine grid, and then the flow-field variables and residuals are transferred to a coarser grid. After the flow is solved on this coarser grid, it is then transferred to the coarsest grid, where the flow is solved again. The correction is then interpolated up to the fine grid with an additional iteration being performed on the fine grid. In this book,two different cycle strategies have been investigated and presented which are the V-cycle and W-cycle as depicted in Fig. 5.1. For one complete V-cycle, the correction is interpolated back up to the two-finer grids once the coarsest grid is reached, with one new evaluation of the flow-field variables on the fine grid. For one complete W-cycle, the flow is interpolated up one grid after performing the iteration on the coarsest grid. Then, the flow-field variables are evaluated on this grid and the variables are transferred down to the coarsest grid again. The flow is solved on the coarsest grid and the correction is then interpolated up to the fine grid. The total computational work per MG cycle relative to the work on the fine grid for both V-cycle and W-cycle using quadrilateral cells are as follows [6]: V-cycle: WMG ,
4 Wfine ; 3
W-cycle: WMG , 2Wfine
where the subscript MG and fine represent MG cycle and one time-step on the fine grid, respectively. Usually, the number of nodes of the coarse grid can be from a half to an eighth of that of the fine grid. The work associated with grid transfer operations has been neglected, which is negligible in comparison with the work in the time-step.
FIGURE 5.1 V-cycle and W-cycle for three-grid levels E denote flow-field computation on a particular level and I denote interpolation of corrections of flow-field variables.
The Multigrid Method Chapter | 5
65
5.2 MATHEMATICAL FORMULATIONS In this section, a two-level MG is used to describe the formulations. Fine grid, h, denotes the first level, whereas coarse grid, h 1 1, denotes the second level of the MG. The solution on the coarse grid is initialized by transferring the flow-field variables from the fine grid using the transfer operator Thh11 : ð0Þ Wh11 5 T h11 h Wh ð0Þ Wh11
ð5:1Þ T h11 h
where is the initial values transferred from the fine grid, is a transfer operator,which will be defined in Section 5.6, and Wh is the solution from the fine grid. And the initial residual on the coarse grid is defined as ð0Þ ~~ R~~h11 5 Qh11 h RðWh Þ
ð0Þ
ð5:2Þ
is the where R~~h11 is the initial residual transferred from the fine grid, Q h11 h residual transfer operator based on weighted average, which will be defined ~~ in Section 5.6 and RðW h Þ is the residual from the fine grid. The transfer operators in Eqs. (5.1) and (5.2) are restriction operators that may not necessarily be the same. To use fine-grid residuals to drive the solution on the coarser grids, a forcing function is calculated at the first stage of the RungeKutta scheme and subsequently, this forcing term was added to the residual at the coarse grid. The forcing function on the coarse grid is defined as ð0Þ ð0Þ Ph11 5 R~~h11 2 R~~ Wh11 ð5:3Þ It should be noted that on the coarse grids, time-dependent terms in the residual containing W at n and n 2 1 time levels are not included for the ease of calculation. Instead, they are only included in the fine-grid residual and directly transferred to the coarse grids. And every stage of the Runge 2 Kutta scheme in Eq. (2.41) for the compressible flow solver and Eq. (3.58) for the incompressible flow solver on the coarse grid is modified to include the forcing term, Ph11 as follows: i Δτ h ~~ ðm21Þ ðmÞ ð0Þ R Wh11 Wh11 5 Wh11 2 αm C 1 Ph11 ΔScv where m denotes each stage of the RungeKutta scheme. Substituting the ð0Þ ~ ~ forcing term into the above equation, it can be seen that R Wh11 and ðm21Þ R~~ Wh11 almost cancel each other if the reduction of residual is small, ~~ Þ is the only term remains. Thus this will ensure that the and Qh11 RðW h
h
residual on the fine grid drives the coarse-grid solution.
66
Computational Fluid-Structure Interaction
The solution and residuals are transferred onto successive independent coarser grids until the coarsest grid is reached. And after calculating the flow-field variables on the coarsest grid, the correction is evaluated and then has to be interpolated back to the fine grid. The correction is the difference 1 between the newly computed value on the coarse grid, Wh11 , and the initial ð0Þ value that was transferred from the finer grid, Wh11 , and this correction is added to the values on the fine grid as follows: ð0Þ 1 2 Wh11 Wh1 5 Wh 1 I hh11 Wh11 ð5:4Þ where I hh11 is an prolongation operator from the coarse grid to the fine grid, which will be defined in Section 5.6. Wh1 is the updated variables of the solution and Wh is the solution from the fine grid. To speed up the calculations at the coarser grids, the viscous terms are only evaluated on the fine grid and not on the coarser grids. Since the coarser grids are used to cancel the dominating low-frequency errors on the fine grid and the solution on the coarser grids are driven by the residual from the fine grid; hence, this treatment does not affect the accuracy of the solution. The upwind-biased discretization scheme is set to first-order in the coarser levels where the left-right states of Eqs. (3.26a) and (3.26b) are taken as the values of the two nodes of the edge to calculate the flux associated with the edge.
5.3 DATA STRUCTURE FOR DIFFERENT MULTIGRID LEVEL In this book, an efficient and robust data structure and dynamic memory allocation for different numbers of MG levels is employed and fully explained. A two-level MG is used as an example to describe the data structure implemented in the solver. The general idea of this data structure is that the data of finer grids is stacked on top of the data of coarser grids. This requires an indexing array pointer to point to the starting location of a particular MG level array and this indexing array is a 1D array. The indexing array pointer is established by getting the number of data for each MG level, for instance, the numbers of nodes and elements for the fine and coarse grids. And the starting location of the fine grid is one, whereas for the coarser grid it is the number of data for the fine grid plus one. The code below describes the idea of this indexing array pointer: nstot 5 0 nttot 5 0 do mg 5 1, nmg
! initialize total number of nodes for all MG level ! initialize total number of elements for all MG level ! mg is the loop index & nmg is the number of MG level
open(nunit,file 5 grdfill,status 5 'old',form 5 'unformatted')
The Multigrid Method Chapter | 5
67
read(nunit, ) ns(mg),nt(mg) ! ns(mg) and nt(mg) is number of nodes & elements for a particular MG level close(nunit) nsdex(mg) 5 nstot 1 1 ! nsdex(mg) is the indexing array pointer for nodes ntdex(mg) 5 nttot 1 1 ! ntdex(mg) is the indexing array pointer for elements nstot 5 nstot 1 ns(mg) ! compute the total nodes for all MG levels nttot 5 nttot 1 nt(mg) ! compute the total elements for all MG levels enddo
With the total number of nodes and elements for all the MG levels known, those arrays for storing the nodes and elements can be dynamically allocated to the memory in a more efficient and memory saving way. The indexing array pointer is used as the parameters of the counting loop; they control the values of the loop index during the execution of the loop. An example of the use of indexing array pointer for different MG levels is shown as follows: do ip 5 nsdex(mg),nsdex(mg) 1 ns(mg)-1 .......... .......... enddo
where ip is the loop index, mg indicates a particular MG level, nsdex (mg) is the indexing array pointer for number of nodes, and ns (mg) is the number of nodes for that grid level. Likewise, the data structure relating to any parameter, for instance, number of edges or boundary elements, can be established in the same manner.
5.4 ZONING OF FLOW FIELD DOMAIN To enhance the efficiency of the flow solver and reduce preprocessing time, the flow field domain is virtually covered by a grid of cubes, forming a background Cartesian mesh, in which searching for a particular node is only done within the cubes, instead of searching the entire whole flow field. This idea is borrowed from the octree search algorithm [7]. In the current method, the flow solver enables the user to specify the number of rows and columns, to get the desired number of cubes. For illustration purposes, Fig. 5.2 shows a 2D flow field domain being covered by a 4 3 4 grid. Increasing the number of cells or zones will increase the efficiency of the solver, but as for
68
Computational Fluid-Structure Interaction FIGURE 5.2 Zoning of flow field domain.
FIGURE 5.3A Negative mapping produced by the algorithm.
odd-shaped geometry with a lot of curvatures, increasing the number of cells or zones will not enhance the efficiency of the solver further. The idea of this algorithm is that if any node of the triangular (2D)/tetrahedral (3D) cell is mapped onto a particular zone as shown in Fig. 5.3A, then this cell is considered mapped to this zone. Fig. 5.3A shows that the cell is mapped to zones 10, 17, and 18. This algorithm will yield a negative mapping if zone 11 of Fig. 5.3A will not be considered to contain this cell, which is not true. Another negative mapping encountered by this algorithm is that the area of the zones may be smaller than the area of the largest element as shown in Fig. 5.3B since the nodes of the cell is used to map onto a particular zone. Fig. 5.3B shows that the cell is mapped to zones 4, 8, and 28, but actually this cell is supposed to be mapped to zones 4, 5, 6, 7, 8, 10, 11, 12, 18, 19, 20, and 28. Both of the above-mentioned problems will give rise to a negative result where node p will not be mapped to these cells
The Multigrid Method Chapter | 5
69
FIGURE 5.3B Area of the zones smaller than the largest element.
during the search for interconnectivity relationship, which will directly affect the accuracy of the flow solution. To eliminate this negative mapping, the corner coordinates of the zones are used to test whether the cell is mapped to this zone. For example, the corner coordinates of zone 11, (xzmin11, yzmax11), are used to test for this mapping, as shown in Figs. 5.3A and 5.3B. The algorithm used to test for this mapping is depicted in Section 5.5. If the corner coordinates are tested to be within the cell, then this triangular/tetrahedral cell is mapped to this zone. By using this technique, the preprocessing work in term of CPU time can be reduced according to the following equation with the assumption of having the same number of nodes in each zone: Total preprocessing work ðCPU timeÞN
Number of nodes NZ
ð5:5Þ
where NZ represents the number of zones virtually decomposed in the flow domain. The developed algorithm is validated with a three-level 3D MG with the fine mesh consisting of 97,336 nodes and 455,625 elements, coarse mesh 24,389 nodes and 109,760 elements, and coarsest mesh 5,832 nodes and 24,565 elements. Fig. 5.4 shows the correlation between number of zones and CPU time. In Fig. 5.4, N denotes the number of dimensions for a problem, which is three here. For example, 3N means the flow domain is virtually covered by grid of 27 cubes. From the figure, it can be seen that the search algorithm takes 50.38 minutes of CPU time to perform the search for the entire flow domain without any virtual decomposition. With the implementation of this zoning algorithm, the searching time has been reduced substantially to 3.18 minutes for 216 zones. Based on the figure, increasing the number of zones will not further reduce the searching time, instead of increasing it. This occurrence is due to the amount of work involved
70
Computational Fluid-Structure Interaction
50
50.38
CPU Time (mins)
40 33.90
30
20 14.29
10 6.98 3.32
3.67
4.62
5N 4N 6N 7N Number of Zones (N = 3)
8N
9N
3.94
3.18
0 1N
2N
3N
FIGURE 5.4 The correlation between number of zones and CPU time.
increases as the number of zones increase beyond an optimum number of zones, which is 216 zones in this case.
5.4.1 Interconnectiveity Relationship Between Meshes Before the flow-field variables and residuals are transferred from the fine grid to the coarse grids or the corrections are interpolated from the coarse grids back to the fine grid, it is necessary to determine in which coarse cell each fine node is located and vice versa. In this book, the MG method makes use of a sequence of independently generated unstructured grids and the new algorithm used to determine in which coarse cell each fine node is located and vice versa is developed here for the purpose of high efficiency. In the preprocessing stage, finding nodes in tetrahedral cells often benefits from having an imaginary bounding box. The node is first tested against this box and having a flag marked as 0 before the full algorithm is employed to detect the node within the cell. The purposes of enclosing the cell within the box and marked the nodes as 0 are to narrow down the search area for the nodes in the same zone and not to repeat the algorithm for a node marked as 1, thus enhance the efficiency of the flow solver. Basically, the algorithm for interconnectivity relationship between meshes is based on the concept of dot product between two vectors. The unit normal vectors oriented inward from the triangular face, n, are used to detect nodes within the cell and outward pointing unit normal vectors are used to detect nodes within the boundary faces for curved surfaces. The vectors p are located at the center of the triangular face of the cell. The dot product,
The Multigrid Method Chapter | 5
71
(pn)nface, between these two vectors, p and n, must be greater than or equal to zero if the angle, θ, subtended between them is less than 90 degrees and less than zero if the angle is greater than 90 degrees, as depicted in the following equation: ðpn Þ
nface
8 9 < nx = 5 p : fng 5 px py pz : ny 5 px nx 1 py ny 1 pz nz nface ; nface 5 1 2 4 : ; nz
ð5:6Þ
where the superscript nface is the triangular face number of the tetrahedral cell. For ease of clarification and illustration, all the figures in this section will be in 2D. With reference to Fig. 5.5A, the dot product between these two vectors must be positive for the three edges for 2D (four faces for 3D) if the fine node is to fall within the coarse cell or vice versa. If one of the dot products between these two vectors is less than zero, then the node being tested is not within the cell, as depicted in Fig. 5.5B. The equation shows that the criterion for the node to be within the tetrahedral cell is: ðpn Þ1 :and: ðpn Þ2 :and: ðpn Þ3 :and: ðpn Þ4 $ 0 ð5:7Þ To reduce the preprocessing time further, those nodes that fall within a particular cell are marked as 1 and they will not be tested for the next cell in the same or next zone. In the present study, the following paragraph will describe the full algorithm for the interconnectivity relationship between meshes implemented in the 3D solver. This algorithm can be applied to both situations where the fine nodes fall within the coarse cells and coarse nodes fall within the fine cells. For illustration purposes, the algorithm will be described as to the fine
FIGURE 5.5A Node falls within a 2D triangle cell using dot product.
72
Computational Fluid-Structure Interaction FIGURE 5.5B Node does not fall within a 2D triangle cell.
nodes falling within the coarse cells. The algorithm will begin from zone 1 to the end of the zone specified by the user. Firstly, the unit normal vectors oriented inward from the triangular face (edge for 2D), n, are established for the four triangular faces of the coarse cells. The coarse cell is then enclosed within a bounding box in the same zone and a search is performed for those fine nodes that are located in the same zone as the coarse cell. The full algorithm as depicted in Fig. 5.5A for 2D is performed if both the criteria for the bounding box and having a flag of 0 are fulfilled. Then the vectors, p, pointing from centers of the four faces of the cell to the node under consideration are established and eventually the dot product between these two vectors, p.n, for the four faces are computed. If the values obtained for the four dot products are greater than or equal to zero, then the fine node is within the coarse cell and the flag for this node is marked as 1. Another subroutine is called to compute the volumes for the four subtetrahedral cells within the coarse cell with the fine node as the common vertex for the transfer operator algorithms and the volume of the coarse cell. On the other hand, if one of the dot product between these two vectors is less than zero, then the node being tested is not within the cell, as depicted in Fig. 5.5B. After all the fine nodes with a flag of 0 in the box of the same zone are tested against this coarse cell, then the next coarse cell in the list will be tested with those fine nodes having a flag of 0 again and those having a flag of 1 are skipped. The process for the zone will stop if all the fine nodes are tested and marked as 1 or if all the coarse cells are tested whichever one is completed earlier. The same procedure will continue for the next zone until the last zone specified by the user. The same searching procedure will apply for the coarse nodes within the coarsest cells for three-grid levels after the search for all the fine nodes within the coarse cell are performed.
The Multigrid Method Chapter | 5
73
FIGURE 5.6A Fine node within the coarse boundary edge.
The search for nodes within the boundary faces (edges for 2D) will continue after the search for nodes within the cells are performed if the user specifies that the physical wall of the geometry is curved. By specifying whether the wall is flat or curved, the program can save time by skipping the computation of interconnectivity relationship at the boundary of the geometry. Basically, the difference between this algorithm and the algorithm depicted in Fig. 5.5 is that the normal vector, n, is pointing outward at the boundary face instead of pointing inward. This algorithm needs to fulfill two criteria in order for a node to be within the boundary face. The first criterion to be fulfilled is depicted in Fig. 5.6A for 2D, where the dot product between the outward normal vector, n, and the vector, p, pointing from the center of the boundary face must be greater than or equal to zero (i.e.,p.n $ 0). The following equation states first criterion for 3D: 8 9 < nx = ðpn Þ 5 fpg : fng 5 px py pz : ny 5 ðpx nx 1 py ny 1 pz nz Þ $ 0 ð5:8Þ : ; nz To ensure that vector p is not pointing to a node that is very far away from the boundary faces, the boundary faces will be enclosed within a bounding box, which has a dimension of adding both the maximum and minimum x-, y-, and z-coordinates of the three nodes and the maximum length of the three edges that form the boundary face. The second criterion to be fulfilled is depicted in Fig. 5.6B for 2D where the dot product of both p1 .v and p2 .v must be different in sign (i.e., (pv)1 $ 0 and (pv)2 # 0 or vice versa). For 3D implementation, it is more complicating than 2D, whereby it needs to compare the sign for the various vectors as stated in Eqs. (5.9a,b,c) and (5.10a,b,c). To fulfill this criterion, Eqs. (5.9a) to (5.9c) or (5.10a) to (5.10c) must be satisfied together.
74
Computational Fluid-Structure Interaction FIGURE 5.6B Different sign computed between p1 U v and p2 U v.
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 : ; vz
.
0
,
0
ð5:9aÞ ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 > vx > < = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 > : > ; vz 8 9 > vx > < = p3z : vy 5 ðp3x vx 1p3y vy 1p3z vz Þ3 > : > ; vz
.
0
,
0
ð5:9bÞ ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 < vx = 2 p2z : vy 5 p2x vx 1p2y vy 1p2z vz : ; 8 vz 9 < vx = 3 p3z : vy 5 p3x vx 1p3y vy 1p3z vz : ; vz
.
0
,
0
ð5:9cÞ
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 : ; vz
,
0
.
0
ð5:10aÞ
The Multigrid Method Chapter | 5
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p3z : vy 5 ðp3x vx 1p3y vy 1p3z vz Þ3 : ; vz
75
,
0
.
0
ð5:10bÞ 8 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 , 0 : ; 8 vz 9 < vx = p3z : vy 5 ðp3x vx 1 p3y vy 1 p3z vz Þ 3 . 0 : ; vz ð5:10cÞ
The purpose of this second criterion is mainly to ensure that the nodes being tested is within the boundary face instead of testing other node belonging to the neighboring faces, as shown in Fig. 5.6C for 2D. This figure shows that the sign for these two dot products are the same, which both the values of dot product are less than zero (i.e., p1 .v , 0 and p2 .v , 0, having the same sign). If the node being tested fulfills these two criteria, then the node is projected downward onto the boundary edge (2D) and the projected lengths between this projected node and the two edge nodes are computed for the transfer operator alogorithms, as depicted in Figure 5.7. This is to ensure that the residual transfer is conservative. The above-described algorithm can be summarized in the flowchart as shown in Fig. 5.8 and the pseudo code for the algorithm is given after the flowchart.
FIGURE 5.6C Node belonging to the neigboring edge (i.e., ~ p 1 U~ v , 0 and ~ p 2 U~ v , 0).
76
Computational Fluid-Structure Interaction FIGURE 5.7 Projected lengths for the transfer operator algorithms (2D).
! ---------------------------------------------------------------------------------------------------! Pseudo-code for finding nodes enclose within cell using dot product of two vectors ! ---------------------------------------------------------------------------------------------------subroutine containabc variables definition initialize all the arrays do loop for number of multigrid levels (loop 1) initialize flag as 0 for all the nodes do loop for number of zones (loop 2) loop for number of cells within the zone (loop 3) establish inward normal vectors ~ n for the 4 faces of the cell enclose cell within an imaginary bounding box do loop for number of nodes within the zone (loop 4) if (flag node.eq.0 .and. within the bounding box) then establish vectors ~ p for the 4 faces of the cell compute ~ p .~ n for the 4 faces if (~ p .~ n $ 0 for the 4 faces) then flag node 5 1 call subroutine to calculate volumes for the 4 sub-tetrahedral & volume of the cell go back to loop 4 go back to loop 3 go back to loop 2 if (geometry got curvature) call subroutine to test for nodes within boundary faces go back to loop 1 return end
77
The Multigrid Method Chapter | 5
5.5 MESH-TO-MESH TRANSFER OPERATORS Following the two-dimensional approach presented in Ref. [7] for data transfer within the domain and incorporate the new technique developed here for a curved boundary, there are two classes of mesh-to-mesh transfer operators being implemented in the method for three-dimensions, which are restriction operators and prolongation operators. List of nodes and cells within the specified zones for three-grid levels
Fine-to-coarse grids (1--2 grid level – Group 1) Coarse-to-coarsest grids (2--3 grid level – Group 2)
Initialize the flag as 0 for all the nodes
List of cells in the array for the zone Go to the next group of grid levels (eg., 2--3 grid level for three-grid levels)
Establish the inward normal vectors, n, for the 4 faces of the cells
Enclose the cell within an imaginary bounding box
Is the node flag 0 and within the bounding box?
No
Yes Establish the vectors, p, for the 4 faces of the cell
Compute p .n for the 4 faces
Is p .n ≥ 0 for the 4
No
Yes
Go to the next zone
Go to the next in the list for the zone
Go to the next cell in the list for the zone
List of nodes in the array for the zone
(1) Compute the volumes for the 4 subtetrahedral within the cell and the volume of the cell. (2) Flag the node as 1
A
B
C
D
E
FIGURE 5.8 Flowchart for interconnectivity relationships between meshes algorithm for threegrid levels.
78
Computational Fluid-Structure Interaction
A
B
C
D
End of the list of nodes for the zone?
E
No
Yes No
End of the list of cells for the zone?
Yes End of zones specified by user?
No
Yes No Curvature in the geometry?
Yes Call subroutine to test for nodes that contain within the boundary faces using algorithm depicted in Fig. 5.6 Yes No
End of three-grid levels?
Yes End of algorithm
FIGURE 5.8 (Continued).
5.5.1 Flow-Field Variables Transfer Operators Once the interconnectivity relationship between meshes has been determined, the restriction transfer operator, Thh11 , that transfers the flow-field values from the fine grid to the coarse grid will be as follows: W1 5
Va Wa 1 Vb Wb 1 Vc Wc 1 Vd Wd Va 1 Vb 1 Vc 1 V d
ð5:11Þ
Fig. 5.9A depicts the operation of this transfer operator for those nodes within a cell. Lower case letters denote fine-grid nodes and arabic numbers
The Multigrid Method Chapter | 5
79
FIGURE 5.9A Transfer of flow-field values from the fine mesh to the coarse mesh using restriction transfer operator.
FIGURE 5.9B Transfer of variables from fine nodes to the coarse node at the boundary.
denote coarse-grid nodes for all the figures in this section. The flow-field values at the coarse node 1, which is contained in the fine cell formed by nodes a, b, c, and d, is a weighted average of the values at those nodes. Va is the volume of the sub tetrahedron cell with vertices b, 1, c, and d. Vb is the volume of the sub tetrahedron cell with vertices a, 1, c, and d. Vc is the sub tetrahedron cell with vertices a, 1, d, and b. Vd is the subtetrahedron cell with vertices a, 1, b, and c. The volume of the corresponding tetrahedron cell is the one opposite to the node. According to Eq. (5.11), if node 1 coincides with node a, then W1 will be equal to Wa. At the boundary with curvature, the flow-field values are transferred from the fine nodes of the triangle formed by vertices a, b, and c onto the coarse node 1 using an area-weighted contribution calculation, as illustrated in Fig. 5.9B. The restriction transfer operator, Thh11 , that transfers the flow field values is as follows: W1 5
Aa W a 1 A b W b 1 Ac W c Aa 1 Ab 1 Ac
ð5:12Þ
Fig. 5.9B illustrates the operation of this transfer operator for those nodes within the boundary triangles. The flow-field values at the coarse node 1,
80
Computational Fluid-Structure Interaction FIGURE 5.10 An illustration of transfer of residuals using prolongation operator will not be conserved.
which is projected downward onto the fine triangle face formed by nodes a, b, and c, is a weighted average of the values at those nodes. Aa is the area of the triangle with vertices b, 1, and c. Ab is the area of the triangle with vertices c, 1, and a. Ac is the area of the triangle with vertices a, 1, and b. The area of the corresponding triangle is the one opposite to the node. Similarly, according to Eq. (5.12), if node 1 coincides with node a, then W1 will be equal to Wa.
5.5.2 Residual Transfer Operators In the current MG method, the restriction transfer operator is used for transferring residuals so that the transferred residuals will be conserved. Fig. 5.10 illustrates the reason why using prolongation operators to transfer residuals will not be conserved. For illustration purposes, Fig. 5.10 is illustrated with a 2D mesh. The coarse-grid nodes 1, 2, and 3 are super-positioned onto the three fine cells as shown in the figure. The first fine cell is formed by node a, b, and f, the second fine cell is formed by node i, j, and n, and the third fine cell is formed by node l, m, and q. If the residuals are transferred using prolongation operators, then the residuals at these nine fine nodes will be transferred to the three coarse nodes, while the residual at fine node k will not be transferred to any of the coarse nodes. It is possible to encounter a situation whereby the residuals of these nine fine nodes are zero and the residual of fine node k is nonzero. A negative transfer will result in the three coarse nodes to have zero residual if using prolongation operators. To ensure that the overall convergence will be conserved and to avoid negative transferring between meshes, the selection of the class of operators to be used in transferring the residuals is very important to ensure that the sum of coarsegrid residuals after the transfer are always conserved, that is equal to the total of the fine-grid residuals. where the transfer of residuals from the fine grid to the coarse grid is based on a weighted average of the residual from the fine node to the four coarse nodes that form the coarse tetrahedral cell. As depicted in Fig. 5.11A,
The Multigrid Method Chapter | 5
81
FIGURE 5.11A Transfer of residuals from the fine mesh to the coarse mesh using restriction transfer operator.
the residual from node a is transferred to coarse nodes 1, 2, 3, and 4 by the transfer operator, Q h11 h . The transferred residual to a particular coarse node is the summation of all the transferred weighted-average residuals from other fine nodes in those tetrahedra surrounding this coarse node, as shown in the following equations, where the contribution from fine-grid node a is calculated and added to the residual of coarse-grid node 1: ~~ 5 R ~~ 1 R 1 1 ~~ 5 R ~~ 1 R 3 3
V1 R~~ a ; V1 1 V2 1 V3 1 V4
V3 R~~ a ; V 1 1 V 2 1 V 3 1 V4
~~ 1 R~~ 2 5 R 2 R~~ 4 5 R~~ 4 1
~~ V2 R a ; V1 1 V2 1 V 3 1 V 4
~~ V4 R a V1 1 V2 1 V3 1 V4
ð5:13Þ
where V1 is the volume of the subtetrahedron cell with vertices 2, a, 3, and 4. V2 is the volume of the subtetrahedron cell with vertices 1, a, 3, and 4. V3 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 4. V4 is the volume of the sub tetrahedron cell with vertices 1, a, 2, and 3. The volume of the corresponding tetrahedron cell is the one opposite to node. It is easy to show that this transfer is conservative in the sense that the total fine mesh residuals are equal to the coarse mesh ones. The residual for a find-grid node on the boundary face of the curvature geometry is transferred to the coarse nodes of the triangle formed by vertices 1, 2, and 3 using an area-weighted contribution calculation as shown in Fig. 5.11B. The transfer operator, Q h11 h , that transfers the residuals from the fine grid to the three coarse-grid nodes is as follows: ~~ 5 R ~~ 1 R 1 1
A1 R~~ a A 1 1 A 2 1 A3
82
Computational Fluid-Structure Interaction FIGURE 5.11B Transfer of residuals from fine node to the coarse nodes at the boundary.
R~~ 2 5 R~~ 2 1
~~ A2 R a ; A1 1 A2 1 A3
~~ 5 R ~~ 1 R 3 3
A3 R~~ a A1 1 A2 1 A3
ð5:14Þ
Fig. 5.11B depicts the operation of this transfer operator typically for those nodes on the boundary faces. The residual at the fine node a, which is projected downward onto the coarse triangle formed by nodes 1, 2, and 3, is a weighted average of the value at this fine node. The transferred residual is the summation of all the transferred residuals from other fine nodes to these coarse nodes. This also ensures that the residual transfer is conservative on the boundary. It has been proven that the total fine-grid residuals are equal to the coarse grid ones and this is achieved by performing a summation of the residuals on both fine and coarse grids and the totals are equal to each other [7]. This proves that this restriction transfer operator is conservative.
5.5.3 Correction Transfer Operators Prolongation operators are used to transfer corrections of the flow-field variables from the coarse mesh to the fine mesh, as depicted in Figs. 5.12A and 5.12B. As shown in Eq. (5.15), the correction, dWh11 , is the difference 1 between the newly computed value on the coarse grid, Wh11 , and the initial ð0Þ value that was transferred from the fine grid, Wh11 . ð0Þ 1 dWh11 5 Wh11 2 Wh11
ð5:15Þ
Corrections are transferred to the fine mesh by the prolongation operator, I
h h11 : h dWh11 vh 5 Ih11
and
Wh1 5 Wh 1 vh
According to Fig. 5.12A, the correction of the flow-field variables transferred from the coarse nodes 1, 2, 3, and 4 to the fine node a is a weighted
The Multigrid Method Chapter | 5
83
FIGURE 5.12A Transfer of corrections from the coarse mesh to the fine mesh using prolongation transfer operator.
FIGURE 5.12B Transfer of corrections from the coarse nodes to the fine node at the boundary.
average of the corrections at these nodes and the expression for the transferred correction is ðva Þh 5
V1 ðdW1 Þh11 1 V2 ðdW2 Þh11 1 V3 ðdW3 Þh11 1 V4 ðdW4 Þh11 V1 1 V2 1 V3 1 V4
ð5:16Þ
where V1 is the volume of the subtetrahedron cell with vertices 2, a, 3, and 5. V2 is the volume of the subtetrahedron cell with vertices 1, a, 3, and 5. V3 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 5. V4 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 3. The volume of the corresponding tetrahedron cell is the one opposite to the node. The corrections of the flow-field variables at the boundary are transferred from the coarse nodes of the triangle formed by vertices 1, 2, and 3 onto the fine node a. The transfer operator, I hh11 , that transfers the corrections is ðva Þh 5
A1 ðdW1 Þh11 1 A2 ðdW2 Þh11 1 A3 ðdW3 Þh11 A 1 1 A2 1 A3
ð5:17Þ
84
Computational Fluid-Structure Interaction
Fig. 5.12B shows the operation of this transfer operator for those nodes on the boundary faces. The correction at the fine node a, which is projected downward onto the coarse face formed by nodes 1, 2, and 3, is a weighted average of the values at those nodes. A1 is the area of the triangle with vertices 2, a, and 3. A2 is the area of the triangle with vertices 1, a, and 3. A3 is the area of the triangle with vertices 1, a, and 2. The area of the corresponding triangle is the one opposite to the node. Similarly, according to Eq. (5.17), if node a coincides with node 1, then ðva Þh will be equal to ðdW1 Þh11 .
5.6 THE MULTIGRID METHOD FOR STRUCTURAL DYNAMIC SOLVER In this section, an efficient unstructured MG scheme developed for the equation of dynamic equilibrium is introduced. The MG algorithm is described as follows. The discretized dynamic equation can be expressed by equation (4.25), which is repeated here for convenience: ΔU ~~ n11;m 5R Δτ
ð5:18Þ
where U is the velocity vector used in CSD solver. This equation is solved iteratively by a dual-time-stepping scheme. An pseudo code is given below to illustrate the basic procedures of the MG scheme. The outer cycles are based on the physical time t, which is numbered from 1 to ktmax, whereas the inner cycles are based on the pseudo time τ, which is numbered from 1 to max_itr_sub. The grid levels range from 1 to nmg, where 1 is the finest level and nmg is the coarsest level. Phh11 is the prolongation operator from level h 1 1 to h, and Qh11 and Thh11 are the residual transfer and restriction h operators from level h to h 1 1. ALGORITHM START call setupinterconnect !---. build up inter-connectivity relationship between levels DO kt 5 1, ktmax !---. start of physical time step DO mg 5 1, nmg call cvvol(mg) !---. calculate the volumes for control volume (CV) and cell call clhaut(mg) !---. calculate CV characteristic length for determination !---. of Δτ and Δt computation ENDDO call dtsize(mg 5 1) !---. computation of the time-step size for the finest level DO itersub 5 1, max_itr_sub !---. start of sub-iteration
The Multigrid Method Chapter | 5
85
DO ialpha 5 1, 5 !---. 5-stage Runge-Kutta time integration process kt call fvsolver(mg 5 1) !---. solve equation (3.1) for ~ U1 following the way !--- . described in forgoing sections ENDDO DO mg 5 2, nmg call transfrsol(mg-1) !---. h to h 1 1 ð0Þ U h ). !---. (~ U h11 5 Thh11 ~ call transfresd(mg-1) !---. h to h 1 1 ~ ~ ð0Þ R~ U !---. (R~ h11 5 Q h11 h h ).
restrict
solution
from
restrict
residual
from
IF (itersub.eq.1) THEN call dtsize(mg) !---. calculate the time-step size for the current level ENDIF DO ialpha 5 1, 5 kt call fvsolver(mg) !---. solve equation (3.1) for ~ U mg based on the initial ð0Þ ~ ð0Þ !---. solution ~ U h11 and initial residual R~ h11 ENDDO ENDDO DO mg 5 nmg, 2, -1 call transfrcorct(mg) !---. prolongate correction from h 1 1 to h 1 1 ð0Þ 1 !---. (~ Uh 5 ~ U h 1 Ihh11 ~ U h11 ). Here ~ U h11 2 ~ U h is !---. the updated solution for the finer grids ENDDO ENDDO !---. end of sub-iteration DO mg 5 2, nmg call transfrsol(mg-1) !---. restrict solution from h to h 1 1 ENDDO DO mg 5 1, nmg call updgrid(mg) !---. finally update the solid mesh ENDDO ENDDO !---. end of physical time step ALGORITHM END
One should note that we always use the initial intermesh connectivity, which is built up before the first time-step starts, to perform the mesh-tomesh transfer operations. In order to enable the solver to tackle more geometrically complex structures, the search for nodes within the boundary faces (edges for 2D) will continue after the search for nodes within the cells are performed if the user specifies that the physical wall of the geometry is
86
Computational Fluid-Structure Interaction Coarser-grid node
Coarser-grid boundary
Finer-grid node
Finer-grid boundary
Projection node for P
FIGURE 5.13 Schematic of boundary node projection algorithm, Ω is the computational domain.
Node P p1 1
θ1 θ 2
p2 2
V Ω
curved (e.g., see Fig. 5.13 ), a normal procedure cannot find a corresponding finer grid cell for the coarser grid node P. It means that the node P cannot get a proper solution vector from the finer grid, which can pose a serious problem while solving the dynamic equation (4.28). Following the boundary node projection algorithm proposed in Section 5.6, we resort to finding a projection node for node P in the nearest finer mesh cell face and then use this projected node to perform the necessary MG operations.
REFERENCES [1] R.P. Fedorenko, The speed of convergence of one iterative process, U.S.S.R., Comput. Math. Math. Phys. 4 (3) (1964) 227235. [2] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput. 31 (1977) 333390. [3] A. Jameson, Solution of the Euler equations for two dimensional transonic flow by a multigrid method, Appl. Math.Comput 13 (1983) 327355. [4] Jameson A., Time Dependent Calculations Using Multigrid, with Applications to Unsteady flows Past Airfoils and Wings, AIAA Paper 91-1596, June 1991. [5] E. Turkel, Preconditioned Methods for Solving the Incompressible and Low Speed Compressible Equations, J. Comput. Physics 72 (1987). [6] Mavriplis D.J., Large-scale Parallel Viscous Flow Computations Using An Unstructured Multigrid Algorithm, NASA/CR-1999-209724, ICASE Report No. 99-44, Nov 1999. [7] C.H. Tai, Y. Zhao, A Finite Volume Unstructured Multigrid Method for Efficient Computation of Unsteady Incompressible Viscous Flows, Int. J. Num. Meth. Fluids 46 (1) (2004) 5984.
The Multigrid Method 5.1 BASIC CONCEPTS The vast majority of computational fluid dynamics (CFD) techniques utilize some type of iterative or time-stepping approach that requires multiple sweeps through the flow field. The convergence of these techniques can be greatly enhanced by the use of a technique called multigrid (MG). The MG method was originally developed by Fedorenko [1] for elliptic equations and the full potential of it was put forward and systemized by Brandt [2]. Later, Jameson [3] applied it to the Euler equations by solving it on the fine grid and then progressively transfers these flow-field variables and residuals to a series of coarser grids. Subsequently, Jameson employs this method to solve the governing equations implicitly in time so as to avoid the explicit time step limitation [4]. In this book, first, the idea of the MG method is applied to both unsteady compressible and incompressible flow solvers, which perform early iterations on a fine grid and then progressively transfers these flow-field variables and residuals to a series of coarser grids. Thus the solutions on the coarser grids are driven by the fine-grid residuals. The advantage of this transfer and subsequent computation on the courser grids is to cancel the dominating lowfrequency error mode in the late stages of convergence and hence accelerate the convergence rate. The flow is then solved on the coarser grids, and the correction, the difference between the newly computed coarse-grid values and the transferred values, is then interpolated back to the fine grid. The purpose of this correction is to correct the solution on the fine grid, thus eliminating its low-frequency errors. The process is repeated over a sufficient number of times until satisfactory convergence on the fine grid is achieved. For the ease of implementation, the nonnested mesh method using independently generated nonnested coarse meshes is adopted [5]. There are two reasons for using the MG method. One is that the control volumes are larger in the coarse grid so that it permits larger time steps when an explicit MG time-stepping scheme is used to advance the solution. A coarse grid has roughly a fourth to one eighth of the number of points as the finer grid for unstructured grid, so the computational effort per time step is dramatically reduced on the coarse grids. The second reason is to cancel the dominant low-frequency error mode in stages of convergence because a low-frequency Computational Fluid-Structure Interaction. DOI: https://doi.org/10.1016/B978-0-12-814770-2.00005-2 © 2019 Elsevier Inc. All rights reserved.
63
64
Computational Fluid-Structure Interaction
error mode becomes a high frequency one once a fine-grid solution is transferred to a coarser grid. Moreover, the MG method is further developed and thus it can be equally applied to solve the equations of dynamic equilibrium, allowing convergence to be achieved with much less CPU time than the single-grid scheme. For a three-grid system, an MG cycle begins by solving the flow on a fine grid, and then the flow-field variables and residuals are transferred to a coarser grid. After the flow is solved on this coarser grid, it is then transferred to the coarsest grid, where the flow is solved again. The correction is then interpolated up to the fine grid with an additional iteration being performed on the fine grid. In this book,two different cycle strategies have been investigated and presented which are the V-cycle and W-cycle as depicted in Fig. 5.1. For one complete V-cycle, the correction is interpolated back up to the two-finer grids once the coarsest grid is reached, with one new evaluation of the flow-field variables on the fine grid. For one complete W-cycle, the flow is interpolated up one grid after performing the iteration on the coarsest grid. Then, the flow-field variables are evaluated on this grid and the variables are transferred down to the coarsest grid again. The flow is solved on the coarsest grid and the correction is then interpolated up to the fine grid. The total computational work per MG cycle relative to the work on the fine grid for both V-cycle and W-cycle using quadrilateral cells are as follows [6]: V-cycle: WMG ,
4 Wfine ; 3
W-cycle: WMG , 2Wfine
where the subscript MG and fine represent MG cycle and one time-step on the fine grid, respectively. Usually, the number of nodes of the coarse grid can be from a half to an eighth of that of the fine grid. The work associated with grid transfer operations has been neglected, which is negligible in comparison with the work in the time-step.
FIGURE 5.1 V-cycle and W-cycle for three-grid levels E denote flow-field computation on a particular level and I denote interpolation of corrections of flow-field variables.
The Multigrid Method Chapter | 5
65
5.2 MATHEMATICAL FORMULATIONS In this section, a two-level MG is used to describe the formulations. Fine grid, h, denotes the first level, whereas coarse grid, h 1 1, denotes the second level of the MG. The solution on the coarse grid is initialized by transferring the flow-field variables from the fine grid using the transfer operator Thh11 : ð0Þ Wh11 5 T h11 h Wh ð0Þ Wh11
ð5:1Þ T h11 h
where is the initial values transferred from the fine grid, is a transfer operator,which will be defined in Section 5.6, and Wh is the solution from the fine grid. And the initial residual on the coarse grid is defined as ð0Þ ~~ R~~h11 5 Qh11 h RðWh Þ
ð0Þ
ð5:2Þ
is the where R~~h11 is the initial residual transferred from the fine grid, Q h11 h residual transfer operator based on weighted average, which will be defined ~~ in Section 5.6 and RðW h Þ is the residual from the fine grid. The transfer operators in Eqs. (5.1) and (5.2) are restriction operators that may not necessarily be the same. To use fine-grid residuals to drive the solution on the coarser grids, a forcing function is calculated at the first stage of the RungeKutta scheme and subsequently, this forcing term was added to the residual at the coarse grid. The forcing function on the coarse grid is defined as ð0Þ ð0Þ Ph11 5 R~~h11 2 R~~ Wh11 ð5:3Þ It should be noted that on the coarse grids, time-dependent terms in the residual containing W at n and n 2 1 time levels are not included for the ease of calculation. Instead, they are only included in the fine-grid residual and directly transferred to the coarse grids. And every stage of the Runge 2 Kutta scheme in Eq. (2.41) for the compressible flow solver and Eq. (3.58) for the incompressible flow solver on the coarse grid is modified to include the forcing term, Ph11 as follows: i Δτ h ~~ ðm21Þ ðmÞ ð0Þ R Wh11 Wh11 5 Wh11 2 αm C 1 Ph11 ΔScv where m denotes each stage of the RungeKutta scheme. Substituting the ð0Þ ~ ~ forcing term into the above equation, it can be seen that R Wh11 and ðm21Þ R~~ Wh11 almost cancel each other if the reduction of residual is small, ~~ Þ is the only term remains. Thus this will ensure that the and Qh11 RðW h
h
residual on the fine grid drives the coarse-grid solution.
66
Computational Fluid-Structure Interaction
The solution and residuals are transferred onto successive independent coarser grids until the coarsest grid is reached. And after calculating the flow-field variables on the coarsest grid, the correction is evaluated and then has to be interpolated back to the fine grid. The correction is the difference 1 between the newly computed value on the coarse grid, Wh11 , and the initial ð0Þ value that was transferred from the finer grid, Wh11 , and this correction is added to the values on the fine grid as follows: ð0Þ 1 2 Wh11 Wh1 5 Wh 1 I hh11 Wh11 ð5:4Þ where I hh11 is an prolongation operator from the coarse grid to the fine grid, which will be defined in Section 5.6. Wh1 is the updated variables of the solution and Wh is the solution from the fine grid. To speed up the calculations at the coarser grids, the viscous terms are only evaluated on the fine grid and not on the coarser grids. Since the coarser grids are used to cancel the dominating low-frequency errors on the fine grid and the solution on the coarser grids are driven by the residual from the fine grid; hence, this treatment does not affect the accuracy of the solution. The upwind-biased discretization scheme is set to first-order in the coarser levels where the left-right states of Eqs. (3.26a) and (3.26b) are taken as the values of the two nodes of the edge to calculate the flux associated with the edge.
5.3 DATA STRUCTURE FOR DIFFERENT MULTIGRID LEVEL In this book, an efficient and robust data structure and dynamic memory allocation for different numbers of MG levels is employed and fully explained. A two-level MG is used as an example to describe the data structure implemented in the solver. The general idea of this data structure is that the data of finer grids is stacked on top of the data of coarser grids. This requires an indexing array pointer to point to the starting location of a particular MG level array and this indexing array is a 1D array. The indexing array pointer is established by getting the number of data for each MG level, for instance, the numbers of nodes and elements for the fine and coarse grids. And the starting location of the fine grid is one, whereas for the coarser grid it is the number of data for the fine grid plus one. The code below describes the idea of this indexing array pointer: nstot 5 0 nttot 5 0 do mg 5 1, nmg
! initialize total number of nodes for all MG level ! initialize total number of elements for all MG level ! mg is the loop index & nmg is the number of MG level
open(nunit,file 5 grdfill,status 5 'old',form 5 'unformatted')
The Multigrid Method Chapter | 5
67
read(nunit, ) ns(mg),nt(mg) ! ns(mg) and nt(mg) is number of nodes & elements for a particular MG level close(nunit) nsdex(mg) 5 nstot 1 1 ! nsdex(mg) is the indexing array pointer for nodes ntdex(mg) 5 nttot 1 1 ! ntdex(mg) is the indexing array pointer for elements nstot 5 nstot 1 ns(mg) ! compute the total nodes for all MG levels nttot 5 nttot 1 nt(mg) ! compute the total elements for all MG levels enddo
With the total number of nodes and elements for all the MG levels known, those arrays for storing the nodes and elements can be dynamically allocated to the memory in a more efficient and memory saving way. The indexing array pointer is used as the parameters of the counting loop; they control the values of the loop index during the execution of the loop. An example of the use of indexing array pointer for different MG levels is shown as follows: do ip 5 nsdex(mg),nsdex(mg) 1 ns(mg)-1 .......... .......... enddo
where ip is the loop index, mg indicates a particular MG level, nsdex (mg) is the indexing array pointer for number of nodes, and ns (mg) is the number of nodes for that grid level. Likewise, the data structure relating to any parameter, for instance, number of edges or boundary elements, can be established in the same manner.
5.4 ZONING OF FLOW FIELD DOMAIN To enhance the efficiency of the flow solver and reduce preprocessing time, the flow field domain is virtually covered by a grid of cubes, forming a background Cartesian mesh, in which searching for a particular node is only done within the cubes, instead of searching the entire whole flow field. This idea is borrowed from the octree search algorithm [7]. In the current method, the flow solver enables the user to specify the number of rows and columns, to get the desired number of cubes. For illustration purposes, Fig. 5.2 shows a 2D flow field domain being covered by a 4 3 4 grid. Increasing the number of cells or zones will increase the efficiency of the solver, but as for
68
Computational Fluid-Structure Interaction FIGURE 5.2 Zoning of flow field domain.
FIGURE 5.3A Negative mapping produced by the algorithm.
odd-shaped geometry with a lot of curvatures, increasing the number of cells or zones will not enhance the efficiency of the solver further. The idea of this algorithm is that if any node of the triangular (2D)/tetrahedral (3D) cell is mapped onto a particular zone as shown in Fig. 5.3A, then this cell is considered mapped to this zone. Fig. 5.3A shows that the cell is mapped to zones 10, 17, and 18. This algorithm will yield a negative mapping if zone 11 of Fig. 5.3A will not be considered to contain this cell, which is not true. Another negative mapping encountered by this algorithm is that the area of the zones may be smaller than the area of the largest element as shown in Fig. 5.3B since the nodes of the cell is used to map onto a particular zone. Fig. 5.3B shows that the cell is mapped to zones 4, 8, and 28, but actually this cell is supposed to be mapped to zones 4, 5, 6, 7, 8, 10, 11, 12, 18, 19, 20, and 28. Both of the above-mentioned problems will give rise to a negative result where node p will not be mapped to these cells
The Multigrid Method Chapter | 5
69
FIGURE 5.3B Area of the zones smaller than the largest element.
during the search for interconnectivity relationship, which will directly affect the accuracy of the flow solution. To eliminate this negative mapping, the corner coordinates of the zones are used to test whether the cell is mapped to this zone. For example, the corner coordinates of zone 11, (xzmin11, yzmax11), are used to test for this mapping, as shown in Figs. 5.3A and 5.3B. The algorithm used to test for this mapping is depicted in Section 5.5. If the corner coordinates are tested to be within the cell, then this triangular/tetrahedral cell is mapped to this zone. By using this technique, the preprocessing work in term of CPU time can be reduced according to the following equation with the assumption of having the same number of nodes in each zone: Total preprocessing work ðCPU timeÞN
Number of nodes NZ
ð5:5Þ
where NZ represents the number of zones virtually decomposed in the flow domain. The developed algorithm is validated with a three-level 3D MG with the fine mesh consisting of 97,336 nodes and 455,625 elements, coarse mesh 24,389 nodes and 109,760 elements, and coarsest mesh 5,832 nodes and 24,565 elements. Fig. 5.4 shows the correlation between number of zones and CPU time. In Fig. 5.4, N denotes the number of dimensions for a problem, which is three here. For example, 3N means the flow domain is virtually covered by grid of 27 cubes. From the figure, it can be seen that the search algorithm takes 50.38 minutes of CPU time to perform the search for the entire flow domain without any virtual decomposition. With the implementation of this zoning algorithm, the searching time has been reduced substantially to 3.18 minutes for 216 zones. Based on the figure, increasing the number of zones will not further reduce the searching time, instead of increasing it. This occurrence is due to the amount of work involved
70
Computational Fluid-Structure Interaction
50
50.38
CPU Time (mins)
40 33.90
30
20 14.29
10 6.98 3.32
3.67
4.62
5N 4N 6N 7N Number of Zones (N = 3)
8N
9N
3.94
3.18
0 1N
2N
3N
FIGURE 5.4 The correlation between number of zones and CPU time.
increases as the number of zones increase beyond an optimum number of zones, which is 216 zones in this case.
5.4.1 Interconnectiveity Relationship Between Meshes Before the flow-field variables and residuals are transferred from the fine grid to the coarse grids or the corrections are interpolated from the coarse grids back to the fine grid, it is necessary to determine in which coarse cell each fine node is located and vice versa. In this book, the MG method makes use of a sequence of independently generated unstructured grids and the new algorithm used to determine in which coarse cell each fine node is located and vice versa is developed here for the purpose of high efficiency. In the preprocessing stage, finding nodes in tetrahedral cells often benefits from having an imaginary bounding box. The node is first tested against this box and having a flag marked as 0 before the full algorithm is employed to detect the node within the cell. The purposes of enclosing the cell within the box and marked the nodes as 0 are to narrow down the search area for the nodes in the same zone and not to repeat the algorithm for a node marked as 1, thus enhance the efficiency of the flow solver. Basically, the algorithm for interconnectivity relationship between meshes is based on the concept of dot product between two vectors. The unit normal vectors oriented inward from the triangular face, n, are used to detect nodes within the cell and outward pointing unit normal vectors are used to detect nodes within the boundary faces for curved surfaces. The vectors p are located at the center of the triangular face of the cell. The dot product,
The Multigrid Method Chapter | 5
71
(pn)nface, between these two vectors, p and n, must be greater than or equal to zero if the angle, θ, subtended between them is less than 90 degrees and less than zero if the angle is greater than 90 degrees, as depicted in the following equation: ðpn Þ
nface
8 9 < nx = 5 p : fng 5 px py pz : ny 5 px nx 1 py ny 1 pz nz nface ; nface 5 1 2 4 : ; nz
ð5:6Þ
where the superscript nface is the triangular face number of the tetrahedral cell. For ease of clarification and illustration, all the figures in this section will be in 2D. With reference to Fig. 5.5A, the dot product between these two vectors must be positive for the three edges for 2D (four faces for 3D) if the fine node is to fall within the coarse cell or vice versa. If one of the dot products between these two vectors is less than zero, then the node being tested is not within the cell, as depicted in Fig. 5.5B. The equation shows that the criterion for the node to be within the tetrahedral cell is: ðpn Þ1 :and: ðpn Þ2 :and: ðpn Þ3 :and: ðpn Þ4 $ 0 ð5:7Þ To reduce the preprocessing time further, those nodes that fall within a particular cell are marked as 1 and they will not be tested for the next cell in the same or next zone. In the present study, the following paragraph will describe the full algorithm for the interconnectivity relationship between meshes implemented in the 3D solver. This algorithm can be applied to both situations where the fine nodes fall within the coarse cells and coarse nodes fall within the fine cells. For illustration purposes, the algorithm will be described as to the fine
FIGURE 5.5A Node falls within a 2D triangle cell using dot product.
72
Computational Fluid-Structure Interaction FIGURE 5.5B Node does not fall within a 2D triangle cell.
nodes falling within the coarse cells. The algorithm will begin from zone 1 to the end of the zone specified by the user. Firstly, the unit normal vectors oriented inward from the triangular face (edge for 2D), n, are established for the four triangular faces of the coarse cells. The coarse cell is then enclosed within a bounding box in the same zone and a search is performed for those fine nodes that are located in the same zone as the coarse cell. The full algorithm as depicted in Fig. 5.5A for 2D is performed if both the criteria for the bounding box and having a flag of 0 are fulfilled. Then the vectors, p, pointing from centers of the four faces of the cell to the node under consideration are established and eventually the dot product between these two vectors, p.n, for the four faces are computed. If the values obtained for the four dot products are greater than or equal to zero, then the fine node is within the coarse cell and the flag for this node is marked as 1. Another subroutine is called to compute the volumes for the four subtetrahedral cells within the coarse cell with the fine node as the common vertex for the transfer operator algorithms and the volume of the coarse cell. On the other hand, if one of the dot product between these two vectors is less than zero, then the node being tested is not within the cell, as depicted in Fig. 5.5B. After all the fine nodes with a flag of 0 in the box of the same zone are tested against this coarse cell, then the next coarse cell in the list will be tested with those fine nodes having a flag of 0 again and those having a flag of 1 are skipped. The process for the zone will stop if all the fine nodes are tested and marked as 1 or if all the coarse cells are tested whichever one is completed earlier. The same procedure will continue for the next zone until the last zone specified by the user. The same searching procedure will apply for the coarse nodes within the coarsest cells for three-grid levels after the search for all the fine nodes within the coarse cell are performed.
The Multigrid Method Chapter | 5
73
FIGURE 5.6A Fine node within the coarse boundary edge.
The search for nodes within the boundary faces (edges for 2D) will continue after the search for nodes within the cells are performed if the user specifies that the physical wall of the geometry is curved. By specifying whether the wall is flat or curved, the program can save time by skipping the computation of interconnectivity relationship at the boundary of the geometry. Basically, the difference between this algorithm and the algorithm depicted in Fig. 5.5 is that the normal vector, n, is pointing outward at the boundary face instead of pointing inward. This algorithm needs to fulfill two criteria in order for a node to be within the boundary face. The first criterion to be fulfilled is depicted in Fig. 5.6A for 2D, where the dot product between the outward normal vector, n, and the vector, p, pointing from the center of the boundary face must be greater than or equal to zero (i.e.,p.n $ 0). The following equation states first criterion for 3D: 8 9 < nx = ðpn Þ 5 fpg : fng 5 px py pz : ny 5 ðpx nx 1 py ny 1 pz nz Þ $ 0 ð5:8Þ : ; nz To ensure that vector p is not pointing to a node that is very far away from the boundary faces, the boundary faces will be enclosed within a bounding box, which has a dimension of adding both the maximum and minimum x-, y-, and z-coordinates of the three nodes and the maximum length of the three edges that form the boundary face. The second criterion to be fulfilled is depicted in Fig. 5.6B for 2D where the dot product of both p1 .v and p2 .v must be different in sign (i.e., (pv)1 $ 0 and (pv)2 # 0 or vice versa). For 3D implementation, it is more complicating than 2D, whereby it needs to compare the sign for the various vectors as stated in Eqs. (5.9a,b,c) and (5.10a,b,c). To fulfill this criterion, Eqs. (5.9a) to (5.9c) or (5.10a) to (5.10c) must be satisfied together.
74
Computational Fluid-Structure Interaction FIGURE 5.6B Different sign computed between p1 U v and p2 U v.
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 : ; vz
.
0
,
0
ð5:9aÞ ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 > vx > < = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 > : > ; vz 8 9 > vx > < = p3z : vy 5 ðp3x vx 1p3y vy 1p3z vz Þ3 > : > ; vz
.
0
,
0
ð5:9bÞ ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 < vx = 2 p2z : vy 5 p2x vx 1p2y vy 1p2z vz : ; 8 vz 9 < vx = 3 p3z : vy 5 p3x vx 1p3y vy 1p3z vz : ; vz
.
0
,
0
ð5:9cÞ
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 : ; vz
,
0
.
0
ð5:10aÞ
The Multigrid Method Chapter | 5
ðpv Þ1 5 fpg1 : fvg 5 p1x
p1y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
ðpv Þ2 5 fpg2 : fvg 5 p2x
p2y
ðpv Þ3 5 fpg3 : fvg 5 p3x
p3y
8 9 < vx = p1z : vy 5 ðp1x vx 1p1y vy 1p1z vz Þ1 : ; 8 vz 9 < vx = p3z : vy 5 ðp3x vx 1p3y vy 1p3z vz Þ3 : ; vz
75
,
0
.
0
ð5:10bÞ 8 9 < vx = p2z : vy 5 ðp2x vx 1p2y vy 1p2z vz Þ2 , 0 : ; 8 vz 9 < vx = p3z : vy 5 ðp3x vx 1 p3y vy 1 p3z vz Þ 3 . 0 : ; vz ð5:10cÞ
The purpose of this second criterion is mainly to ensure that the nodes being tested is within the boundary face instead of testing other node belonging to the neighboring faces, as shown in Fig. 5.6C for 2D. This figure shows that the sign for these two dot products are the same, which both the values of dot product are less than zero (i.e., p1 .v , 0 and p2 .v , 0, having the same sign). If the node being tested fulfills these two criteria, then the node is projected downward onto the boundary edge (2D) and the projected lengths between this projected node and the two edge nodes are computed for the transfer operator alogorithms, as depicted in Figure 5.7. This is to ensure that the residual transfer is conservative. The above-described algorithm can be summarized in the flowchart as shown in Fig. 5.8 and the pseudo code for the algorithm is given after the flowchart.
FIGURE 5.6C Node belonging to the neigboring edge (i.e., ~ p 1 U~ v , 0 and ~ p 2 U~ v , 0).
76
Computational Fluid-Structure Interaction FIGURE 5.7 Projected lengths for the transfer operator algorithms (2D).
! ---------------------------------------------------------------------------------------------------! Pseudo-code for finding nodes enclose within cell using dot product of two vectors ! ---------------------------------------------------------------------------------------------------subroutine containabc variables definition initialize all the arrays do loop for number of multigrid levels (loop 1) initialize flag as 0 for all the nodes do loop for number of zones (loop 2) loop for number of cells within the zone (loop 3) establish inward normal vectors ~ n for the 4 faces of the cell enclose cell within an imaginary bounding box do loop for number of nodes within the zone (loop 4) if (flag node.eq.0 .and. within the bounding box) then establish vectors ~ p for the 4 faces of the cell compute ~ p .~ n for the 4 faces if (~ p .~ n $ 0 for the 4 faces) then flag node 5 1 call subroutine to calculate volumes for the 4 sub-tetrahedral & volume of the cell go back to loop 4 go back to loop 3 go back to loop 2 if (geometry got curvature) call subroutine to test for nodes within boundary faces go back to loop 1 return end
77
The Multigrid Method Chapter | 5
5.5 MESH-TO-MESH TRANSFER OPERATORS Following the two-dimensional approach presented in Ref. [7] for data transfer within the domain and incorporate the new technique developed here for a curved boundary, there are two classes of mesh-to-mesh transfer operators being implemented in the method for three-dimensions, which are restriction operators and prolongation operators. List of nodes and cells within the specified zones for three-grid levels
Fine-to-coarse grids (1--2 grid level – Group 1) Coarse-to-coarsest grids (2--3 grid level – Group 2)
Initialize the flag as 0 for all the nodes
List of cells in the array for the zone Go to the next group of grid levels (eg., 2--3 grid level for three-grid levels)
Establish the inward normal vectors, n, for the 4 faces of the cells
Enclose the cell within an imaginary bounding box
Is the node flag 0 and within the bounding box?
No
Yes Establish the vectors, p, for the 4 faces of the cell
Compute p .n for the 4 faces
Is p .n ≥ 0 for the 4
No
Yes
Go to the next zone
Go to the next in the list for the zone
Go to the next cell in the list for the zone
List of nodes in the array for the zone
(1) Compute the volumes for the 4 subtetrahedral within the cell and the volume of the cell. (2) Flag the node as 1
A
B
C
D
E
FIGURE 5.8 Flowchart for interconnectivity relationships between meshes algorithm for threegrid levels.
78
Computational Fluid-Structure Interaction
A
B
C
D
End of the list of nodes for the zone?
E
No
Yes No
End of the list of cells for the zone?
Yes End of zones specified by user?
No
Yes No Curvature in the geometry?
Yes Call subroutine to test for nodes that contain within the boundary faces using algorithm depicted in Fig. 5.6 Yes No
End of three-grid levels?
Yes End of algorithm
FIGURE 5.8 (Continued).
5.5.1 Flow-Field Variables Transfer Operators Once the interconnectivity relationship between meshes has been determined, the restriction transfer operator, Thh11 , that transfers the flow-field values from the fine grid to the coarse grid will be as follows: W1 5
Va Wa 1 Vb Wb 1 Vc Wc 1 Vd Wd Va 1 Vb 1 Vc 1 V d
ð5:11Þ
Fig. 5.9A depicts the operation of this transfer operator for those nodes within a cell. Lower case letters denote fine-grid nodes and arabic numbers
The Multigrid Method Chapter | 5
79
FIGURE 5.9A Transfer of flow-field values from the fine mesh to the coarse mesh using restriction transfer operator.
FIGURE 5.9B Transfer of variables from fine nodes to the coarse node at the boundary.
denote coarse-grid nodes for all the figures in this section. The flow-field values at the coarse node 1, which is contained in the fine cell formed by nodes a, b, c, and d, is a weighted average of the values at those nodes. Va is the volume of the sub tetrahedron cell with vertices b, 1, c, and d. Vb is the volume of the sub tetrahedron cell with vertices a, 1, c, and d. Vc is the sub tetrahedron cell with vertices a, 1, d, and b. Vd is the subtetrahedron cell with vertices a, 1, b, and c. The volume of the corresponding tetrahedron cell is the one opposite to the node. According to Eq. (5.11), if node 1 coincides with node a, then W1 will be equal to Wa. At the boundary with curvature, the flow-field values are transferred from the fine nodes of the triangle formed by vertices a, b, and c onto the coarse node 1 using an area-weighted contribution calculation, as illustrated in Fig. 5.9B. The restriction transfer operator, Thh11 , that transfers the flow field values is as follows: W1 5
Aa W a 1 A b W b 1 Ac W c Aa 1 Ab 1 Ac
ð5:12Þ
Fig. 5.9B illustrates the operation of this transfer operator for those nodes within the boundary triangles. The flow-field values at the coarse node 1,
80
Computational Fluid-Structure Interaction FIGURE 5.10 An illustration of transfer of residuals using prolongation operator will not be conserved.
which is projected downward onto the fine triangle face formed by nodes a, b, and c, is a weighted average of the values at those nodes. Aa is the area of the triangle with vertices b, 1, and c. Ab is the area of the triangle with vertices c, 1, and a. Ac is the area of the triangle with vertices a, 1, and b. The area of the corresponding triangle is the one opposite to the node. Similarly, according to Eq. (5.12), if node 1 coincides with node a, then W1 will be equal to Wa.
5.5.2 Residual Transfer Operators In the current MG method, the restriction transfer operator is used for transferring residuals so that the transferred residuals will be conserved. Fig. 5.10 illustrates the reason why using prolongation operators to transfer residuals will not be conserved. For illustration purposes, Fig. 5.10 is illustrated with a 2D mesh. The coarse-grid nodes 1, 2, and 3 are super-positioned onto the three fine cells as shown in the figure. The first fine cell is formed by node a, b, and f, the second fine cell is formed by node i, j, and n, and the third fine cell is formed by node l, m, and q. If the residuals are transferred using prolongation operators, then the residuals at these nine fine nodes will be transferred to the three coarse nodes, while the residual at fine node k will not be transferred to any of the coarse nodes. It is possible to encounter a situation whereby the residuals of these nine fine nodes are zero and the residual of fine node k is nonzero. A negative transfer will result in the three coarse nodes to have zero residual if using prolongation operators. To ensure that the overall convergence will be conserved and to avoid negative transferring between meshes, the selection of the class of operators to be used in transferring the residuals is very important to ensure that the sum of coarsegrid residuals after the transfer are always conserved, that is equal to the total of the fine-grid residuals. where the transfer of residuals from the fine grid to the coarse grid is based on a weighted average of the residual from the fine node to the four coarse nodes that form the coarse tetrahedral cell. As depicted in Fig. 5.11A,
The Multigrid Method Chapter | 5
81
FIGURE 5.11A Transfer of residuals from the fine mesh to the coarse mesh using restriction transfer operator.
the residual from node a is transferred to coarse nodes 1, 2, 3, and 4 by the transfer operator, Q h11 h . The transferred residual to a particular coarse node is the summation of all the transferred weighted-average residuals from other fine nodes in those tetrahedra surrounding this coarse node, as shown in the following equations, where the contribution from fine-grid node a is calculated and added to the residual of coarse-grid node 1: ~~ 5 R ~~ 1 R 1 1 ~~ 5 R ~~ 1 R 3 3
V1 R~~ a ; V1 1 V2 1 V3 1 V4
V3 R~~ a ; V 1 1 V 2 1 V 3 1 V4
~~ 1 R~~ 2 5 R 2 R~~ 4 5 R~~ 4 1
~~ V2 R a ; V1 1 V2 1 V 3 1 V 4
~~ V4 R a V1 1 V2 1 V3 1 V4
ð5:13Þ
where V1 is the volume of the subtetrahedron cell with vertices 2, a, 3, and 4. V2 is the volume of the subtetrahedron cell with vertices 1, a, 3, and 4. V3 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 4. V4 is the volume of the sub tetrahedron cell with vertices 1, a, 2, and 3. The volume of the corresponding tetrahedron cell is the one opposite to node. It is easy to show that this transfer is conservative in the sense that the total fine mesh residuals are equal to the coarse mesh ones. The residual for a find-grid node on the boundary face of the curvature geometry is transferred to the coarse nodes of the triangle formed by vertices 1, 2, and 3 using an area-weighted contribution calculation as shown in Fig. 5.11B. The transfer operator, Q h11 h , that transfers the residuals from the fine grid to the three coarse-grid nodes is as follows: ~~ 5 R ~~ 1 R 1 1
A1 R~~ a A 1 1 A 2 1 A3
82
Computational Fluid-Structure Interaction FIGURE 5.11B Transfer of residuals from fine node to the coarse nodes at the boundary.
R~~ 2 5 R~~ 2 1
~~ A2 R a ; A1 1 A2 1 A3
~~ 5 R ~~ 1 R 3 3
A3 R~~ a A1 1 A2 1 A3
ð5:14Þ
Fig. 5.11B depicts the operation of this transfer operator typically for those nodes on the boundary faces. The residual at the fine node a, which is projected downward onto the coarse triangle formed by nodes 1, 2, and 3, is a weighted average of the value at this fine node. The transferred residual is the summation of all the transferred residuals from other fine nodes to these coarse nodes. This also ensures that the residual transfer is conservative on the boundary. It has been proven that the total fine-grid residuals are equal to the coarse grid ones and this is achieved by performing a summation of the residuals on both fine and coarse grids and the totals are equal to each other [7]. This proves that this restriction transfer operator is conservative.
5.5.3 Correction Transfer Operators Prolongation operators are used to transfer corrections of the flow-field variables from the coarse mesh to the fine mesh, as depicted in Figs. 5.12A and 5.12B. As shown in Eq. (5.15), the correction, dWh11 , is the difference 1 between the newly computed value on the coarse grid, Wh11 , and the initial ð0Þ value that was transferred from the fine grid, Wh11 . ð0Þ 1 dWh11 5 Wh11 2 Wh11
ð5:15Þ
Corrections are transferred to the fine mesh by the prolongation operator, I
h h11 : h dWh11 vh 5 Ih11
and
Wh1 5 Wh 1 vh
According to Fig. 5.12A, the correction of the flow-field variables transferred from the coarse nodes 1, 2, 3, and 4 to the fine node a is a weighted
The Multigrid Method Chapter | 5
83
FIGURE 5.12A Transfer of corrections from the coarse mesh to the fine mesh using prolongation transfer operator.
FIGURE 5.12B Transfer of corrections from the coarse nodes to the fine node at the boundary.
average of the corrections at these nodes and the expression for the transferred correction is ðva Þh 5
V1 ðdW1 Þh11 1 V2 ðdW2 Þh11 1 V3 ðdW3 Þh11 1 V4 ðdW4 Þh11 V1 1 V2 1 V3 1 V4
ð5:16Þ
where V1 is the volume of the subtetrahedron cell with vertices 2, a, 3, and 5. V2 is the volume of the subtetrahedron cell with vertices 1, a, 3, and 5. V3 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 5. V4 is the volume of the subtetrahedron cell with vertices 1, a, 2, and 3. The volume of the corresponding tetrahedron cell is the one opposite to the node. The corrections of the flow-field variables at the boundary are transferred from the coarse nodes of the triangle formed by vertices 1, 2, and 3 onto the fine node a. The transfer operator, I hh11 , that transfers the corrections is ðva Þh 5
A1 ðdW1 Þh11 1 A2 ðdW2 Þh11 1 A3 ðdW3 Þh11 A 1 1 A2 1 A3
ð5:17Þ
84
Computational Fluid-Structure Interaction
Fig. 5.12B shows the operation of this transfer operator for those nodes on the boundary faces. The correction at the fine node a, which is projected downward onto the coarse face formed by nodes 1, 2, and 3, is a weighted average of the values at those nodes. A1 is the area of the triangle with vertices 2, a, and 3. A2 is the area of the triangle with vertices 1, a, and 3. A3 is the area of the triangle with vertices 1, a, and 2. The area of the corresponding triangle is the one opposite to the node. Similarly, according to Eq. (5.17), if node a coincides with node 1, then ðva Þh will be equal to ðdW1 Þh11 .
5.6 THE MULTIGRID METHOD FOR STRUCTURAL DYNAMIC SOLVER In this section, an efficient unstructured MG scheme developed for the equation of dynamic equilibrium is introduced. The MG algorithm is described as follows. The discretized dynamic equation can be expressed by equation (4.25), which is repeated here for convenience: ΔU ~~ n11;m 5R Δτ
ð5:18Þ
where U is the velocity vector used in CSD solver. This equation is solved iteratively by a dual-time-stepping scheme. An pseudo code is given below to illustrate the basic procedures of the MG scheme. The outer cycles are based on the physical time t, which is numbered from 1 to ktmax, whereas the inner cycles are based on the pseudo time τ, which is numbered from 1 to max_itr_sub. The grid levels range from 1 to nmg, where 1 is the finest level and nmg is the coarsest level. Phh11 is the prolongation operator from level h 1 1 to h, and Qh11 and Thh11 are the residual transfer and restriction h operators from level h to h 1 1. ALGORITHM START call setupinterconnect !---. build up inter-connectivity relationship between levels DO kt 5 1, ktmax !---. start of physical time step DO mg 5 1, nmg call cvvol(mg) !---. calculate the volumes for control volume (CV) and cell call clhaut(mg) !---. calculate CV characteristic length for determination !---. of Δτ and Δt computation ENDDO call dtsize(mg 5 1) !---. computation of the time-step size for the finest level DO itersub 5 1, max_itr_sub !---. start of sub-iteration
The Multigrid Method Chapter | 5
85
DO ialpha 5 1, 5 !---. 5-stage Runge-Kutta time integration process kt call fvsolver(mg 5 1) !---. solve equation (3.1) for ~ U1 following the way !--- . described in forgoing sections ENDDO DO mg 5 2, nmg call transfrsol(mg-1) !---. h to h 1 1 ð0Þ U h ). !---. (~ U h11 5 Thh11 ~ call transfresd(mg-1) !---. h to h 1 1 ~ ~ ð0Þ R~ U !---. (R~ h11 5 Q h11 h h ).
restrict
solution
from
restrict
residual
from
IF (itersub.eq.1) THEN call dtsize(mg) !---. calculate the time-step size for the current level ENDIF DO ialpha 5 1, 5 kt call fvsolver(mg) !---. solve equation (3.1) for ~ U mg based on the initial ð0Þ ~ ð0Þ !---. solution ~ U h11 and initial residual R~ h11 ENDDO ENDDO DO mg 5 nmg, 2, -1 call transfrcorct(mg) !---. prolongate correction from h 1 1 to h 1 1 ð0Þ 1 !---. (~ Uh 5 ~ U h 1 Ihh11 ~ U h11 ). Here ~ U h11 2 ~ U h is !---. the updated solution for the finer grids ENDDO ENDDO !---. end of sub-iteration DO mg 5 2, nmg call transfrsol(mg-1) !---. restrict solution from h to h 1 1 ENDDO DO mg 5 1, nmg call updgrid(mg) !---. finally update the solid mesh ENDDO ENDDO !---. end of physical time step ALGORITHM END
One should note that we always use the initial intermesh connectivity, which is built up before the first time-step starts, to perform the mesh-tomesh transfer operations. In order to enable the solver to tackle more geometrically complex structures, the search for nodes within the boundary faces (edges for 2D) will continue after the search for nodes within the cells are performed if the user specifies that the physical wall of the geometry is
86
Computational Fluid-Structure Interaction Coarser-grid node
Coarser-grid boundary
Finer-grid node
Finer-grid boundary
Projection node for P
FIGURE 5.13 Schematic of boundary node projection algorithm, Ω is the computational domain.
Node P p1 1
θ1 θ 2
p2 2
V Ω
curved (e.g., see Fig. 5.13 ), a normal procedure cannot find a corresponding finer grid cell for the coarser grid node P. It means that the node P cannot get a proper solution vector from the finer grid, which can pose a serious problem while solving the dynamic equation (4.28). Following the boundary node projection algorithm proposed in Section 5.6, we resort to finding a projection node for node P in the nearest finer mesh cell face and then use this projected node to perform the necessary MG operations.
REFERENCES [1] R.P. Fedorenko, The speed of convergence of one iterative process, U.S.S.R., Comput. Math. Math. Phys. 4 (3) (1964) 227235. [2] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput. 31 (1977) 333390. [3] A. Jameson, Solution of the Euler equations for two dimensional transonic flow by a multigrid method, Appl. Math.Comput 13 (1983) 327355. [4] Jameson A., Time Dependent Calculations Using Multigrid, with Applications to Unsteady flows Past Airfoils and Wings, AIAA Paper 91-1596, June 1991. [5] E. Turkel, Preconditioned Methods for Solving the Incompressible and Low Speed Compressible Equations, J. Comput. Physics 72 (1987). [6] Mavriplis D.J., Large-scale Parallel Viscous Flow Computations Using An Unstructured Multigrid Algorithm, NASA/CR-1999-209724, ICASE Report No. 99-44, Nov 1999. [7] C.H. Tai, Y. Zhao, A Finite Volume Unstructured Multigrid Method for Efficient Computation of Unsteady Incompressible Viscous Flows, Int. J. Num. Meth. Fluids 46 (1) (2004) 5984.