@prefix vivo: .
@prefix edm: .
@prefix ns0: .
@prefix dcterms: .
@prefix skos: .
vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Mechanical Engineering, Department of"@en ;
edm:dataProvider "DSpace"@en ;
ns0:degreeCampus "UBCV"@en ;
dcterms:creator "Malik, Shahzaib"@en ;
dcterms:issued "2019-05-14T14:37:57Z"@en, "2019"@en ;
vivo:relatedDegree "Master of Applied Science - MASc"@en ;
ns0:degreeGrantor "University of British Columbia"@en ;
dcterms:description """In this thesis, we present a new procedure for mesh adaptation for wakes. The approach
starts by tracking the wake centerline with an initial isotropic unstructured mesh. A vertex-centered finite volume method is used, and the velocity field is obtained from solution reconstruction. The velocity data is integrated numerically using an adaptive fourth-order
Runge-Kutta method. We insert the wake centerline into the existing unstructured mesh
as an internal boundary and use a metric-based anisotropic mesh adaptation to generate
anisotropic cells in regions with large second derivatives of flow variables. In the second
step, the problem is solved on adapted mesh and a new wake centerline is tracked. We
then move the previous wake centerline (which is now a part of adapted mesh) to match
the centerline obtained from the adapted mesh. To move the wake centerline, a solid mechanics analogy is used and the linear elasticity equation is solved on the adapted mesh.
As a result, the displacement is propagated throughout the mesh and the already adapted
regions along the wake centerline are preserved. The process is then followed for subsequent
cycles of anisotropic mesh adaptation to obtain a more accurate approximation of the wake
centerline. As an alternate strategy for obtaining an anisotropic mesh in the wake, we take the first geometry, together with the captured wake centerline from an unstructured triangular mesh,
as an initial geometry to produce a quad dominant mesh, using an advancing layer method. The correctness of the streamline tracking algorithm is verified using an analytical velocity field. The mesh morphing approach is tested using the method of manufactured solutions,
demonstrating that the linear finite element solution is second-order accurate. The results
of laminar flow test cases for the attached and separated flow are presented and compared with some well-established numerical results in the literature. Our results show that the advancing layer mesh is more efficient in resolving the wake. In the end, one case for turbulent
subsonic flow is considered. For turbulent flow, a cell-centered finite volume method is used
and we only track the wake centerline at different angles of attack."""@en ;
edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/70155?expand=metadata"@en ;
skos:note "Mesh Adaptation for Wakes viaSurface InsertionbyShahzaib MalikBEng., Queen Mary University of London (UK), 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)May 2019© Shahzaib Malik, 2019The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:Mesh Adaptation for Wakes via Surface Insertionsubmitted by Shahzaib Malik in partial fulfillment of the requirements forthe degree of Master of Applied Sciencein Mechanical EngineeringExamining Committee:Dr. Carl Ollivier-GoochSupervisorDr. Dana GrecovSupervisory Committee MemberDr. Steven RogakSupervisory Committee MemberiiAbstractIn this thesis, we present a new procedure for mesh adaptation for wakes. The approachstarts by tracking the wake centerline with an initial isotropic unstructured mesh. A vertex-centered finite volume method is used, and the velocity field is obtained from solution re-construction. The velocity data is integrated numerically using an adaptive fourth-orderRunge-Kutta method. We insert the wake centerline into the existing unstructured meshas an internal boundary and use a metric-based anisotropic mesh adaptation to generateanisotropic cells in regions with large second derivatives of flow variables. In the secondstep, the problem is solved on adapted mesh and a new wake centerline is tracked. Wethen move the previous wake centerline (which is now a part of adapted mesh) to matchthe centerline obtained from the adapted mesh. To move the wake centerline, a solid me-chanics analogy is used and the linear elasticity equation is solved on the adapted mesh.As a result, the displacement is propagated throughout the mesh and the already adaptedregions along the wake centerline are preserved. The process is then followed for subsequentcycles of anisotropic mesh adaptation to obtain a more accurate approximation of the wakecenterline.As an alternate strategy for obtaining an anisotropic mesh in the wake, we take the firstgeometry, together with the captured wake centerline from an unstructured triangular mesh,as an initial geometry to produce a quad dominant mesh, using an advancing layer method.The correctness of the streamline tracking algorithm is verified using an analytical veloc-ity field. The mesh morphing approach is tested using the method of manufactured solutions,demonstrating that the linear finite element solution is second-order accurate. The resultsof laminar flow test cases for the attached and separated flow are presented and comparediiiAbstractwith some well-established numerical results in the literature. Our results show that the ad-vancing layer mesh is more efficient in resolving the wake. In the end, one case for turbulentsubsonic flow is considered. For turbulent flow, a cell-centered finite volume method is usedand we only track the wake centerline at different angles of attack.ivLay SummaryWith the advancements in computing resources, it has become important to develop CFDmethods in such a way that they are efficient in terms of computational cost while main-taining and improving the solution accuracy at the same time. The errors in numericalsimulations can influence the flow physics and solution accuracy on different scales and arehighly dependent on mesh spacing. Mesh refinement can help to reduce these errors but uni-form mesh refinement is not computationally efficient. Many researchers focus on selectivelyrefining the mesh in areas where the impact of errors on the solution is greater and such aprocess is known as solution adaptive refinement. This research introduces a new approachto anisotropic mesh adaptation for wakes. The results from different meshes are comparedto identify the more efficient mesh in tracking the wake.vPrefaceThe research work presented in this thesis is the outcome of a close working relationshipbetween Dr. Carl Ollivier Gooch and Shahzaib Malik. The algorithm development, conceptformulation, data analysis, and manuscript preparation were done by the author, ShahzaibMalik, under valuable guidance from Dr. Carl Ollivier Gooch in the department of Mechan-ical Engineering at the University of British Columbia.Some of the results in Chapter 4 of this thesis were presented by the author at the AIAASciTech 2019 Forum in January 2019. Parts of Chapter 1, Chapter 3 and Chapter 4 werepublished in the paper [31]: “Shahzaib Malik and Carl F Ollivier Gooch. Mesh adaptationfor wakes via surface insertion. In AIAA SciTech 2019 Forum, AIAA Paper 2019-1996,2019”. https://arc.aiaa.org/doi/abs/10.2514/6.2019-1996viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Finite Volume Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Mesh Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Structured and Unstructured Meshes . . . . . . . . . . . . . . . . . . 8viiTable of Contents2.2.2 Isotropic and Anisotropic Meshes . . . . . . . . . . . . . . . . . . . . 102.3 Solution Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.4 Metric-Based Anisotropic Mesh Adaptation . . . . . . . . . . . . . . . . . . . 152.4.1 Aspect ratio and orientation from metric . . . . . . . . . . . . . . . . 182.5 Operations for Mesh Quality Improvement . . . . . . . . . . . . . . . . . . . 192.5.1 Face Swapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5.2 Vertex Insertion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.5.3 Vertex Removal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.5.4 Vertex Movement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.6 Streamlines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.7 Runge-Kutta Method and Adaptive Stepsize Control . . . . . . . . . . . . . 233 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.1 Wake Centerline Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2 Wake Centerline Insertion into Unstructured Mesh . . . . . . . . . . . . . . . 283.3 Anisotropic Mesh Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.4 Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.5 Algorithm Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.6 Advancing Layer Mesh Generation . . . . . . . . . . . . . . . . . . . . . . . . 334 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374.1 Verification Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.1.1 Streamline Tracking Algorithm . . . . . . . . . . . . . . . . . . . . . . 384.1.2 Mesh Morphing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 404.2 Laminar Flow Test Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.2.1 Viscous flow at Re = 2100, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . 424.2.2 Viscous flow at Re = 2100, Ma = 0.5 and α = 2◦ . . . . . . . . . . . . 504.2.3 Viscous flow at Re = 5000, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . 574.3 Turbulent Subsonic Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65viiiTable of Contents5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.3 Recommended Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73AppendicesA Command-line Options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A.1 Wake Centerline Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80A.2 Wake Centerline Insertion into the Mesh . . . . . . . . . . . . . . . . . . . . 81A.3 Obtaining the metric file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81A.4 Anisotropic Mesh Adaptation (Metric-Based) . . . . . . . . . . . . . . . . . . 81A.5 Mesh Morphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82A.6 Advancing Layer Mesh Generation . . . . . . . . . . . . . . . . . . . . . . . . 82ixList of Tables4.1 Error in solution (values at t = 2.20) . . . . . . . . . . . . . . . . . . . . . . . 394.2 Comparison of lift and drag coefficients for NACA 0012 airfoil at Re = 5000,Ma = 0.5 and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60xList of Figures1.1 Boeing 717 (engines installed at the rear) [47] . . . . . . . . . . . . . . . . . 42.1 Structured mesh around leading edge of a NACA 0012 airfoil . . . . . . . . . 82.2 Unstructured mesh around leading edge of NACA 0012 airfoil . . . . . . . . 92.3 Isotropic and anisotropic mesh fragments . . . . . . . . . . . . . . . . . . . . 102.4 Reconstruction stencils [39] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.5 Deformation of tensor Lp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.6 Anisotropy components for a certain point . . . . . . . . . . . . . . . . . . . . 172.7 Face swapping in metric space . . . . . . . . . . . . . . . . . . . . . . . . . . 192.8 Vertex insertion for an interior edge . . . . . . . . . . . . . . . . . . . . . . . 202.9 Vertex insertion for a boundary edge . . . . . . . . . . . . . . . . . . . . . . . 212.10 Vertex removal for an interior edge . . . . . . . . . . . . . . . . . . . . . . . . 212.11 Vertex removal for a boundary edge . . . . . . . . . . . . . . . . . . . . . . . 222.12 Streamlines around NACA 0012 airfoil at Re = 5000, Ma = 0.5 and α = 0◦ . 232.13 Adaptive stepsize control for RK4 [44] . . . . . . . . . . . . . . . . . . . . . . 253.1 Separated flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.2 Algorithm overview (mesh adaptation loop) . . . . . . . . . . . . . . . . . . . 323.3 Initial geometry to generate advancing layer mesh . . . . . . . . . . . . . . . 333.4 Development of wake behind an airfoil [59] . . . . . . . . . . . . . . . . . . . 333.5 Number of layers across thickness δ (x) . . . . . . . . . . . . . . . . . . . . . 36xiList of Figures4.1 Isotropic unstructured mesh around a NACA 0012 airfoil with 7310 cells and3871 vertices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.2 Streamline tracking for an analytical velocity field . . . . . . . . . . . . . . . 394.3 Error in solution at each point along the curve . . . . . . . . . . . . . . . . . 404.4 L2 norm of error for linear elasticity . . . . . . . . . . . . . . . . . . . . . . . 414.5 Mesh morphing to displace horizontal line from 0◦ to 3◦ . . . . . . . . . . . . 424.6 Streamline tracking at Re = 2100, Ma = 0.5 and α = 0◦ . . . . . . . . . . . 434.7 First adapted mesh (black) and morphed mesh (red) before second adaptation 434.8 Adapted meshes and flow solutions with wake centerline at Re = 2100, Ma =0.5 and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.9 Wake centerline position - initial (black) and final (red) at Re = 2100, Ma =0.5 and α = 0◦- Closeup view at the trailing edge . . . . . . . . . . . . . . . 454.10 Initial geometry for advancing layer mesh at Re = 2100, Ma = 0.5 and α = 0◦ 454.11 Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 2100, Ma = 0.5 and α = 0◦ , showing better solution behavior in thewake for the advancing layer mesh . . . . . . . . . . . . . . . . . . . . . . . . 454.12 Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 2100, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . 474.13 Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 2100, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . . . . . . 484.14 Velocity profiles in the wake for triangular and advancing layer quad dominantmeshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.15 Streamline tracking algorithm result at Re = 2100, Ma = 0.5 and α = 2◦ . . 504.16 First adapted mesh (black) and morphed mesh (red) before second adaptation 504.17 Adapted meshes and flow solutions with wake centerline at Re = 2100, Ma =0.5 and α = 2◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.18 Wake centerline position - initial (black) and final (red) at Re = 2100, Ma =0.5 and α = 2◦- Closeup view at the trailing edge . . . . . . . . . . . . . . . . 52xiiList of Figures4.19 Initial geometry for quad dominant mesh at Re = 2100, Ma = 0.5 and α = 2◦ 524.20 Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 2100, Ma = 0.5 and α = 2◦ , showing better solution behavior in thewake for the advancing layer mesh . . . . . . . . . . . . . . . . . . . . . . . . 524.21 Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 2100, Ma = 0.5 and α = 2◦ . . . . . . . . . . . . . . . . . . . . . . . 544.22 Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 2100, Ma = 0.5 and α = 2◦ . . . . . . . . . . . . . . . . . 554.23 Velocity profiles in the wake for triangular and advancing layer quad dominantmeshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.24 Wake centerline for separated flow . . . . . . . . . . . . . . . . . . . . . . . . 584.25 Wake centerline position - initial (black) and final (red) at Re = 5000, Ma =0.5 and α = 0◦- Closeup view at the trailing edge . . . . . . . . . . . . . . . 584.26 First adapted mesh (black) and morphed mesh (red) before second adaptation 584.27 Adapted meshes and flow solutions with wake centerline at Re = 5000, Ma =0.5 and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.28 Initial geometry for quad dominant mesh at Re = 5000, Ma = 0.5 and α = 0◦ 604.29 Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 5000, Ma = 0.5 and α = 0◦ , showing better solution behavior in thewake for the advancing layer mesh . . . . . . . . . . . . . . . . . . . . . . . . 614.30 Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 5000, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . . . . . . . . . . . . 624.31 Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 5000, Ma = 0.5 and α = 0◦ . . . . . . . . . . . . . . . . . 634.32 Velocity profiles in the wake for triangular and advancing layer quad dominantmeshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.33 Mapping from reference triangle to cubic triangular cell . . . . . . . . . . . . 664.34 Unstructured hybrid mesh (n.DoF = 25, 088) . . . . . . . . . . . . . . . . . . 67xiiiList of Figures4.35 Wake centerlines at Re = 6× 106, Ma∞ = 0.15 and α = 0°, 4°, 10° using thestreamline-tracking algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 67xivNomenclatureRoman SymbolsAR aspect ratioc chord lengthCD drag coefficientCL lift coefficientd distance function, linear elasticity coefficientsdymin off-wall spacingE Young’s modulusF flux dyadH Hessian matrixL deformation tensorlM metric lengthlmaxMmaximum metric lengthlminMminimum metric lengthM metricMa Mach numberxvNomenclatureRe Reynolds numberR rotation matrixS source termU solution vectoru x-velocityv y-velocityGreek Symbolsα angle of attackΛ size matrixν Poisson’s ratioxviAcknowledgementsFirst of all, I would like to thank my supervisor Dr. Carl Ollivier Gooch for his valuabletime, supervision and constant support throughout my research and graduate studies. Tobe honest, Carl’s enthusiasm for research and involvement with research students fueled meduring my time at UBC. I am glad that I worked with him and learned so many things underhis supervision.I would like to express my gratitude to ANSYS Canada, the Natural Sciences and En-gineering Research Council of Canada and the University of British Columbia for providingfinancial support.I am thankful to the members of the supervisory committee: Dr. Steven Rogak and Dr.Dana Grecov, for generously offering their time to review my thesis.Also, I am very thankful to all my colleagues at Advanced Numerical Simulation Lab:Gary, Reza, Jasmeet, Ehsan and Akhil for their encouragement and valuable discussionsregarding the challenges faced during the research work.Finally, thanks to my friends: Faisal, Zhihao, Somesh, Harsh, Claire, Juuso, Alondra,Miguel, Adriana, Ratul and Sarah for being with me.xviiDedicationI would like to thank my family especially my mom, dad and sister as their support, trustand encouragement has always made be strong enough to overcome all the challenges.To my family.xviiiChapter 1IntroductionComputational fluid dynamics (CFD) has become an increasingly valuable tool in the lastfew decades with many industrial and non-industrial applications. These applications includeaerodynamics of motor vehicles and aircraft, oil recovery, turbo-machinery, hydrodynamicsof ships and astrophysics [9, 11, 57]. Therefore, it is immensely important to further developCFD methods in such a way that they are efficient in terms of computational cost whilemaintaining and improving the solution accuracy at the same time. One of the main objec-tives of CFD analysis is to study the physical processes that occur around different objects.These processes can be the result of phenomena such as shock waves, boundary layers, flowseparation and turbulence. The processes are governed by different equations depending onthe type of problem under investigation. It is possible to obtain an analytical solution tomany conservation equations but only under certain simple conditions. Analytical solutionsare usually possible when the equation is linear and the geometry and boundary conditionsare not very complicated. In reality, most of the problems in fluid dynamics are governedby the Navier-Stokes equations as the effect of viscosity cannot be ignored very close to theboundary when drag or heat transfer is important. The Navier-Stokes equations are non-linear and very complicated to solve, thus making analytical solutions impossible especiallyfor complex geometries and this is where numerical methods become very useful. Differentnumerical methods exist that can be used to obtain numerical solutions to the ordinary andpartial differential equations related to the problem [10, 28].When numerical methods are applied to differential equations, the original equations arediscretized and the resulting algebraic equations are then solved through numerical methods.As a result of this discretization, the numerical solution obtained is an approximate solution11.1. Motivation and Objectivesrather than the exact solution. The difference between the exact and numerical solution isthe numerical error. These errors can effects flow physics and solution on different scalesand are highly dependent on the mesh spacing. One way to reduce these errors is to reducethe mesh spacing uniformly and ideally refine the mesh until the solution is grid converged.However, uniform mesh refinement is not computationally efficient and usually, the focus isto identify those areas where the impact of errors on the solution is greater and then reducethe error in those areas by refining the mesh locally. Refining the mesh in such a way isknown as solution adaptive refinement [17, 43].This thesis deals with anisotropic mesh adaptation which is one step in the relentlesseffort to make more efficient use of computing resources.1.1 Motivation and ObjectivesCFD analysis to obtain a numerical solution involves three essential steps: (1) modeling,(2) mesh generation and (3) solution. Modeling refers to specifying the problem geometry,governing equations and the boundary conditions. In the second step, to solve the governingequations numerically on the specified geometry, the continuous domain is divided into aset of discrete regions known as the mesh. The third step involves discretization of theequations on the mesh and solution of the resulting system of algebraic equations to obtainthe numerical solution [54]. For real aerodynamics applications and industry level problems,the mesh generation process is considered to be the most time consuming out of all processesin terms of human time and can be up to 45 times longer compared to the solver time [2, 33].Therefore, it is very important to explore ways to reduce the mesh generation time for agiven level of solution accuracy.A mesh is a discrete representation of the geometry of a CFD problem and has significanteffects on convergence, computational cost and accuracy of the numerical solution [13, 54].Therefore, it is immensely important to obtain a good quality mesh and take into consid-eration the geometric properties of the mesh such as aspect ratio, smoothness, maximum21.1. Motivation and Objectivesand minimum angle, local mesh size and orientation. Any poor quality cells can adverselyaffect the numerical approximation and thus lead to unreliable results and instability in thesolution. There is generally a trade-off between the accuracy and computational cost in com-putational fluid dynamics. In many aerodynamics applications, accuracy can be improvedfor a given computational cost by using high aspect ratio cells to resolve the flow properties.A typical example is a boundary layer, where there are high velocity gradients perpendic-ular to the wall and small gradients parallel to the wall. Obtaining a reasonable accuracyin the boundary layer using fine isotropic cells leads to additional memory usage and un-necessary computational effort. Therefore, it is more efficient to use anisotropic cells withanisotropy aligned along the desired direction as the anisotropic cells help to improve theinverse relationship between solution accuracy and computational cost. Many researcherssuch as Bottaso [6], Huang [18], Apel [4] and Loseille [29] have provided a detailed view onthe advantages of using anisotropic elements and metric-based anisotropic mesh adaptationschemes. Moreover, Rippa [45] showed that for linear reconstruction, the cells should havelonger edge length in the direction where the second derivatives of solution variables aresmall compared to the direction with large second derivatives [49].Why do we care about resolving the wake?A wake is the region of decelerated flow behind a body immersed in a fluid [59]. It can bethought of as the boundary layer separated from the rear of the body. Quite often in aero-dynamics, we study the flow over individual components such as wing, fuselage or tailplaneof the aircraft but when the flow over the entire aircraft is studied, all these components arein the vicinity of each other and the flow from one component will interfere with flow pastanother component and cause interference effects. There are several applications where weneed to take into account the aerodynamics of wake due to these interference effects. Theprime examples to resolve the wake downstream of an airfoil or a wing is an interactionwith the fuselage and the tailplane. Another important application is the wake and flapinteraction noise which is one of the main sources of noise during landing as mentioned in31.1. Motivation and Objectivesthe experimental study by V. Hutcheson, J. Stead and E. Plassman [19]. Moreover, if theengines are installed on the rear of the aircraft as shown in Figure 1.1, the unsteady wakebehavior of the flow from the wing can potentially effect the engine performance to somedegree. Therefore, it is important to resolve the wake so that an optimal configuration canbe obtained which will reduce drag and noise due to interference effects.Figure 1.1: Boeing 717 (engines installed at the rear) [47]The main objectives of this research are:• To resolve the wake in an efficient way and keep a track of where the wake is movingas the mesh is refined. It is possible to track the wake position by tracking the wakecenterline based on the velocity data (see Section 3.1).• Compare the results of mesh adaptation with and without the wake centerline forunstructured triangular meshes using an existing metric-based anisotropic mesh adap-tation scheme developed by Zuniga Vazquez [55].• Solve the linear elasticity problem to morph the adapted mesh to match the new wakecenterline in subsequent cycles of mesh adaptation.• Use the advancing layer mesh generator in our in-housing meshing code to generatean anisotropic quad dominant advancing layer mesh and compare the results withunstructured triangular meshes with and without the wake centerline.41.2. Thesis Outline1.2 Thesis OutlineThe overall structure of this thesis is as follows:Chapter 2 provides background on different topics needed throughout this research. Section2.1 briefly summarizes the finite volume solver. In Section 2.2, different mesh propertiesare discussed. Section 2.3 discusses the solution reconstruction as the velocity data neededto track the streamline (wake centerline) is obtained from solution reconstruction. Section2.4 talks about the metric-based anisotropic mesh adaptation scheme. In Section 2.5, theprincipal operations used for mesh quality improvement by the mesh adaptation scheme arediscussed. The last two sections of Chapter 2 provide an overview of streamlines and theadaptive fourth-order Runge-Kutta method used to numerically integrate the velocity data.In Chapter 3, the overall methodology of anisotropic mesh adaptation for wakes is dis-cussed. First, it includes wake centerline tracking using an adaptive fourth-order Runge-Kutta method from an initial isotropic unstructured mesh. The second step is to insert thewake centerline as an internal boundary into the initial unstructured mesh. Next, it talksabout anisotropic mesh adaptation and re-computing the solution on the adapted mesh toobtain the new wake centerline. Then, the mesh morphing technique used to morph the oldmesh to match the new wake centerline in subsequent adaptation cycles is discussed. Finally,this chapter concludes by discussing the advancing layer mesh generation scheme used togenerate quad dominant meshes.Chapter 4 of this thesis illustrates the results. The verification cases for the streamlinetracking algorithm and the mesh morphing approach are included in Section 4.1. Section4.2 contains three laminar flow test cases to demonstrate the mesh adaptation algorithm.Finally, Section 4.3 includes the streamline tracking results for the turbulent subsonic flow.Chapter 5 ends the thesis with concluding remarks. It contains a summary of the researchwork, conclusions drawn based on the results and the recommendations for the future work.5Chapter 2BackgroundThis chapter provides details of the background material needed to do this research. In thebeginning, a brief overview of finite volume solver is given. In Section 2.2, different meshproperties are discussed. Section 2.3 discusses the solution reconstruction as the velocity dataneeded to track the streamlines is obtained from solution reconstruction. Sections 2.4 and2.5 include details about the metric-based anisotropic mesh adaptation and the principaloperations used for mesh quality improvement respectively. The last two sections of thechapter include discussion on streamlines and the adaptive fourth-order accurate Runge-Kutta method used to numerically integrate the velocity data.2.1 Finite Volume SolverThe equations governing the behavior of fluid flow are a set of partial differential equations(PDE’s) that describe the physical conservation laws. The CFD discretization of these PDE’sneeds to be stable, convergent, accurate and conservative. In the finite volume method,when the flow equations are solved, a conservative discretization is satisfied automaticallythrough the use of integral conservative laws. This means the finite volume discretizationis conservative by nature and a natural choice for solving fluid related problems [8, 25, 35].The finite volume solver begins with governing equations in conservation form, representedas:∂U∂t+ ∂F∂x+ ∂G∂y+ ∂H∂z= S∂U∂t+∇.−→F = S,(2.1)62.2. Mesh Propertieswith U being the solution vector, −→F the flux dyad and S is the source term.Equation 2.1 is integrated over each sub-domain (known as a control volume); applyingGauss’s theorem, we getˆCVi∂U∂tdV +˛∂(CV )i−→F .−→n dA =ˆCViSdV.Assuming that the mesh is not moving, the above expression can be simplified to givedU idt= − 1V˛∂(CV )i−→F .−→n dA+ Si, (2.2)where U i ≡ 1V´CViUdV is the average value of the solution in the control volume i andSi ≡ 1V´CViSdV is the average source term contribution.Equation 2.2 implies that the finite volume solution for any time varying problem can beadvanced from one time level to the next with following steps:• Compute the flux −→F at the surface of the control volume.• Integrate the normal flux −→F .−→n around the control volume boundary.• Compute and integrate the source term S in the control volume.• Advance the average solution value U in time.2.2 Mesh PropertiesThe partial differential equations governing the behavior of fluid flow are defined over a con-tinuum domain. In CFD, before we solve these PDE’s numerically, we divide the continuousdomain into a discrete domain, called the mesh. Depending on the connectivity of cells andvertices in the mesh, the mesh is classified as structured or unstructured. The type of meshto use for CFD analysis of a certain problem depends on the discretization technique, flowphysics and problem geometry. For example, the finite element and finite volume method72.2. Mesh Propertiescan be used with both structured and unstructured meshes while the finite difference methodcan be used only with structured meshes. Apart from distinguishing the mesh based on theconnectivity, another way of categorizing the mesh is by means of shape, orientation andsize of cells throughout the computational domain. This classifies the mesh into categoriesof isotropic and anisotropic. In an isotropic mesh, the edge lengths of cells are roughly equalin all directions whereas, in an anisotropic mesh, the cells have longer edge lengths in onedirection compared to the other. Each of these mesh types are discussed briefly.2.2.1 Structured and Unstructured MeshesIn a structured mesh, there is regular connectivity between the mesh cells and vertices.All the interior cells and vertices in the mesh have a fixed number of neighbors. Figure 2.1represents the topology for a structured mesh in which the interior cell represented by indices(i, j) has four neighbors represented by indices (i − 1, j), (i + 1, j), (i, j + 1) and (i, j − 1).These meshes usually consist of quadrilaterals in 2D and hexahedral elements in 3D.Figure 2.1: Structured mesh around leading edge of a NACA 0012 airfoilContrariwise, in an unstructured mesh, the mesh cells and vertices are not connected in anyregular pattern. The interior cells and vertices in the mesh do not have a fixed numberof neighbors. Figure 2.2 shows the topology for an unstructured mesh. Mesh resolutionand grading are important factors for unstructured meshes. The mesh typically needs tofiner in regions close to the solid boundary and coarser in the far-field. This is because the82.2. Mesh Propertiesgradients of flow variables are high close to the solid boundary as in the boundary layerand wake region and so finer mesh resolution is required in those regions of the domain forbetter solution accuracy. The coarser mesh in the far-field where changes are small helps toreduce the overall computational cost. Unstructured meshes usually consist of triangles orquadrilaterals in 2D and tetrahedra, prisms or pyramids in 3D [23, 54].Figure 2.2: Unstructured mesh around leading edge of NACA 0012 airfoilFor structured meshes, there is no need to store the data for each mesh vertex explicitly sincethe neighboring vertices in the physical space are also neighboring elements in the connectiv-ity matrix. This reduces the memory usage which is an advantage from the computationalpoint of view. For unstructured meshes, there is more memory usage due to the irregulartopology of unstructured mesh: the connectivity of each vertex in a mesh to other verticesneeds to be stored explicitly [27, 52].In structured meshes, due to the fixed topology, less CPU time is required by the finitevolume solver for evaluation of numerical fluxes, because the solution needed for flux evalua-tion can be computed by one-dimensional interpolation. For unstructured meshes, accurateflux evaluation requires polynomial reconstruction of solution (See Section 2.3) over eachcontrol volume which increases the CPU time.However, many industrial level problems involve complex geometries and generation ofstructured meshes around complex geometries is a major challenge. Although less CPU timeis required to compute the numerical flux, the drawback of fully structured meshes is thatwe don’t have the ability to put all the resolution where it is needed and we can end up92.2. Mesh Propertieswith excess resolution in some places. More specifically, we can’t independently control cellspacing in different parts of the domain and excess resolution increases CPU time. Also,generating a structured mesh around complex shapes can require time-consuming humaninterventions such as splitting the domain into multiple blocks depending on the complexityof the geometry. Such factors can offset the CPU time advantage gained in the evaluationof numerical fluxes by the solver. Unstructured meshes are much easier to generate aroundcomplex geometries and usually a preferred choice for complex configurations. Althoughthe polynomial reconstruction of the solution over each control volume increases CPU time,unstructured meshes have the ability to capture the flow features around complex geometrieswith fewer cells. This leads to a considerable reduction in overall CPU time. Also, it is easierto do adaptive mesh refinement since cell regularity is not an issue [32, 33].2.2.2 Isotropic and Anisotropic MeshesThe accuracy of numerical solutions is greatly influenced by the characteristics such as theshape, orientation and size of mesh elements. Figure 2.3a shows an isotropic mesh fragmentwith triangular elements in which the edge lengths of cells are approximately equal. In otherwords, triangular elements are roughly equilateral triangles. Contrariwise, an anisotropicmesh has cells with higher aspect ratio as depicted in Figure 2.3b.(a) Isotropic mesh (b) Anisotropic meshFigure 2.3: Isotropic and anisotropic mesh fragmentsThe choice of using an isotropic or an anisotropic mesh depends on the nature of the phys-102.3. Solution Reconstructionical processes occurring in the problem. If the physical changes are isotropic, it is desir-able to use an isotropic mesh with a cell size that equidistributes the error. However, theflow field in aerodynamics applications is usually accompanied by phenomena such as tur-bulence, boundary layers, wakes and shocks that have anisotropic physical behavior. Asmentioned earlier, resolving these flow features using isotropic elements is computationallyinefficient [14]. Moreover, isotropic meshes can provide only the proper mesh resolution whileanisotropic mesh adaptation gives both mesh resolution and proper cell alignment with thesolution to capture anisotropic flow features. Therefore, such complex flow features can beresolved more efficiently using anisotropic cells.2.3 Solution ReconstructionThe accuracy of the finite volume solution is highly dependent on the accuracy of fluxintegrals which in turn requires accurate approximation of the unknown variables in thecontrol volume. In order to approximate the unknown variables, represented by φ, in theflow field, the flow solver aims to replace the control volume average φi for each interiorcontrol volume in the mesh with a Taylor series about a control volume reference point(xi, yi) [39]:φi = φRi (x, y) = φ∣∣i+ ∂φ∂x∣∣∣∣i(x− xi) + ∂φ∂y∣∣∣∣i(y − yi)+ ∂2φ∂x2∣∣∣∣i(x− xi)22 +∂2φ∂x∂y∣∣∣∣i(x− xi)(y − yi)+ ∂2φ∂y2∣∣∣∣i(y − yi)22 + ..., (2.3)whereφi is the average value of the solution in control volume iφRi (x, y) is the value of reconstructed solution at a point (x, y) for control volume i∂k+1φi∂xk∂yiare the derivatives of the reconstructed solution for control volume iThe values φi and ∂k+1φi∂xk∂yicorrespond to coefficients of Taylor polynomials. The accuracy of112.3. Solution Reconstructionthe solution is one order higher than the polynomial degree. These coefficients are chosenso that we conserve the mean value of the solution in the control volume i, that is, we mustsatisfy the following criterion:1AiˆViφRi dA = φi (2.4)Combining Equations 2.3 and 2.4 leads to:φi =1AiˆViφRi dA = φ∣∣i+ ∂φ∂x∣∣∣∣ixi +∂φ∂y∣∣∣∣iyi +∂2φ∂x2∣∣∣∣ix2i2+ ∂2φ∂x∂y∣∣∣∣ixyi +∂2φ∂y2∣∣∣∣iy2i2 + ..., (2.5)with xnymi being the moments of area such that:xnymi =1AiˆVi(x− xi)n(y − yi)mdA. (2.6)The accuracy of a least-squares reconstruction scheme for smooth functions can be said,equivalently, to be k-exact or (k+ 1) order accurate and requires us to expand the modifiedTaylor series, represented by Equation 2.5, up to the k − th derivatives. While evaluatingthe derivatives, we try to minimize errors in the average value of φRi (reconstructed solution)for a nearby set of control volumes called the stencil {Vj}i.(a) Vertex-centered control volume (b) Cell-centered control volumeFigure 2.4: Reconstruction stencils [39]122.3. Solution ReconstructionIn the reconstruction stencil, we require at least as many control volumes as derivativeterms that we need to compute. However, practically we exceed this minimum number ofcontrol volumes for each order of accuracy, thus increasing the robustness of the least-squaresreconstruction scheme. For example, for second-order accuracy, we have three neighboringcontrol volumes and for third-order accuracy, we have nine neighboring control volumes. Theapproach to adding control volumes for both vertex-centered and cell-centered simulations isto add all the control volumes at the same topological distance to the reconstruction stencilat once. Figure 2.4 above shows the control volumes for the vertex-centered and cell-centeredcases with control volume labeled R representing the control volume to be reconstructed andthe numbers represent the order of accuracy at which each neighboring control volume isadded.The average value of φRi (reconstructed function) for each control volume in the stencilwill be:1AjˆVjφRi (−→x −−→x i)dA = φ∣∣i+ ∂φ∂x∣∣∣∣i1AjˆVj(x− xi)dA+∂φ∂y∣∣∣∣i1AjˆVj(y − yi)dA+∂2φ∂x2∣∣∣∣i12AjˆVj(x− xi)2dA+ ∂2φ∂x∂y∣∣∣∣i1AjˆVj(x− xi)(y − yi)dA+ ∂2φ∂y2∣∣∣∣i12AjˆVj(y − yi)2dA+ ... (2.7)In Equation 2.7, x−xi and y−yi are replaced with (x− xj)+(xj − xi) and (y − yj)+(yj − yi)to avoid computing moments for each control volume of {Vj}i about the reference point ofcontrol volume i. Equation 2.7 then becomes132.3. Solution Reconstruction1AjˆVjφRi (−→x −−→x i)dA = φi +∂φ∂x∣∣∣∣i(xj + (xj − xi)) + ∂φ∂y∣∣∣∣i(yj + (yj − yi))+ ∂2φ∂x2∣∣∣∣∣ix2j + 2xj (xj − xi) + (xj − xi)22+ ∂2φ∂x∂y∣∣∣∣∣i(xyj + xj (yj − yi) + (xj − xi) yj + (xj − xi) (yj − yi))+ ∂2φ∂y2∣∣∣∣∣iy2j + 2yj (yj − yi) + (yj − yi)22 + ... (2.8)The geometric terms in above equation are mesh dependent and have the following formx̂nymij ≡1AjˆVj((x− xj) + (xj − xi))n · ((y − yj) + (yj − yi))m dA=n∑l=0m∑k=0n!l! (n− l)!m!k! (m− k)! (xj − xi)k · (yj − yi)l · xn−kym−ljEquation 2.8 then becomes1AjˆVjφRi (−→x −−→x i)dA = φ|i +∂φ∂x∣∣∣∣ix̂ij +∂φ∂y∣∣∣∣iŷij (2.9)+ ∂2φ∂x2∣∣∣∣∣ix̂2ij2 +∂2φ∂x ∂y∣∣∣∣∣ix̂yij +∂2φ∂y2∣∣∣∣∣iŷ2 ij2 + · · ·Equation 2.9 gives the average value of the reconstructed solution φRi (−→x − −→x i) for controlvolume j. The difference between the actual control volume average and the average valueof the reconstructed solution can be easily accessed. The derivatives at −→x i are chosen tominimize the magnitude of error in predicting the control volume average values in a least-squares sense [39].142.4. Metric-Based Anisotropic Mesh AdaptationThe resulting least-squares system to compute the derivatives is1 xi yi x2i xyi y2i · · ·wi1 wi1x̂i1 wi1ŷi1 wi1x̂2i1 wi1x̂yi1 wi1ŷ2i1 · · ·wi2 wi2x̂i2 wi2ŷi2 wi2x̂2i2 wi2x̂yi2 wi2ŷ2i2 · · ·wi3 wi3x̂i3 wi3ŷi3 wi3x̂2i3 wi3x̂yi3 wi3ŷ2i3 · · ·.................. . . .wiN wiN x̂iN wiN ŷiN wiN x̂2iN wiN x̂yiN wiN ŷ2iN · · ·φ∂φ∂x∂φ∂y12∂2φ∂x2∂2φ∂x ∂y12∂2φ∂y2...i=φ¯iwi1φ¯1wi2φ¯2wi3φ¯3...wiN φ¯NIn the least-squares system, wij represents the geometric weights that are based on the dis-tance between control volume reference points. These geometric weights can be used tospecify the relative importance of neighboring control volumes. N is the number of nearbycontrol volumes in the reconstruction stencil. This constrained least-squares system is con-verted into an unconstrained least-squares system by eliminating the first row representingthe mean constraints by using Gauss elimination. The remaining unconstrained least-squaresproblem is then solved for each control volume using the singular value decomposition (SVD)method [15]. This is then followed by flux evaluation and flux integration, the details of whichcan be found in reference [3, 39].2.4 Metric-Based Anisotropic Mesh AdaptationAs stated in Chapter 1, the purpose of generating anisotropic meshes is to reduce the com-putational cost for a given level of solution accuracy. Anisotropic mesh cells help to achievethis goal by capturing the solution features with fewer cells. Moreover, unlike isotropic meshadaptation, anisotropic adaptation also provides good cell alignment in addition to the meshresolution in areas of interest.A robust anisotropic mesh adaptation scheme demands good error estimation and highquality anisotropic cells. Many researchers [4, 6, 12, 30] have illustrated that selectively refin-152.4. Metric-Based Anisotropic Mesh Adaptationing/coarsening the mesh cells and producing highly anisotropic meshes is effective. Differentmethods exist for error estimation. This thesis uses the anisotropic mesh error estimate im-plemented by Pagnutti [41]. This scheme makes use of local reconstruction error and givesan error measure based on the overall solution vector rather than any specific functional asin adjoint error estimation. Zuniga Vazquez’s [55] anisotropic mesh adaptation scheme wasused in this research. The robustness of the mesh adaptation scheme is further demonstratedwith different test cases in this thesis.The concept behind metric-based adaptation is to produce a unit mesh in Riemannianmetric space. The error estimate gives the desired anisotropy which is stored in a metrictensor at each vertex in the computational domain. The desired anisotropy is then commu-nicated to the mesh adaptation code to produce an anisotropic mesh. For 2D problems, themetric is a 2×2 symmetric positive definite tensor. Geometrically, the metric at a point tellshow distance and angles are measured as seen from that point. If we have a metric tensorMp, then the deformation tensor Lp can be defined as:Mp = LTp Lp and det (Lp) > 0 (2.10)with Lp being an invertible matrix that maps physical space into metric space as depictedin Figure 2.5 and the distances and angles are measured in the Euclidean way as seen bypoint p′.Lp'Physical space Metric spaceFigure 2.5: Deformation of tensor Lp162.4. Metric-Based Anisotropic Mesh AdaptationThe ellipses around point p′ in Figure 2.5 are equidistant from each other in metric spaceand define the anisotropy using three components: radius r1, radius r2 and angle θ as shownin Figure 2.6. Then, the metric is represented by these three components as:M = RΛRT =cos θ − sin θsin θ cos θ1/r21 00 1/r22cos θ sin θ− sin θ cos θ (2.11) θ r r21Figure 2.6: Anisotropy components for a certain pointA metric is a distance function d : X × X → R that defines the positive distance measurebetween any two distinct points in a vector space X and fulfills the following conditions:d(x, y) > 0 ⇐⇒ x 6= y non-negativityd(x, y) = d(y, x) symmetryd(x, z) ≤ d(x, y) + d(y, z) triangle inequality (2.12)In the case of a uniform isotropic mesh, the metric is the standard Cartesian distancedstd(x, y) =√(x− y)T (x− y). The metric for anisotropic adaptation dani(x, y) can beobtained by the invertible linear transformation L of dstd(x, y).172.4. Metric-Based Anisotropic Mesh Adaptationdstd(x, y) =√(x− y)T (x− y)dani(x, y) =√(L(x− y))T (L(x− y))dani(x, y) =√(x− y)T (LTL) (x− y)dani(x, y) =√(x− y)TM(x− y) (2.13)This illustrates that the metric dani(x, y) can be represented by a positive definite matrixM = LTL. Further details on the complete metric computation and its suitability to generatemeshes can be found in references [41, 49].2.4.1 Aspect ratio and orientation from metricThe metric is given by:M = Mxx MxyMxy MyyFrom Mohr’s circle, we have:Mmax =Mxx +Myy +√(Mxx −Myy)2 + 4M2xy2Mmin =Mxx +Myy −√(Mxx −Myy)2 + 4M2xy2The orientation is given by:θ = 12arctan(2MxyMxx −Myy)(2.14)The aspect ratio is given by:AR =√MmaxMmin(2.15)182.5. Operations for Mesh Quality Improvement2.5 Operations for Mesh Quality ImprovementThe anisotropic mesh adaptation scheme used in this research uses four principal operationsfor mesh quality improvement. These include face swapping in metric space, vertex insertionbased on quality and maximum metric length, vertex removal (coarsening) based on theminimum metric length and vertex movement (smoothing) based on quality. The first threeoperations are implemented by GRUMMP1 [38], our in-house meshing code, while vertexmovement is done using Mesquite [24].2.5.1 Face SwappingFace swapping changes the connectivity of vertices in the mesh but the number of meshvertices and cells remain unchanged. The new connectivity is decided such that the maximumangle in the metric space is minimized and new triangles with better quality are obtained.Figure 2.7 shows the triangles before and after face swapping. The necessary condition forface swapping is:Quality (V1V2V3) +Quality (V1V3V4) < Quality (V1V2V4) +Quality (V2V3V4)V3V4V2 V1(a) Initial configuration V3V4V2 V1(b) Modified configurationFigure 2.7: Face swapping in metric space1Generation and Refinement of Unstructured Mixed-Element Meshes in Parallel192.5. Operations for Mesh Quality Improvement2.5.2 Vertex InsertionThe vertex insertion or edge splitting is done based on two criteria and an insertion queue isused for this operation. The first criterion is based on quality and so the vertex is insertedonto the edge if it will improve the quality and higher priority in the insertion queue is givento the cells with the worst quality. The second criterion is related to maximum metric lengthand so edges longer than a certain length (based on the metric) are split and higher priorityis given to the edge with longest metric length. However, to prevent the creation of extraedges, an edge is not placed in the insertion queue if it will result in an edge smaller thanthe desired minimum length in the metric space. This means that the desired length is suchthat:lminM ≤ lM ≤ lmaxMThe vertex insertion for an interior and boundary edge is shown in Figures 2.8 and 2.9respectively. This operation increases mesh vertices by one and cells by two for interioredges while vertices and cells both by one for boundary edges. The same criterion is usedfor both the boundary and interior edges. After vertex insertion, face swapping is done asrequired to improve the mesh quality.V3 V2V1 V4(a) Initial configurationV3 V2 V1 V4Vnew(b) Modified configurationFigure 2.8: Vertex insertion for an interior edge202.5. Operations for Mesh Quality Improvement V3V2 V1 (a) Initial configuration V3 V2 V1Vnew(b) Modified configurationFigure 2.9: Vertex insertion for a boundary edge2.5.3 Vertex RemovalThe third mesh modification technique is vertex removal which coarsens the mesh andthereby helps to reduce the computational cost. We get rid of edges that have a smallmetric length. Again a priority queue is used and higher priority is given to the edge withthe shortest metric length. This means that the desired length is such that:lM ≥ lminMThis operation is again followed by face swapping as required to improve mesh quality.Figures 2.10 and 2.11 shows the process of vertex removal for an interior and boundary edgerespectively.V3 V2V1 V4V5V6(a) Initial configurationV3 V2V1 V4V5Removed(b) Modified configurationFigure 2.10: Vertex removal for an interior edge212.6. StreamlinesV3V2 V1 V4(a) Initial configurationV3V2 V1 Removed(b) Modified configurationFigure 2.11: Vertex removal for a boundary edge2.5.4 Vertex MovementThe vertex movement operation is based on the Target-Matrix Optimization Paradigm(TMOP) as used in Mesquite [7, 24]. The mesh vertices are moved or repositioned in such away that an optimal mesh of superior quality is obtained based on the target matrix. Thetarget matrix depends on the shape of the target element. In the case of anisotropic meshes,the target element should be a right triangle to produce quasi-structured meshes [49]. Theshape (aspect ratio and orientation) of each cell in an anisotropic mesh is determined fromthe metric and so the metric can also be used to calculate the target matrix. Since themetric information is stored at each vertex of the mesh element, an average metric is usedto calculate the target matrix for each cell. More details on TMOP and average metriccalculation can be found in Sharbatdar’s thesis [49].Further details on above mesh quality operations can be found in [41, 49, 55].2.6 StreamlinesA streamline can be defined as a path traced out by a massless particle propagating with thefluid velocity. The streamline is tangent everywhere to the local velocity field or instanta-neous velocity. Therefore, we need to consider one instant in time; for unsteady flow, thesestreamlines can look very different at each instant in time as they represent the direction ofthe flow (there is no flow across streamlines) [1]. This is what makes streamlines immenselyuseful for investigating different types of flow. The velocity fields evaluated numerically can222.7. Runge-Kutta Method and Adaptive Stepsize Controlbe called CFD velocity fields and these discrete CFD velocity fields, defined on a discretecomputational domain, are actually an approximation to the exact mass conservative veloc-ity fields [26]. Figure 2.12 shows streamlines around the NACA 0012 airfoil at Re = 5000,Ma = 0.5 and α = 0°.Figure 2.12: Streamlines around NACA 0012 airfoil at Re = 5000, Ma = 0.5 and α = 0◦It is possible to obtain a streamline equation from the velocity field by using velocity vectorsin each direction in the flow. For a two-dimensional velocity field, this can be expressed indifferential form as:dxdt= u (x, y) (2.16)dydt= v (x, y) (2.17)Equations 2.16 and 2.17 then need to be solved simultaneously. We can integrate themnumerically to obtain x (t) and y (t) along the streamline.2.7 Runge-Kutta Method and Adaptive Stepsize ControlThe velocity data obtained from solution reconstruction is integrated numerically to obtainthe streamline. The Runge-Kutta method is one of the classical methods used to numeri-cally integrate ordinary differential equations (ODE’s). Several versions of the Runge-Kuttamethod are available but in this thesis, we use the fourth-order accurate Runge-Kutta method232.7. Runge-Kutta Method and Adaptive Stepsize Control(RK4) to numerically integrate the ODE representing the velocity field. For an initial valueproblem such as:x˙ = f (x, t) x (t0) = x0 (2.18)RK4 can be expressed as:xi+1 = xi +h6 (k1 + 2k2 + 2k3 + k4) +O(h5)(2.19)ti+1 = ti + hwhere:k1 = f (xi, ti)k2 = f(xi +k1h2 , ti +h2)k3 = f(xi +k2h2 , ti +h2 ,)k4 = f (xi + k3h, ti + h)and xi+1 is the numerical approximation of x (ti+1) [50].For each step, Equation 2.19 requires four computations of the right-hand side. Thesecomputations are the values of the derivative at a point: k1 is the value at the starting pointof the interval, k2 and k3 are the values at the first and second trial midpoints respectivelyand k4 at the trial endpoint. The final function in Equation 2.19 is then obtained by takingthe average of the four derivatives with higher weight given to the derivatives at the trialmidpoints.An efficient differential equation integrator should have some stepsize control based on anerror estimate. Adaptive stepsize control can be implemented for the RK4 scheme based onan estimate of the local truncation error. In this case, each step is taken twice, once as a full242.7. Runge-Kutta Method and Adaptive Stepsize Controlstep and then as two half steps and the difference between the two numerical estimates is thetruncation error. Figure 2.13 illustrates the adaptive stepsize control for RK4, with filledsquares representing the points at which the derivatives are evaluated. The big step andtwo small steps procedure share the starting point represented by the empty square for thesmall step so no extra elevation is required at that point and the total number of evaluationsfor two small steps will be 11. Since reducing the step to half also drops the error, the 11evaluations per two small steps should be compared with 8 and not 4 evaluations of the bigstep. Therefore, the factor for the additional cost is 1.375.If we denote the exact solution by x (t+ 2h) while advancing from t to t + 2h and thenumerical estimates by x1 and x2, then [44]:x (t+ 2h) = x1 + (2h)5 φ +O(h6)x (t+ 2h) = x2 + 2 (h)5 φ +O(h6)Truncation error ≡ x2 − x1If the estimated error is above the user-specified tolerance for the truncation error, then thestep is repeated and stepsize is cut in half until the tolerance criteria is satisfied. Similarly, ifthe error is too small compared to the tolerance, then the step is repeated and the stepsize isdoubled to save the computational cost. This process provides optimal stepsize automaticallyand the user does not need to spend time on finding the optimal stepsize each time.big steptwo small stepsFigure 2.13: Adaptive stepsize control for RK4 [44]25Chapter 3MethodologyThe approach for anisotropic mesh adaptation is divided into several steps. This chaptergives an overview of each step and summarizes how our mesh adaptation loop works. In theend, it discusses the advancing layer mesh generation which is used as an alternate techniqueto generate an anisotropic quad dominant mesh.3.1 Wake Centerline TrackingThe algorithm to track the streamlines is implemented in our in-house solver code, ANSLib[40]. For 2D airfoil flows, we can find the wake centerline by tracking a streamline from thetrailing edge. A massless particle propagating with the fluid velocity is released in the flowfield and its location tracked in time to see where the streamline goes, as the streamline istangent to the velocity vector. Note that we do not integrate the flow in time but rather themotion of particle within the steady flow. The velocity data is obtained from the solutionreconstruction using the vertex-centered finite volume solver. For linear reconstruction, theODE for the streamline location is:ddt xy = uv = u0 + a1 (x− xi) + a2 (y − yi)v0 + b1 (x− xi) + b2 (y − yi) (3.1)In the above equation, u0 and v0 are the x and y components of velocity respectively at theinitial point (x0, y0) where the streamline enters the control volume. Moreover, a1, a2, b1,and b2 are the derivatives of the reconstructed solution at any point (xi, yi) in the controlvolume i.263.1. Wake Centerline TrackingWe integrate this system of ODE’s in time using an adaptive fourth-order Runge-Kuttamethod. A geometric search tree is used to identify which control volume the particle isin based on its position at a given time. The adaptive time step scheme for numericalintegration ensures that the local truncation error remains within a user-specified tolerance(see Section 2.7).Furthermore, we get a steady separation in the wake at different flow conditions. Figure3.1 shows the streamlines for a separated flow. In case of flow separation, numericallyintegrating Equation 3.1 forward in time is not sufficient. To track the wake centerline,we need to march downstream from the trailing edge and so x (t+ ∆t) should be alwayshigher than the previous value x (t). Simply doing forward integration in time from thetrailing edge can sometimes take us inside in the airfoil (for example at Reynolds numberRe = 5000, Mach number Ma = 0.5 and angle of attack α = 0°) rather than downstreamfrom the trailing edge. In such a case, we need to start with backward integration in time tomarch in the right direction. Therefore, in the streamline tracking algorithm, we always startwith backward integration in time and if it takes us upstream (the next value of x (t) willbe less than 1), we exit the backward integration cycle immediately at the trailing edge anddo only forward integration in time to move downstream. However, if backward integrationtakes us downstream from the trailing edge, then there is a saddle point along the wakecenterline and the backward integration cycle is terminated once we hit the saddle point.An offset is added to escape the separation bubble, followed by forward integration in timeto move downstream from the separated region. The offset that needs to be added is usuallybetween 1% and 3% of the chord length but can be slightly higher — 5% to 6% of chordlength — when separation is large at higher angles of attack. Moreover, one exception is thatthere can be a saddle point in the flow-field along the wake centerline even when forwardintegration in time initially takes us downstream from the trailing edge. In that case (forexample at Re = 5000, Ma = 0.5 and α = 10°) two cycles of forward integration in timeare required. We exit the first cycle of forward integration in time as soon as we hit thesaddle point and the value of x (t) starts to decreases after that point. An offset is added273.2. Wake Centerline Insertion into Unstructured Meshto escape the separation bubble and we then do second forward integration in time to movedownstream again.Because streamline integration does not exactly hit the saddle point, we identify thesaddle point as the first point where |v| > |u| and the value of x (t) obtained from numericalintegration is higher than the previous value of x (t). While inexact, this is a very closeapproximation to the saddle point location. As illustrated by Figure 3.1, for separated flows,the wake centerline can be tracked in several steps: (1) backward or forward integration intime from the trailing edge, (2) find the saddle point along the wake centerline, and (3) addan offset from the saddle point to escape the separation bubble and do forward integrationin time to continue downstream.1 2 3trailing edgesaddle pointFigure 3.1: Separated flow3.2 Wake Centerline Insertion into Unstructured MeshThe obtained streamline is inserted as an internal boundary into the existing isotropic un-structured mesh, using the surface insertion algorithm implemented in GRUMMP by Zaideand Ollivier-Gooch [60]. Inserting the streamline as an internal boundary helps to deal withphysical behavior in the wake and enables the pre-existing mesh to adapt to that behavior.The surface insertion approach starts with two initial inputs: an initial unstructured meshand information about the geometry to be inserted. This approach has no requirement onthe initial mesh and rather than moving the existing vertices onto the surface, new vertexlocations on the boundary curve are computed and points are inserted to match the surface.The algorithm samples the boundary curve on the mesh, using the length scale of the meshto determine edge lengths. Vertices in the original mesh near the curve are removed, the283.3. Anisotropic Mesh Adaptationnew vertices on the sampled curve are inserted into the mesh and surface recovery is done.The surface insertion technique gives freedom to work with an arbitrary initial mesh butthe main point of this approach is to recover the surface rather than moving points whilenot changing the mesh at all more than a cell size or two from the wake surface. Once thesurface is recovered after insertion, the mesh quality is improved with swapping, smoothingand refinement. Further details on the surface insertion technique can be found in the paperby Zaide and Ollivier-Gooch [60].3.3 Anisotropic Mesh AdaptationAfter inserting the wake centerline into the existing unstructured mesh, the next step is todo metric-based anisotropic mesh adaptation. As explained earlier, our metric is based onsolution approximation error, and the error estimation is done by interpolation of the solutionbetween known values. For a second-order accurate solution, the solution approximationerror12fxx (x, y)x2 + fxy (x, y)xy +12fyy (x, y) y2,can be related to the Hessian of the solutionH =∂2f∂x2∂2f∂x∂y∂2f∂y∂x∂2f∂y2A metric function is chosen such that it gives a good approximation of the solution error. Itis approximated by doing a Fourier transform around a unit circle of the maximum error inany solution component and the Fourier coefficients are obtained with the help of numericalintegration. The combinations of Fourier coefficients are then used as metric terms [42, 41].293.4. Mesh MorphingAs explained in Section 2.5, the anisotropic mesh adaptation scheme uses four operationsfor mesh quality improvement in the metric space but it preserves the wake centerline. Theadapted mesh is used to compute the flow solution again, using the second-order accuratevertex-centered finite volume scheme and the velocity field is integrated to obtain the newwake centerline.3.4 Mesh MorphingIn the past, researchers [21, 58] have used the linear elasticity method to deform the meshin order to take into account the effect of curvature especially for higher-order solvers toaccurately represent the boundary faces. Jalali and Ollivier-Gooch [21] implemented thisapproach in our in-house code for curving the boundary faces and solved the linear elasticityproblem using the finite element discretization. In this thesis, we extend their approachfurther and take into account linear elements and change the boundary condition to displacethe wake centerline.In the second and subsequent mesh adaptation cycles, we begin by morphing the meshto match the new wake centerline. We project the mesh points along the old wake centerlineonto the new wake centerline and solve a linear elasticity problem with displacements betweenthe two wake centerlines as constraints to morph the rest of the mesh. After this, we applymetric-based adaptation to the morphed mesh. Because the morphing process preservesexisting adaptation in the wake region, the metric adaptation is more efficient. The linearelasticity equation is given by:∂∂x(d11∂δx∂x+ d12∂δy∂y)+ ∂∂y(d33(∂δx∂y+ ∂δy∂x))= 0∂∂x(d33(∂δx∂x+ ∂δy∂y))+ ∂∂y(d21∂δx∂y+ d22∂δy∂x)= 0 (3.2)where δ = (δx, δy) are the nodal displacement vectors. The coefficients d are based on303.5. Algorithm OverviewYoung’s modulus E and Poisson’s ratio ν as follow:d11 = d22 =E (1− ν)(1 + ν) (1− 2ν)d12 = d21 =Eν(1 + ν) (1− 2ν) (3.3)d33 =E2 (1 + ν)We used Poisson’s ratio equal to 0.25. Young’s modulus was assumed to be constant through-out the computational domain and therefore cancels in Equation 3.2. The complete detailson the finite element discretization of Equation 3.2 can be found in reference [20].3.5 Algorithm OverviewThe metric-based anisotropic mesh adaptation algorithm is summarized here and the meshadaptation loop is depicted in Figure 3.2.1. We start by generating an initial isotropic unstructured mesh using the isotropic meshgenerator in GRUMMP.2. The second step is to compute the flow solution on the initial unstructured mesh. Thevelocity field is obtained from solution reconstruction and integrated numerically usingthe adaptive fourth-order Runge-Kutta method to obtain the wake centerline.3. The wake centerline is then inserted as an internal boundary into the initial unstruc-tured mesh using the surface insertion algorithm as explained in Section 3.2.4. The metric is then computed at each mesh vertex.5. The mesh and metric are passed to the anisotropic mesh adaptation algorithm in themesh generator to produce an anisotropic mesh. The mesh adaptation scheme usesface swapping, vertex insertion, vertex removal and smoothing operations to improvethe mesh quality but preserves the wake centerline without moving it at all.313.5. Algorithm Overview6. The adapted mesh is used to re-compute the flow solution and a new wake centerlineis obtained.7. We then morph the adapted mesh to match the new wake centerline by solving thelinear elasticity equation as discussed in Section 3.4. The process from step 4 to 6in Figure 3.2 is then repeated. If the aerodynamic coefficients (CL, CD) are not gridconverged, we move to step 7 and repeat steps 4, 5 and 6 until the aerodynamiccoefficients are converged. Step 2: Solve the problem to obtain wake centerline Step 3: Insert wake centerline into the mesh Step 5: Metric-based anisotropic mesh adaptation Step 6: Re-compute solution Yes Final mesh Repeat from step 4 for next cycle No Step 1: Create an initial mesh Is solution converged? Step 4: Compute metric at each vertex Step 7: Morph old mesh to match new wake centerline Figure 3.2: Algorithm overview (mesh adaptation loop)323.6. Advancing Layer Mesh Generation3.6 Advancing Layer Mesh GenerationFinally, another scheme is used to compare with the results obtained on the final adaptedmesh using the aforementioned method. In this scheme, we use the advancing layer meshgenerator of Numerow and Ollivier-Gooch [37]. We define our initial surface to include theairfoil and the discretized wake centerline as shown in Figure 3.3, and then march off thesurface into the domain interior layer by layer. The wake centerline used is obtained bynumerically integrating the velocity data from the initial solution on isotropic unstructuredmesh. The advancing layer mesh generator produces mixed element meshes that are quaddominant and have a few triangular cells. This leads to several advantages such as: (1) lowercell count, (2) better alignment and orthogonality, (3) more geometrically similar cells, and(4) excellent resolution of the boundary layer and wake.Figure 3.3: Initial geometry to generate advancing layer meshThe advancing layer mesh is generated taking into account the flow physics. We consider theBlasius boundary layer solution for a flat plate and the two-dimensional asymptotic solutionin the wake.(x) cDeveloped wakeWakeDeveloping wakeU0U0U0U0U0uu1< 3c). As a result,Equation 3.6 is also good only for the far-wake. Therefore, to generate the quad dominantadvancing layer mesh, we take the average value of Equations 3.4 and 3.6 and consider aReynolds number based on x, giving us Equation 3.7. Using Equation 3.7 ensures that δ (x)increases as we move from the leading edge of the airfoil to the trailing edge and then in thewake region downstream.343.6. Advancing Layer Mesh Generationδ (x) = 3.5x√Rex(3.7)Equation 3.7 gives an approximate thickness δ (x) for the shear layer both on the airfoiland in the wake. In order to properly resolve the boundary layer and wake, many layers(between 20 and 30 layers) of mesh cells are required across the thickness δ (x) and sothe off-wall spacing in the mesh should be much smaller than δ (x). Figure 3.5 shows thenumber of layers across the thickness δ (x) with variable spacing. Increasing the number oflayers increases the resolution within that region and so we can derive a formula to relatethe thickness δ (x), number of layers and the stretching factor for variable spacing betweenlayers. The thickness δ is given by:δ =N∑n=1hn= dymin(1 + s+ s2 + ...+ sN−1)= dyminN∑n=1sn−1= dymin(1 + s+ s2 + ...+ sN−1)(1− s)(1− s)δ = dymin(sN − 1)(s− 1) (3.8)Then from Equations 3.7 and 3.8, we havedymin =3.5x√Rex(s− 1)(sN − 1) (3.9)where dymin is the off-wall spacing, s is the stretching factor (s > 1) and N is the numberof layers across the thickness δ. As per Equation 3.9, off-wall spacing is based on x and theReynolds number Rex.353.6. Advancing Layer Mesh Generationdymindymindymins2dymins3sFigure 3.5: Number of layers across thickness δ (x)36Chapter 4ResultsIn this chapter, we first present results to verify the correctness of the streamline tracking andmesh morphing algorithms. Later, different laminar flow cases are included to demonstrateour approach for mesh adaptation. For each of these test cases, the actual velocity data fromthe flow solver was used and integrated numerically. The solution was computed based on thesecond-order accurate vertex-centered finite volume solver as discussed by Ollivier-Gooch inhis paper [40]. This flow solver makes use of least squares reconstruction [39], Roe’s scheme[46] and the Newton-GMRES method [34, 36] for fast steady-state convergence. Moreover,the results of test cases from adapted triangular meshes with and without the wake centerlineand advancing layer meshes are compared.The laminar flow results include both attached and separated flow cases. The flowconditions for the first two cases are: (1) Re = 2100, Ma = 0.5, α = 0◦ and (2) Re = 2100,Ma = 0.5, α = 2◦ and for these cases, we choose Re = 2100 as the flow is attached at theseconditions. Moreover, the third test is chosen to show that our algorithm works for separatedflows as well and also there are numerical results available by well-established CFD codesfor Re = 5000, Ma = 0.5 and α = 0◦ .In the end, we also present a test case for turbulent subsonic flow. For turbulent flow,the solution was computed using the second-order accurate cell-centered finite volume solverand the velocity field was integrated to track just the wake centerline at different angles ofattack, without doing any mesh adaptation.374.1. Verification Cases4.1 Verification Cases4.1.1 Streamline Tracking AlgorithmThe streamline tracking algorithm was tested using an analytical velocity field given by:ddt(xy)=(−2yx)(4.1)Integrating the above velocity field analytically gives streamlines that are ellipses. Theproblem was solved numerically on an isotropic unstructured mesh around the NACA 0012airfoil. Figure 4.1 shows the mesh, which was generated using the isotropic mesh generator inGRUMMP [38]. The far-field boundary was a circle with radius 500c centered at the leadingedge of the airfoil. The trailing edge of the airfoil, (1, 0), was chosen as the starting pointto track the streamline. The velocity data was control volume averaged onto the mesh, thenreconstructed. Because the velocity is linear, this reconstruction was exact. Nevertheless,this test verifies our ability to track streamlines through a mesh using reconstructed velocitydata.(a) Far-field view (b) Closeup viewFigure 4.1: Isotropic unstructured mesh around a NACA 0012 airfoil with 7310 cells and3871 verticesThe streamline was integrated halfway around the ellipse. Finally, to validate the result, thecase was compared to the Paraview streamtracer result as shown in Figure 4.2.384.1. Verification CasesFigure 4.2: Streamline tracking for an analytical velocity fieldNext, the order of accuracy of the numerical method was tested based on the streamlinedata obtained from numerical integration. The velocity field was integrated numericallyusing three different time-steps to obtain three sets of data for streamline location. Then,the error between the actual and numerically predicted location on the streamline at aninstant in time was calculated as:error =√(x− xn)2 + (y − yn)2where x and y are the analytical values while xn and yn are numerical values.Time-Step error ratio( (error)4t=h(error)4t=h/2)0.05 4.601× 10−7 −0.025 2.870× 10−8 16.0310.0125 1.800× 10−9 15.944Table 4.1: Error in solution (values at t = 2.20)Table 4.1 illustrates that the error drops by 16 times when the time-step is reduced by half sothe Runge-Kutta method is fourth-order accurate as expected. Moreover, Figure 4.3 showsthe error in solution at each point along the curve for three different time-steps mentionedabove. The ratio of error between the green and red curves and between the red and bluecurves at each point was found to be almost exactly 16.394.1. Verification CasesFigure 4.3: Error in solution at each point along the curve4.1.2 Mesh Morphing AlgorithmTo verify the correctness of the mesh morphing approach, we used the method of manufac-tured solutions. The nodal displacement vectors were selected as:δx = sin(pix+ pi2)sin(piy + pi2)(4.2)δy = sin(pix+ pi2)sinh(piy + pi2)The manufactured solution is substituted into the linear elasticity equations (Equation 3.2)and a source term vector was identified (Equation 4.3).S = SxSy (4.3)whereSx = d11sin(pix+ pi2)pi2sin(piy + pi2)− d12cos(pix+ pi2)pi2cosh(piy + pi2)−d33(−sin(pix+ pi2)pi2sin(piy + pi2)+ cos(pix+ pi2)pi2cosh(piy + pi2))404.1. Verification CasesandSy = −d33(cos(pix+ pi2)pi2cos(piy + pi2)− sin(pix+ pi2)pi2sinh(piy + pi2))−d21cos(pix+ pi2)pi2cos(piy + pi2)− d22sin(pix+ pi2)pi2sinh(piy + pi2)The discrete problem was solved with this source term. The source term vector was inte-grated using the Gauss quadrature rule. The L2 norm of error was then computed as thedifference between the manufactured solution of Equation 4.2 and the linear finite elementnumerical solution. Figure 4.4 shows the variation of L2 norm of error in solution withthe mesh size and confirms that the numerical solution is second-order accurate for linearelements, as expected.xxy+ yorder = 1.97order = 1.94order = 2.00Figure 4.4: L2 norm of error for linear elasticityThe correct implementation of boundary condition to displace the wake centerline was alsotested by inserting an imaginary horizontal line at the trailing edge of the airfoil in anunstructured mesh and then displacing it from 0° to 3° as shown in Figure 4.5. This testconfirms that the mesh morphing algorithm displaces the inserted line as required and morphs414.2. Laminar Flow Test Casesthe mesh accordingly.(a) Line at 0° (black) and 3° (red)(b) Mesh before (black) and after (red) mesh morphingFigure 4.5: Mesh morphing to displace horizontal line from 0◦ to 3◦4.2 Laminar Flow Test Cases4.2.1 Viscous flow at Re = 2100, Ma = 0.5 and α = 0◦Our first adaptation test is subsonic viscous flow at Reynolds number Re = 2100, Machnumber Ma = 0.5 and angle of attack α = 0°. The viscous terms are discretized as discussedby Ollivier-Gooch and Van Altena [39]. Using the mesh in Figure 4.1, we computed the initialflow solution and integrated the velocity field to track the streamline from the trailing edgeof NACA 0012 airfoil. The streamline was then inserted as an internal boundary into thesame unstructured mesh using the surface insertion algorithm. It can be seen from Figures4.6a2 and 4.6b below that the result of Paraview streamtracer matches with the streamlinetracking algorithm result. The wake centerline is not exactly straight but the oscillations areof the order 10−3 to 10−5 chords from the symmetry axis; these oscillations become smalleron more refined meshes and as we go further downstream from the trailing edge.2In Figure 4.6a, red line is the position of streamline obtained from Paraview streamtracer and so it is nota part of the mesh.424.2. Laminar Flow Test Cases(a) Paraview streamtracer result (Mesh before streamline insertion with 7310 cells) - Viscous flow - Red line isthe streamline(b) Streamline tracking algorithm result (Mesh after streamline insertion with 7468 cells) - Viscous flowFigure 4.6: Streamline tracking at Re = 2100, Ma = 0.5 and α = 0◦A metric was computed for the mesh in Figure 4.6b and passed to the anisotropic meshgeneration algorithm in GRUMMP to do mesh adaptation. The adapted mesh was used tocompute the flow solution again and obtain the new wake centerline. We then morphed theadapted mesh to match the new wake centerline by solving the linear elasticity problem asdiscussed in mesh morphing section in Chapter 3. Figure 4.7 shows the mesh before andafter mesh morphing. This process was repeated for subsequent cycles of anisotropic meshadaptation and helped to improve the solution accuracy. The flow solution in Figure 4.8aon the initial mesh shows that the mesh is not well-suited to compute the viscous flow andresolve the boundary layer and wake. However, after each adaptation, better resolution andmesh alignment is obtained in the boundary layer and wake as shown in Figure 4.8.The position of the wake centerline on the initial and final adapted mesh is depicted inFigure 4.9. The wake centerline from the final adapted mesh is closer to the symmetry axis,illustrating that the final mesh can track the position of wake better for this symmetric flow.The small variations are still present since the mesh is not completely symmetric.Figure 4.7: First adapted mesh (black) and morphed mesh (red) before second adaptation434.2. Laminar Flow Test Cases(a) Solution on initial mesh (3871 vertices)(b) Solution after first adaptation (6048 vertices)(c) Solution after second adaptation (8268 vertices)(d) Solution on final mesh (12121 vertices)Figure 4.8: Adapted meshes and flow solutions with wake centerline at Re = 2100, Ma = 0.5and α = 0◦444.2. Laminar Flow Test CasesFigure 4.9: Wake centerline position - initial (black) and final (red) at Re = 2100, Ma = 0.5and α = 0◦- Closeup view at the trailing edgeAfter obtaining the solutions on triangular meshes, an anisotropic quad dominant mesh wasgenerated using the advancing the layer mesh generator in GRUMMP. The initial geometryused to generate the advancing layer mesh included the NACA 0012 airfoil and the discretizedwake centerline obtained from an initial solution in Figure 4.8a.Figure 4.10: Initial geometry for advancing layer mesh at Re = 2100, Ma = 0.5 and α = 0◦(a) Solution on final adapted triangular mesh (12121 vertices)(b) Solution on advancing layer quad dominant mesh (5558 vertices)Figure 4.11: Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 2100, Ma = 0.5 and α = 0◦ , showing better solution behavior in the wake for theadvancing layer mesh454.2. Laminar Flow Test CasesFigure 4.11 shows the solution on the final adapted triangular mesh and the quad dominantadvancing layer mesh, illustrating that the advancing layer mesh with many fewer degreesof freedom can track the wake for a longer distance compared to adapted triangular meshes.The convergence history of aerodynamic coefficients with respect to degrees of freedomfor triangular meshes with and without the wake centerline and for the advancing layer quaddominant meshes is displayed in Figure 4.12. For triangular meshes, it can be observed thatthe lift coefficient for meshes with the wake centerline included converges slightly faster,while the convergence of drag coefficients is effectively identical. Moreover, the solution isconverging to a single value on meshes with and without the wake centerline. By compari-son, for the advancing layer meshes, the aerodynamic coefficients converge faster comparedto both triangular meshes. Especially, the lift coefficient on the advancing layer meshesconverges much faster and has a smaller error as we expect to get CL = 0 for this symmetricflow. The difference in drag coefficient from the refined advancing layer quad dominant meshand the grid converged value using triangular meshes is just 0.26%.In Figure 4.13, we present the convergence history against total CPU time for meshingand solver. The total CPU time for triangular meshes with a wake centerline is somewhatless compared to meshes without the wake centerline. There is clearly a CPU time advantagefor the advancing layer quad dominant meshes compared to both triangular meshes. Themain reason is that it is much cheaper to generate advancing layer mesh compared to doingmesh adaptation and this reduces CPU time required for meshing. Also, the advancing layermesh can track the wake much better with fewer degrees of freedom which reduces the CPUtime needed by the solver.464.2. Laminar Flow Test Cases(a) Lift coefficient(b) Drag coefficientFigure 4.12: Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 2100, Ma = 0.5 and α = 0◦474.2. Laminar Flow Test Cases(a) Lift coefficient(b) Drag coefficientFigure 4.13: Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 2100, Ma = 0.5 and α = 0◦484.2. Laminar Flow Test Cases(a) Two chords downstream of the trailing edge (b) Five chords downstream of the trailing edge(c) Eight chords downstream of the trailing edgeFigure 4.14: Velocity profiles in the wake for triangular and advancing layer quad dominantmeshesIn Figure 4.14, we present the comparison of velocity profiles in the wake for triangularmeshes with wake centerline and the advancing layer quad dominant meshes at differentlocations downstream from the trailing edge. Figures 4.14a to 4.14c illustrate that there is avelocity deficit in the wake as expected. Furthermore, we expect to get symmetric velocityprofiles in the wake for this case. At each location (2, 5 and 8 chords downstream the trailing494.2. Laminar Flow Test Casesedge), the velocity profiles from the advancing layer meshes are almost symmetric. Thevelocity profile from the initial triangular mesh is not physical due to poor mesh resolutionbut as we do mesh adaptation, the velocity profile improves with each mesh refinement andapproaches the velocity profiles from the advancing layer meshes. At two chords downstreamof the trailing edge, the final adapted mesh solution is nearly as smooth and symmetric as thesolution on the advancing layer meshes. Note that compared to the advancing layer meshes,the final adapted triangular mesh has much more the degrees of freedom. Also, the velocityprofile from the adapted meshes gets worse in comparison as we move further downstreamfrom the trailing edge from two chords to eight chords. Therefore, this comparison showsthat the advancing layer meshes can resolve the wake more efficiently. The velocity profilesfrom the advancing layer quad dominant meshes with 5558 and 7806 vertices are identical,illustrating that increasing the mesh resolution does not change the velocity profiles further.4.2.2 Viscous flow at Re = 2100, Ma = 0.5 and α = 2◦The second test case is viscous flow around the NACA 0012 airfoil at a non-zero angleof attack. All other flow conditions and initial mesh used were similar to first test case.Figure 4.15 shows the streamline tracking algorithm result for this case. The steps used foranisotropic mesh adaptation were same but here four cycles of adaptation were required forgrid convergence of aerodynamic coefficients. Figure 4.16 shows the mesh morphing results.Figure 4.15: Streamline tracking algorithm result at Re = 2100, Ma = 0.5 and α = 2◦Figure 4.16: First adapted mesh (black) and morphed mesh (red) before second adaptation504.2. Laminar Flow Test Cases(a) Solution on initial mesh (3871 vertices)(b) Solution after first adaptation (6038 vertices)(c) Solution after second adaptation (8584 vertices)(d) Solution after third adaptation (10667 vertices)(e) Solution on final mesh (13849 vertices)Figure 4.17: Adapted meshes and flow solutions with wake centerline at Re = 2100,Ma = 0.5and α = 2◦ 514.2. Laminar Flow Test CasesThe meshes and flow solutions are depicted in Figure 4.17 and show improvement in theresolution of the boundary layer and wake after each adaptation. Figure 4.18 displays thewake centerline position on the initial and final adapted mesh.Figure 4.18: Wake centerline position - initial (black) and final (red) at Re = 2100, Ma = 0.5and α = 2◦- Closeup view at the trailing edgeNext, an anisotropic quad dominant mesh was generated using the advancing layer meshgenerator. The initial geometry is depicted in Figure 4.19 and consists of the NACA 0012airfoil and the discretized wake centerline obtained from an initial solution.Figure 4.19: Initial geometry for quad dominant mesh at Re = 2100, Ma = 0.5 and α = 2◦(a) Solution on final adapted triangular mesh (13849 vertices)(b) Solution on advancing layer quad dominant mesh (5024 vertices)Figure 4.20: Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 2100, Ma = 0.5 and α = 2◦ , showing better solution behavior in the wake for theadvancing layer mesh524.2. Laminar Flow Test CasesThe flow solutions on the final adapted triangular and advancing layer quad dominant meshesare depicted in Figure 4.20, indicating that the advancing layer mesh can track the wake forlonger distance with many fewer degrees of freedom.The variation in aerodynamic coefficients for each mesh type is displayed in Figure 4.21.For triangular meshes, the convergence of lift and drag coefficients is effectively identicaland the coefficients converge to a single value with and without the wake centerline. Usingthe advancing layer meshes, the aerodynamic coefficients converge faster. For example,the converged lift coefficient obtained from the advancing layer meshes is the same as thegrid converged value of lift coefficient from triangular meshes, although the advancing layermeshes have many fewer degrees of freedom. For the drag coefficient, the difference betweenthe grid converged value from triangular meshes and the advancing layer quad dominantmeshes is just 0.52%.Figure 4.22 shows the convergence of aerodynamic coefficients with respect to total CPUtime for meshing and flow solution. In the case of triangular meshes, there is a CPU timeadvantage for the mesh with wake centerline compared to the mesh without wake centerlinebut just for the last cycle of mesh adaptation. However, there is clearly a CPU time advantagein using the advancing layer mesh. This is because we increased the mesh density in thedesired regions and so with much fewer degrees of freedom overall, we can obtain almostthe same values for aerodynamic coefficients as the grid converged values from triangularmeshes.534.2. Laminar Flow Test Cases(a) Lift coefficient(b) Drag coefficientFigure 4.21: Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 2100, Ma = 0.5 and α = 2◦544.2. Laminar Flow Test Cases(a) Lift coefficient(b) Drag coefficientFigure 4.22: Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 2100, Ma = 0.5 and α = 2◦554.2. Laminar Flow Test Cases(a) Two chords downstream of the trailing edge (b) Five chords downstream of the trailing edge(c) Eight chords downstream of the trailing edgeFigure 4.23: Velocity profiles in the wake for triangular and advancing layer quad dominantmeshesIn Figure 4.23, we compare the velocity profiles in the wake from triangular meshes withwake centerline and the advancing layer meshes at different positions downstream from thetrailing edge. As expected, there is a velocity deficit in the wake. With many fewer degreesof freedom, the velocity profiles from the advancing layer meshes are better for all three po-564.2. Laminar Flow Test Casessitions downstream the trailing edge. Also, the velocity profiles from advancing layer mesheswith 5024 vertices and 7096 vertices are identical. However, with each mesh adaptation,the velocity profiles from the triangular meshes get closer to matching the advancing layerquad dominant meshes. The advancing layer meshes have much fewer degrees of freedomcompared to the fourth adapted mesh. Once again the velocity profiles from the triangularmeshes get worse in comparison as we move further downstream from the trailing edge. Thistest case further confirms that the quad dominant advancing layer meshes can track thewake more efficiently.4.2.3 Viscous flow at Re = 5000, Ma = 0.5 and α = 0◦The third adaptation case is Re = 5000, Mach number Ma = 0.5 and angle of attack α = 0°.The initial mesh used was an isotropic unstructured mesh around NACA 0012 airfoil with6654 cells. At these flow conditions, we get a steady separation in the wake but due topoor resolution of the initial mesh, the flow separation was not observed until the firstadaptation. The wake centerline for the separated flow was tracked using three steps: (1)backward integration in time from the trailing edge, (2) find the saddle point, and (3) add ahorizontal offset from saddle point to escape separation bubble and do forward integrationin time. The flow separation observed is depicted in Figure 4.24a and Figure 4.24b showsthe streamline obtained from streamline tracking algorithm. The streamline, in this case, issplit into two curves and these curves were joined together to give the single curve depictedin Figure 4.24c so that the adapted mesh (which already has the wake centerline insertedinto it) can be morphed to match the new wake centerline. The process was repeated foreach cycle of mesh adaptation.The position of the wake centerline on the initial and final mesh is depicted in Figure4.25. The comparison shows that the wake centerline position from the final adapted mesh ismore accurate since it coincides with the symmetry axis for this symmetric flow. The meshmorphing result for this case is shown in Figure 4.26.574.2. Laminar Flow Test Cases(a) Paraview streamtracer result showing flow separation(b) Streamline tracking algorithm result(c) Streamline after joining two curves in Figure 4.24bFigure 4.24: Wake centerline for separated flowFigure 4.25: Wake centerline position - initial (black) and final (red) at Re = 5000, Ma = 0.5and α = 0◦- Closeup view at the trailing edgeFigure 4.26: First adapted mesh (black) and morphed mesh (red) before second adaptationFigure 4.27 shows the solution on each mesh for this case and shows improvement in theresolution of the boundary layer and wake with each mesh refinement.584.2. Laminar Flow Test Cases(a) Solution on initial mesh (3543 vertices)(b) Solution after first adaptation (5515 vertices)(c) Solution after second adaptation (8851 vertices)(d) Solution on final mesh (12903 vertices)Figure 4.27: Adapted meshes and flow solutions with wake centerline at Re = 5000,Ma = 0.5and α = 0◦594.2. Laminar Flow Test CasesMesh nDOF CL CDPresent resultsInitial 3543 0.003618 0.0642221st adapted 5515 0.000360 0.0568522nd adapted 8851 0.000173 0.0558263rd adapted 12903 0.000155 0.055724Advancing layer 5006 0.000050 0.056000Advancing layer 6766 0.000040 0.056000Zuniga Vazquez [55]Initial 2543 0.014490 0.0611801st adapted 4235 0.016440 0.0348442nd adapted 6688 0.001424 0.0549503rd adapted 10335 0.004885 0.054430Sharbatdar [49]Initial 2543 0.015074 0.9034031st adapted 3899 0.009244 0.763832nd adapted 9140 0.000203 0.0565553rd adapted 19964 0.000108 0.0552304thorder Spectral Volume [22] 34560 - 0.054672ARC2D (320× 128) [16] 40960 - 0.054200R. C. Swanson and S. Langer (4096× 2048) [51] 8388608 - 0.055649Table 4.2: Comparison of lift and drag coefficients for NACA 0012 airfoil at Re = 5000,Ma = 0.5 and α = 0◦Table 4.2 summarizes some of the published results for this test case. The lift coefficient in ourscheme converges faster. For example, the second adapted mesh with wake centerline gives asmaller error in lift coefficient compared to the second adapted mesh of Sharbatdar [49] andthird adapted mesh results from Zuniga Vazquez [55], both of which have more vertices. Thedrag coefficient on the final meshes is similar to these published results. However, our resultof drag coefficient on the third adapted mesh is closer to results obtained on an extremelyfine mesh by R. C. Swanson and S. Langer of NASA Langley Research Center [51].Once again the advancing layer mesh was generated for this case based on the initialgeometry shown in Figure 4.28. The flow solutions on final adapted triangular and advancinglayer quad dominant meshes are depicted in Figure 4.29, showing that the advancing layermesh can track the wake for longer distance with many fewer degrees of freedom.Figure 4.28: Initial geometry for quad dominant mesh at Re = 5000, Ma = 0.5 and α = 0◦604.2. Laminar Flow Test Cases(a) Solution on final adapted triangular mesh (12903 vertices)(b) Solution on advancing layer quad dominant mesh (5006 vertices)Figure 4.29: Flow solutions on final adapted triangular mesh and advancing layer mesh atRe = 5000, Ma = 0.5 and α = 0◦ , showing better solution behavior in the wake for theadvancing layer meshThe plots for the convergence history of aerodynamic coefficients against degrees of freedomare presented in Figure 4.30. The lift coefficient converges faster on advancing layer meshesand has smaller error compared to both triangular meshes since we expect to get CL = 0 forthis symmetric flow. However, comparing the triangular meshes only, the lift coefficient onmesh with wake centerline converges faster as it better symmetry compared to mesh withoutwake centerline. The convergence of the drag coefficient is roughly similar for triangularmeshes but faster on the quad dominant advancing layer meshes.Figure 4.31 shows the convergence history with respect to total CPU time for meshingand solver. For triangular meshes, the total CPU time for meshes with wake centerlineis less compared to the meshes without the wake centerline. Nevertheless, the total CPUtime for the advancing layer quad dominant meshes is much less compared to both thetriangular meshes with advancing layer meshes giving a better result for lift coefficient androughly same drag coefficient as grid converged value on final adapted triangular mesh.The difference between the grid converged value of drag coefficient from adapted triangular614.2. Laminar Flow Test Casesmeshes and the advancing layer quad dominant meshes is only 0.49%.(a) Lift coefficient(b) Drag coefficientFigure 4.30: Convergence of lift and drag coefficients against degrees of freedom (vertices)at Re = 5000, Ma = 0.5 and α = 0◦624.2. Laminar Flow Test Cases(a) Lift coefficient(b) Drag coefficientFigure 4.31: Convergence of lift and drag coefficients against total CPU time for meshingand solver at Re = 5000, Ma = 0.5 and α = 0◦634.2. Laminar Flow Test Cases(a) Two chords downstream of the trailing edge (b) Five chords downstream of the trailing edge(c) Eight chords downstream of the trailing edgeFigure 4.32: Velocity profiles in the wake for triangular and advancing layer quad dominantmeshesThe comparison of velocity profiles in the wake between the advancing layer meshes and theadapted triangular meshes for this case is given in Figure 4.32. Once again, the advancinglayer quad dominant meshes give the expected symmetrical velocity profile in the wake atdifferent positions. The velocity profiles from advancing layer meshes with 5006 and 6766644.3. Turbulent Subsonic Flowvertices are identical at two and five chords downstream but slightly smoother for the morerefined advancing layer mesh at eight chords downstream from the trailing edge. The velocityprofile from the triangular meshes improves with each cycle of mesh adaptation but in thisseparated flow case, the adapted triangular meshes represent the wake velocity profile morepoorly compared to previous test cases especially for five and eight chords downstream fromthe trailing edge where the wake velocity profile is poorly resolved even after the thirdadaptation. This comparison indicates that the advancing layer mesh tracks the wake betterfor separated flow as well.4.3 Turbulent Subsonic FlowFinally, we considered a case for turbulent subsonic flow. The free stream Mach number usedwas Ma∞ = 0.15 while the Reynolds number based on the airfoil chord was Rec = 6× 106.This test case is among of the validation case on NASA Turbulence Modeling Resource(TMR) website [48], which documents results from well known CFD codes. An unstructuredhybrid mesh shown in Figure 4.34 was used. No-slip adiabatic boundary conditions wereapplied on the airfoil surface while the far-field boundary conditions (outflow/inflow basedon the direction of velocity vector) were used at the outer boundary located 500c away fromthe airfoil. We also took into account the effect of curvature and used higher-order cubiccells for a more accurate representation of the boundary faces. However, streamline trackingin curved cells requires using the finite element mapping from straight cells to curved cells.We used the finite element mapping function implemented by Jalali [20] in ANSLib [40] andsolved the non-linear system of equations using Newton’s method to find the curved controlvolume containing the massless particle based on its position at a given time.The system of equations solved using Newton’s method to find the particle in curved cellwas:x = ∑xiφi (ξ, η)y = ∑ yiφi (ξ, η) (4.4)654.3. Turbulent Subsonic Flowwhere φi is the basis function at node i. We then solved for ξ and η such that:x (ξ, η) = x∗y (ξ, η) = y∗The procedure is as follow:• First, find the straight cell in the straight mesh that contains the particle. This givesus a good initial guess.• Then grab the neighbors around that cell and search them one by one. For each cell,solve for ξ and η in the reference element and ensure that ξ and η values are suchthat the point falls inside the triangle. Using the values of ξ and η, re-calculate thevalues of x and y. If x and y are the same as what we started with, then we are in thecorrect cell. Otherwise, we go to other neighbors and repeat the process until we findthe desired cell.(0,0) (1,0)(0,1)Figure 4.33: Mapping from reference triangle to cubic triangular cellWe computed the flow solution at different angles of attack, that is α = 0°, 4° and 10°,using the second-order accurate cell-centered finite volume solver and the velocity field wasintegrated to track the wake centerline in each case. Figure 4.35 shows the streamline trackingresults. These results demonstrate the robustness of the streamline tracking algorithm andthe ability to track the streamlines through mixed element meshes and high-order cubic664.3. Turbulent Subsonic Flowcurved cells. The flow, in this case, is attached and no flow separation was observed fromα = 0° to 10°.Figure 4.34: Unstructured hybrid mesh (n.DoF = 25, 088)Figure 4.35: Wake centerlines at Re = 6 × 106, Ma∞ = 0.15 and α = 0°, 4°, 10° using thestreamline-tracking algorithm67Chapter 5Concluding Remarks5.1 SummaryA new approach for mesh adaptation for wakes has been presented in this research in aneffort to improve the solution accuracy and resolve the wake efficiently. The first step ofthis approach involved the development of streamline tracking algorithm in order to trackthe wake centerline. The velocity data was obtained from solution reconstruction using thesecond-order accurate vertex-centered finite volume solver and integrated numerically usingthe adaptive fourth-order Runge-Kutta method to obtain the wake centerline. The obtainedstreamline was then successfully inserted into the existing isotropic unstructured mesh usingthe surface insertion algorithm in GRUMMP. The surface insertion algorithm takes an initialunstructured mesh and geometry of wake centerline as two initial inputs and its main pointis to recover the wake surface while not changing the mesh size at all more than a cell sizeor two from the wake surface.In the next step, an existing anisotropic mesh adaptation scheme developed by Pagnutti[42] and further extended by Sharbatdar [49] and Zuniga Vazquez [56] was used. The metricbased on solution approximation error was computed at each mesh vertex. The metricrepresents the desired anisotropy for each vertex and is a 2 × 2 symmetric positive definitematrix. The desired anisotropy was then communicated through the metric to the meshadaptation code in GRUMMP to generate an anisotropic mesh. The mesh adaptation schemeproduced highly anisotropic and quasi-structured cells where needed based on derivatives offlow solution variables. Unlike isotropic mesh adaptation which only provides mesh resolutionin desired areas, anisotropic mesh adaptation also provides proper cell alignment with the685.2. Conclusionssolution in addition to mesh resolution to capture anisotropic flow features. The high aspectratio of anisotropic cells and regular structure of quasi-structured cells allows to capture theflow information much more efficiently compared to isotropic cells. Refining the mesh insuch a way based on solution features helps to produce an optimal mesh for computing thesolution. The mesh adaptation scheme produces high quality meshes using the four principaloperations for mesh quality improvement as discussed in Section 2.5. The mesh modificationoperations preserve the wake centerline and do not move the vertices on it at all. Afterobtaining the adapted mesh, the solution was re-computed and a new wake centerline wasobtained.In the second and subsequent mesh adaptation cycles, the mesh was morphed and meshpoints from old wake centerline were projected to match the new wake centerline by solvingthe linear elasticity problem, using displacement between two wake centerlines as constraints.This helped to move the wake centerline to a more accurate position and preserved thealready adapted regions, thus making the adaptation process more efficient.An advancing layer mesh generator was used as an alternative technique to produce ananisotropic quad dominant mesh and the results are compared with triangular meshes withand without the wake centerline. We also presented the comparison of velocity profiles inthe wake at different locations downstream from the trailing edge for the advancing layerquad dominant and adapted triangular meshes.5.2 ConclusionsTo verify the correctness of streamline tracking, results are presented for an analytical veloc-ity field which confirmed the order of accuracy of the numerical method and the ability totrack streamline through a mesh using the reconstructed velocity data. The accuracy of themesh morphing approach is verified by using the method of manufactured solutions whichproved that the linear elasticity solver is second-order accurate.Test cases for viscous flow at different flow conditions have been considered and the695.2. Conclusionsrobustness of the mesh adaptation scheme is demonstrated by improvement in solution ac-curacy and resolution of the boundary layer and wake as the adaptation process is repeated.The flow separation not observed on the initial triangular mesh in the third test case high-lights the fact that it is immensely important to have an appropriate resolution of flowsolution in the boundary layer and wake. The wake centerline for the first and third testcases should coincide with the symmetry axis; comparison on initial and final adapted tri-angular meshes show that wake centerline on the final mesh is closer to the symmetry line.Comparing just the triangular meshes, the convergence of aerodynamic coefficients demon-strated that there is a small CPU time advantage in case of meshes with wake centerlineand solution converges slightly faster especially for the third test case. We expect to get liftcoefficient equal to zero at zero angle of attack but this is not the case due to asymmetry inthe mesh. The triangular meshes with wake centerline have better symmetry compared tomeshes without wake centerline which plays an important role in faster convergence of thesolution.Nevertheless, the comparison of results between the triangular meshes and the advancinglayer quad dominant meshes showed that the aerodynamic coefficients converge much fasteron the advancing layer meshes in all the test cases. There is a smaller error in lift coeffi-cient on the advancing layer meshes for symmetric flow cases compared to both triangularmeshes since the advancing layer quad dominant meshes have better symmetry comparedto triangular meshes. Also, there is a clearly noticeable CPU time advantage in using theadvancing layer meshes. This is because it is much cheaper to generate the advancing layermesh rather than doing mesh adaptation which saves CPU time required for meshing. Theadvancing layer mesh can also track the wake better with many fewer degrees of freedomwhich helps in reducing the CPU time for the solver.The velocity profiles from the advancing layer meshes provide evidence that it can capturethe wake more efficiently even up to several chords downstream from the trailing edge. Thevelocity profile on the initial triangular meshes is not physical in all the cases due to poormesh resolution in the wake. The velocity profiles from triangular meshes improve with mesh705.3. Recommended Future Workadaptation and approach the velocity profiles represented by the quad dominant advancinglayer meshes. The final adapted triangular meshes already have almost twice the degreesof freedom compared to the quad dominant advancing layer meshes. However, the velocityprofiles from triangular meshes get worse in comparison as we move further downstreamfrom the trailing edge such as at five and eight chords especially for the third test case withflow separation. This comparison study illustrates that the advancing layer mesh is moreefficient in resolving the wake.Moreover, the robustness of the streamline tracking algorithm is further demonstratedby tracking the wake centerline at different angles of attack for a fully turbulent flow. Theturbulent flow case demonstrates our ability to track streamlines through mixed elementmeshes and high-order cubic curved cells.5.3 Recommended Future WorkThis thesis demonstrated a new approach for anisotropic mesh adaptation for wakes. Thereare several areas to which the overall methodology can be extended and improved further.• The streamline tracking algorithm tracks the streamlines for 2-D flows. One possibleaddition in the long run can be to extend the algorithm to 3-D flows and trace a streamsurface.• In the future, it will be good to do the anisotropic mesh adaptation for turbulent flowcase as well and follow the mesh adaptation loop discussed in Section 3.5. It will alsobe good to solve the drag polar for turbulent flow.• For quad dominant advancing layer mesh, we produced the mesh by using an average ofthe laminar boundary layer and laminar far-wake profiles. One potential improvementfor this can be to explore ways to get a better estimate of the overall profile that canresolve the wake even better especially in the near-wake region.• In addition to resolving the boundary layer and wake, we also need to resolve the shock715.3. Recommended Future Workif the flow is transonic. Considering this additional feature of transonic flow, it willbe good to identify the shock surface and insert it into the mesh along with the wakecenterline and do anisotropic mesh adaptation. Also, then the initial geometry for theadvancing layer approach can be modified to include both the wake centerline and theshock and march off both the surfaces to generate the advancing layer mesh. Thismight help to resolve the shock more efficiently in addition to the wake.72Bibliography[1] National Aeronautics and Space Administration (NASA). Definition of stream-lines. https://www.grc.nasa.gov/www/k-12/airplane/stream.html, May 2015. Accessed:2017-12-30.[2] Zaib Ali, Paul G Tucker, and Shahrokh Shahpar. Optimal mesh topology generation forCFD. Computer Methods in Applied Mechanics and Engineering, 317:431–457, 2017.[3] Michael Van Altena. High-order finite-volume discretisations for solving a modifiedadvection-diffusion problem on unstructured triangular meshes. Master’s thesis, TheUniversity of British Columbia, Department of Mechanical Engineering, 1999.[4] Thomas Apel. Anisotropic Finite Elements: Local Estimates and Applications. 1999.[5] Stanley A Berger. Laminar Wakes. American Elsevier Publishing Company, Inc., NewYork, 1971.[6] Carlo L Bottasso. Anisotropic mesh adaption by metric-driven optimization. Interna-tional Journal for Numerical Methods in Engineering, 60(3):597–639, 2004.[7] Michael L Brewer, Lori Freitag Diachin, Patrick M Knupp, Thomas Leurent, and Dar-ryl J Melander. The Mesquite Mesh Quality Improvement Toolkit. In Proceedings ofthe 12th International Meshing Roundtable, pages 239-250, 2003.[8] Zhiqiang Cai. On the finite volume element method. Numerische Mathematik,58(1):713–735, 1990.73Bibliography[9] Anil W Date. Introduction to Computational Fluid Dynamics. Cambridge UniversityPress, 2005.[10] D Davidson. The role of computational fluid dynamics in process industries. In Frontiersof Engineering: Reports on Leading-Edge Engineering from the 2002 NAE Symposiumon Frontiers of Engineering, page 21. National Academies Press, 2003.[11] MN Dhaubhadel. CFD applications in the automotive industry. Journal of FluidsEngineering, 118(4):647–653, 1996.[12] Cécile Dobrzynski and Pascal Frey. Anisotropic delaunay mesh adaptation for unsteadysimulations. In Proceedings of the 17th International Meshing Roundtable, pages 177–194. Springer, 2008.[13] Lori A Freitag and Carl Ollivier-Gooch. The effect of mesh quality on solution efficiency.In Proceedings of the Sixth International Meshing Roundtable, page 249, 1997.[14] Pascal-Jean Frey and Frédéric Alauzet. Anisotropic mesh adaptation for CFD com-putations. Computer Methods in Applied Mechanics and Engineering, 194:5068–5082,2005.[15] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins Uni-versity Press, Baltimore, 3rd edition, 1983.[16] Carl F. Gooch. Solution of the Navier-Stokes Equations on Locally-Refined CartesianMeshes Using State-Vector Splitting. PhD thesis, Stanford University, 1993.[17] D Holmes and S Connell. Solution of the 2D Navier-Stokes equations on unstructuredadaptive grids. In 9th Computational Fluid Dynamics Conference, AIAA Paper 89-1932,1989.[18] Weizhang Huang. Metric tensors for anisotropic mesh generation. Journal of Compu-tational Physics, 204(2):633–665, 2005.74Bibliography[19] Florence V Hutcheson, Daniel Stead, and Gerald Plassman. Experimental studyof wake/flap interaction noise and the reduction of flap side edge noise. In 22ndAIAA/CEAS Aeroacoustics Conference, AIAA paper 2016-2955, 2016.[20] Alireza Jalali. An adaptive higher-order unstructured finite volume solver for turbulentcompressible flows. PhD thesis, The University of British Columbia, Department ofMechanical Engineering, 2017.[21] Alireza Jalali and Carl F Ollivier Gooch. Higher-order finite volume solution recon-struction on highly anisotropic meshes. In 21st AIAA Computational Fluid DynamicsConference, AIAA Paper 2013-2565, 2013.[22] Ravi Kannan and Zhijian Wang. Improving the high order spectral volume formulationusing a diffusion regulator. Communications in Computational Physics, 12(1):247–260,2012.[23] Abhishek Khare, Ashish Singh, and Kishor Nokam. Best practices in grid generationfor CFD applications using hypermesh. Computational Research Laboratories, page 2,2009.[24] P Knupp, L Freitag-Diachin, and B Tidwell. Mesquite: Mesh Quality ImprovementToolkit User’s Guide. Sandia National Laboratories, 2013.[25] Randall J LeVeque. Finite Volume Methods for Hyperbolic Problems (Cambridge Textsin Applied Mathematics). Cambridge University Press, 2002.[26] Zhenquan Li and Gordon Mallinson. An adaptive streamline tracking method for twodimensional CFD velocity fields. In Proceedings of the 7th Asian Symposium on Visu-alization, pages 1–6, Singapore, November 2003.[27] Vladimir D Liseikin. A computational differential geometry approach to grid generation.Springer Science & Business Media, 2006.75Bibliography[28] Harvard Lomax, Thomas H Pulliam, and David W Zingg. Fundamentals of Computa-tional Fluid Dynamics. Springer Science & Business Media, 2013.[29] Adrien Loseille. Metric-orthogonal anisotropic mesh generation. Procedia Engineering,82:403 – 415, 2014. 23rd International Meshing Roundtable (IMR23).[30] Adrien Loseille and Rainald Lohner. Anisotropic adaptive simulations in aerodynam-ics. In 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum andAerospace Exposition, page 169, 2010.[31] Shahzaib Malik and Carl F Ollivier Gooch. Mesh Adaptation for Wakes via SurfaceInsertion. In AIAA SciTech 2019 Forum, AIAA Paper 2019-1996, 2019.[32] DJ Mavriplis. Unstructured grid techniques. Annual Review of Fluid Mechanics,29(1):473–514, 1997.[33] Fotis Mavriplis. CFD in aerospace in the new millennium. Canadian Aeronautics andSpace Journal, 46(4):167–177, 2000.[34] Krzysztof Michalak and Carl Ollivier-Gooch. Matrix-explicit GMRES for a higher-orderaccurate inviscid compressible flow solver. In 18th AIAA Computational Fluid DynamicsConference, AIAA Paper 2007-3943, 2007.[35] F Moukalled, L Mangani, M Darwish, et al. The finite volume method in computationalfluid dynamics - An Advanced Introduction with OpenFOAM and Matlab. Springer,2016.[36] Amir Nejat and Carl Ollivier-Gooch. A high-order accurate unstructured finite volumeNewton–Krylov algorithm for inviscid compressible flows. Journal of ComputationalPhysics, 227(4):2582–2609, 2008.[37] Logan Numerow and Carl F Ollivier Gooch. A mixed-element 3D advancing layer meshgenerator for complex geometries. In 23rd AIAA/CFD Conference, AIAA Paper 2017-3453, 2017.76Bibliography[38] Carl Ollivier-Gooch. An unstructured mesh improvement toolkit with ap-plication to mesh improvement, generation and (de-) refinement. In 36thAIAA Aerospace Sciences Meeting and Exhibition, AIAA Paper 98-0218, 1998.http://tetra.mech.ubc.ca/GRUMMP/.[39] Carl Ollivier-Gooch and Michael Van Altena. A high-order-accurate unstructured meshfinite-volume scheme for the advection–diffusion equation. Journal of ComputationalPhysics, 181(2):729–752, 2002.[40] Carl F. Ollivier-Gooch. ANSLib: A scientific computing toolkit supporting rapid numer-ical solution of practically any PDE. In Proceedings of the Eighth Annual Conference ofthe Computational Fluid Dynamics Society of Canada, pages 21–28. Societe canadiennede CFD / CFD Society of Canada, June 2000.[41] Doug Pagnutti. Anisotropic adaptation: Metrics and meshes. Master’s thesis, TheUniversity of British Columbia, Department of Mechanical Engineering, 2008.[42] Doug Pagnutti and Carl Ollivier-Gooch. A generalized framework for high orderanisotropic mesh adaptation. Computers and Structures, 87(11-12):670–679, 2009.[43] Michael A Park, Joshua A Krakos, Todd Michal, Adrien Loseille, and Juan J Alonso.Unstructured grid adaptation: Status, potential impacts, and recommended investmentstoward CFD Vision 2030. 2016.[44] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery. Nu-merical Recipes 3rd edition: The Art of Scientific Computing. Cambridge UniversityPress, 2007.[45] Shmuel Rippa. Long and thin triangles can be good for linear interpolation. SIAMJournal on Numerical Analysis, 29(1):257–270, 1992.[46] Philip L Roe. Approximate Riemann solvers, parameter vectors, and difference schemes.Journal of Computational Physics, 43(2):357–372, 1981.77Bibliography[47] Rolls-Royce. BR715 for the Boeing 717. https://www.rolls-royce.com/products-and-services/civil-aerospace/airlines/br715.aspx. Accessed: 2018-12-25.[48] C Rumsey. Turbulence Modeling Resource, NASA Langley Research Center, 2018.[49] Mahkame Sharbatdar. Anisotropic mesh adaptation: Recovering quasi-structuredmeshes. Master’s thesis, The University of British Columbia, Department of MechanicalEngineering, 2012.[50] Justin Solomon. Numerical algorithms: methods for computer vision, machine learning,and graphics. CRC Press, 2015.[51] R Charles Swanson and S Langer. Comparison of NACA 0012 Laminar Flow Solutions:Structured and Unstructured Grid Methods. Technical Memorandum NASA-TM-2016-219003, NASA Langley Research Center, 2016.[52] Joe F Thompson, Bharat K Soni, and Nigel P Weatherill. Handbook of Grid Generation.CRC press, 1998.[53] W Tollmien. Handbuch der Experimental Physik. Vol. 4, Hydro- und Aero-Dynamik,Part 1. Akademie-Verlag, Leipzig, 1931.[54] Jiyuan Tu, Guan Heng Yeoh, and Chaoqun Liu. Computational Fluid Dynamics: Apractical approach. Butterworth-Heinemann, 2018.[55] Zuniga Vazquez. Adaptive creation of orthogonal anisotropic triangular meshes usingtarget matrices. Master’s thesis, The University of British Columbia, Department ofMechanical Engineering, 2015.[56] Zuniga Vazquez and Carl F Ollivier Gooch. Target-edge based orthogonal anisotropicmesh. In 22nd AIAA Computational Fluid Dynamics Conference, AIAA Paper 2015-2294, 2015.[57] Henk Kaarle Versteeg and Weeratunge Malalasekera. An Introduction to ComputationalFluid Dynamics: The Finite Volume Method. Pearson Education, 2nd edition, 2007.78[58] Li Wang, W Kyle Anderson, and Jon Erwin. Solutions of high-order methods for three-dimensional compressible viscous flows. In 42nd AIAA Fluid Dynamics Conference andExhibit, AIAA Paper 2012-2836, 2012.[59] Frank M White. Viscous Fluid Flow. McGraw-Hill Book Company, 3rd edition, 2006.[60] Daniel W Zaide and Carl F Ollivier-Gooch. Inserting a surface into an existing unstruc-tured mesh. International Journal for Numerical Methods in Engineering, 106(6):484–500, 2016.79Appendix ACommand-line OptionsThe command line options used to obtain the results in Chapter 4 are included below.A.1 Wake Centerline TrackingThe executable file for wake centerline tracking in ANSLib is streamline-tracking. It canbe executed through the Linux shell with the following command line options:./streamline-trackingOptions-f Address of .mesh file without extension-physics () RoeVisc2D for laminar and RoeTurbSA2D for turbulent flow-mesh_type v Vertex-centered mesh-reynolds () Reynolds number-mach () Mach number-angle () Angle of attack-initial_step Initial step size (for forward and backward integration)-CASE 2 CASE 2 is for actual velocity data from flow solverAlternative Options-initial_step_backward Initial step size for backward integration only (optional)-CASE () Use CASE 1 for analytical velocity field-mesh_type c Cell-centered mesh80A.2. Wake Centerline Insertion into the MeshA.2 Wake Centerline Insertion into the MeshThe executable file in GRUMMP for two-dimensional curve insertion is surfins2d. Thefollowing command line options should be used:./surfins2dOptions-i Address of .mesh file with extension-z Address of .bdry file without extension-I Melts internal boundary faces and gives .mesh.mesh fileA.3 Obtaining the metric fileThe metric file can be obtained from ANSLib using the Driver executable. The commandline options are:./DriverOptions-f Address of .mesh file without extension-physics RoeVisc2D Physics-mesh_type v Vertex-centered mesh-reynolds () Reynolds number-mach () Mach number-angle () Angle of attack-AA Anisotropic adaptation-r 2 Adaptation orderA.4 Anisotropic Mesh Adaptation (Metric-Based)The executable file used for metric-based anisotropic mesh adaptation is tmop_adapt2dand can be found in GRUMMP. The following command line options should be used:81A.5. Mesh Morphing./tmop_adapt2dOptions-i Address of .mesh file or .grmesh file with extension-B Address of .bdry file without extension-I Melts internal boundary faces if .grmesh file is usedAdditional file that must be providedA metric file must be provided for the corresponding mesh file and should have same extension asmesh file with additional .metric at the end. For example, if the mesh file is 0012.mesh, then themetric file should be 0012.mesh.metric.NoteFor .grmesh, the boundary file should have coordinates of the curve (e.g wake centerline) insertedinto the mesh.A.5 Mesh MorphingThe executable file used for mesh morphing is LinearElements and can be found in ANSLib.The following command line options should be used:./LinearElementsOptions for general cases-i Address of .mesh file without extension-B Address of .bdry file without extension (Initial wake centerline)-C Address of .bdry file without extension (New wake centerline)Options to test mesh morphing code-i Address of .mesh file without extension-MMS Method of manufactured solutionsA.6 Advancing Layer Mesh GenerationThe executable file used for advancing layer mesh generation is edam2d and can be foundin GRUMMP. The following command line options should be used:82./edam2dOptions-i Address of .bdry file without extension-o