Latest revision as of 11:58, 8 October 2020

Abstract

This thesis develops a numerical tool (the Virtual Wind Tunnel, VWT) for the resolution of problems involving fluid flow around structures. Due to the limitations that traditional methods may have in this context, the VWT is based on the use of fixed mesh technologies (CutFEM-type) combined with an implicit representation of the embedded bodies.

One of the main contributions of the thesis is the use of such fixed mesh methods to solve lightweight thin-walled structures problems. Hence, two embedded formulations capable of representing the flow around bodies with and without internal volume are proposed. The first one results in a simpler implementation and lower computational effort but can only represent a slip behavior of the wall. The second one gets rid of such limitation by including a Nitsche imposition of the Navier-slip condition, thus allowing modelling any wall behavior as a wall law would do.

The applicability range of the VWT includes the fluid–structure interaction problem (FSI). To that purpose an improvement for the boundary condition imposition of the FM-ALE algorithm mesh motion problem is also proposed. Moreover, the implementation, which has been conceived to be easily extended to any other coupled problem, is also treated.

The validation of the technologies within the VWT includes multiple theoretical test cases as well as feasible industrial applications. Among these, the FSI analysis of a 4-point tent during a strong wind episode deserves to be highlighted as it showcases the achievement of the initial objective of the thesis.

Resumen

spanish En esta tesis se desarrolla una herramienta numérica (el Virtual Wind Tunnel, VWT) para la resolución de problemas que involucran el flujo de un fluido alrededor de una estructura. Debido a las limitaciones que los métodos tradicionales pueden tener en este contexto, el VWT se basa en el empleo de técnicas de malla fija (tipo CutFEM) combinadas con una descripción implícita de los cuerpos embebidos.

Una de las principales contribuciones de la tesis es el empleo de dichos métodos de malla fija para resolver problemas de estructuras ligeras de pared delgada. Así pues, se proponen dos formulaciones embebidas capaces de representar el flujo alrededor de un cuerpo con o sin volumen interno. La primera de ellas resulta en una implementación más sencilla así como en un menor coste computacional pero únicamente puede reprentar un comportamiento deslizante de la pared. La segunda elimina esta limitación incluyendo una imposición mediante el método de Nitsche de la condición de Navier-slip, permitiendo así modelar cualquier comportamiento del mismo modo que lo haría una ley de pared.

El rango de aplicabilidad del VWT incluye el problema de interacción fluido–estructura (FSI). A tal propósito se plantea una mejora para la imposición de las condiciones de contorno del problema de movimiento de la malla del algoritmo FM-ALE. Asimismo, también se hace especial hincapié en la implementación, que ha sido concebida para ser fácilmente extensible a cualquier otro problema acoplado.

La validación de las tecnologías implementadas en el VWT incluyen múltiples casos teóricos así como posibles aplicaciones industriales. Entre éstas se destaca el análisis FSI de una 4-point tent durante un episodio de viento severo ya que demuestra la consecución del objetivo inicial de la tesis.

Acknowledgements

First of all, I would like to acknowledge the Spanish Minister of Science, Innovation and Universities (Ministerio de Ciencia, Innovación y Universidades) and the International Center for Numerical Methods in Engineering (CIMNE) for the financial support during these four years. In this regard, the International Graduate School of Science and Engineering (IGSSE) of the Technische Universität München (TUM) is also acknowledged. Furthermore, I also thank the Oskar von Miller Forum (OvMF) for the accommodation during my international stay.

My gratitude extends to my two directors. I thank Prof. Eugenio Oñate for allowing me to work in such an interesting topic and for the confidence in me to achieve the objectives of this thesis. This work would not have been completed without the help of Prof. Rossi. Riccardo, many thanks for always having time to invest on this. I hope this is just the very fist one of many upcoming research successes.

The entire Kratos team deserve to be acknowledged by its own merits. Among all the developers, I particularly thank Dr. Jordi Cotela, Dr. Eduardo Soudah, Dr. Antonia Larese, Dr. Pooyan Dadvand, Carlos A. Roig (aka Charlie), Philipp Bucher and (soon Dr.) Vicente Mataix. Jordi, thanks for the valuable discussions (possibly lessons) about fluid dynamics and turbulence as well as for the coding advises. I also thank you for the fun during my stay in Munich. Edu, thanks for giving me the chance of applying my research in the biomedical field and for constantly try to find new opportunities for the group (by the way, you still owe me some hours of conversation). Antonia, thanks for all the help since I became your student for the first time. I want to also thank the KratosCore-man Pooyan for the enormous help while developing the level set algorithms as well as for the C++ and design patterns guidelines. The guy who silently makes everything work must be acknowledged too. Charlie, thanks for all the help with the compilation and debugging weird issues. I also thank the guy who has most probably reviewed the entire Kratos many times to force all of us to do things properly. Philipp, thanks for all the code reviews and brilliant suggestions. At this point, I must also thank the coolest office mate I ever had. Vicente, thanks for allowing me to constantly tease you, for the laughs during the tough moments and for being always disposed to help. I want to also thank the rest of CIMNE staff, specially the GiD team and the rookies that brought fresh air to the office a few years ago.

My gratitude also goes for the Bavarian colleagues in the Statik Chair. I personally thank Prof. Roland Wüchner for all the help and the guys of the Raum 0101.Z1.026, which I can confirm is the loudest office in the Chair. Thanks Andreas, Altug, Íñigo and Tobi for making me feel like at home during my stay.

Furthermore, I would like to express my gratitude to Prof. Miguel Cervera for the confidence in me to lecture one of his subjects. In this regard, I also acknowledge my teaching colleague Alejandro Cornejo as well as my colleague in the Department Gabriel Barbat.

There are people that has contributed to make this thesis possible without realizing. Here I want to thank my friends Toni (aka Kaborza), Ana, Enric (aka Epi), Eva María (aka Kpo), Dani Rey, Laura Cuni and Carlos (aka Carlutus). I also thank my long term teammates Víctor, Ortuño, Samu (aka Chino), Fer and Javi, as well as the last seasons joins Nachete, Nacho (aka Bibi), Andreu, Antonio and Lluís. These lists lack the important former teammates Álex, Adri, Joaquín (aka Maka), Bistu, Aitor and Manolo. All of them made me (temporally) forget about the troubles that raised during the prosecution of this thesis.

Last but not least, I want to also thank my family. Gracias Elena por acompañarme durante estos años y conseguir siempre hacerme reír, especialmente en los momentos más duros. Agradecer a mi madre María Luisa Martínez por su amor incondicional y por constantemente cuidar de nosotros. Gustaríame tamén dedicar unas palabras a o recordo da miña avoa Hortensia. Finalmente, mi más sentido agradecimiento es para mi padre Juan Francisco Zorrilla. Gracias Papa por enseñarme como luchar frente a las adversidades, por siempre confiar ciegamente en mí y por constantemente empujarme a dar un paso más allá. Esta tesis es una consecuencia de todo lo que he aprendido de ti.

input-ch2-StateArt input-ch3-NavierStokes input-ch4-CMAME_emb_1 input-ch5-AusasNavierSlip input-ch6-BIO input-ch7-CMAME_fsi input-ch8-Conclusions

1 Introduction

1.1 Fluid dynamics. A brief historical survey

Since the ancient era, humans have tried to tame the fluids for their own benefit. In the very beginning, this was limited to take profit of the wind thrust for navigation purposes. The earliest evidence of the use of a sail appears on a painted ceramic disk dating between 5500 and 5000 BC which was found in the current Kuwait. This discovery suggests that primitive wind sailing was already known in the Ubaid 3 period of the Mesopotamian empire [1]. Later on, the watermill appeared during the Hellenistic period (early 3rd century BC) in Byzantium and Alexandria to harness power from water [2,3].

These technologies kept evolving during centuries. Hence, the advances in sail navigation made possible the goods commerce as well as the culture exchange among all the civilizations in the Mediterranean sea and, after the 15th century, around the world. Watermills also evolved, to be more efficient or to become windmills, which appeared in the current Iran in the mid 9th century to generate power from the wind [4]. Although these are only two examples it is clear that these and other fluid-related technologies helped to write the history of mankind.

Furthermore, the technological advances nourished the interest on the theoretical comprehension of the mechanics of fluids. In this context, the Archimedes' principle (250 BC) can be considered as the very first theoretical milestone, which became the base for the hydrostatics studies of the ancient Greek philosophers. A few centuries later, the middle age Islamic physicists took the lead on this. This allowed them to reach a deep comprehension about the hydraulics of many engineering structures.

The breakthrough into modern fluid dynamics came in the 16th century from the hand of Leonardo da Vinci (1452-1519). Based on the observation of nature, da Vinci achieved an empirical understanding of some of the most complex fluid dynamics phenomena. Some of these are the first description of a pressure lift force, which can be understood as the first theory of flight, or the first explanation of turbulence, which he did from the observation of the water flow patterns that generate behind obstacles in rivers.

Nonetheless, the 17th and 18th centuries can be considered as the most prolific times for the theory of fluid dynamics thanks to the contributions of many brilliant mathematicians. Among all of them, Isaac Newton (1643-1727), Daniel Bernoulli (1700-1782), Leonhard Euler (1707-1783) and Jean le Rond d'Alembert (1717-1783) stand out. Isaac Newton settled the foundations with his well-known laws of motion and studies about the friction and viscosity (Principia, 1687). Later on, Daniel Bernoulli proved that the gradient of pressure is proportional to the acceleration of the fluid in his theory of the motion of fluids (Hydrodynamica, 1738). In 1757 Leonhard Euler posed his famous Euler differential equations for inviscid flows. Just one year later, Jean le Rond d'Alembert proved that the drag of a body of any shape moving through a fluid with no viscosity is zero (1758). This phenomenon was known from then on as the d'Alembert paradox.

During the 19th century the effort was put on adding a friction term to Euler's equations in order to overcome the d'Alembert paradox and obtain the results observed in nature. This was achieved in 1822 by Claude-Louis Navier (1785-1836), who proposed the addition of an extra molecular viscosity term. Although several authors published works in this regard, it took two decades to find the first mathematical reversible derivation of the new viscous Euler equations, which was presented by George Gabriel Stokes (1819-1903) in 1845. Hereafter, these equations were known as the Navier-Stokes equations. In 1851 George Gabriel Stokes also derived the well-known Stokes' law for the calculation of the drag force exerted on spherical objects in high viscous flows.

Viscosity and turbulence have been possibly the major theoretical concerns of contemporary fluid dynamics. In 1883 Osborne Reynolds (1842-1912) demonstrated with his famous experiment how the transition from laminar to turbulent flow occurs. 20 years later, Ludwig Prandtl (1875-1953) presented his viscous boundary layer theory (1904) which has been the basis for most of the aerospace technological advances of the past century [5]. Besides these, the contributions of Theodore von Karman (1881-1963), Geoffrey Ingram Taylor (1886-1975) and Andrey Nikolaevich Kolmogorov (1903-1987) aiming at understanding the turbulence also deserve to be mentioned.

The eruption of the World Ward in the beginning of the 20th century came with an astonishing evolution of the aerospace engineering, which promoted the construction of experimental wind tunnel facilities. The wind tunnels proliferation continued until the end of the World War II and allowed scientists and engineers to properly visualize and accurately measure the flow passing through objects for first time in history. This capability turned into a deep empirical knowledge that made possible the amazing aerodynamic advances of the past century.

Although numerical methods were used from long time ago, it was the apparition of modern computers during the second half of the 20th century what effectively enabled the virtual resolution of real-life problems by using computational techniques. Hence, the constantly growing computational power paced the apparition of a widespread range of numerical techniques since the 50s. It was however during the 60s and the 70s when the aircraft industry, mainly Boeing and NASA, effectively drove the Computational Fluid Dynamics (CFD) forward by developing numerical techniques to solve simplified inviscid models, namely the potential and Euler equations. Thanks to this, CFD methods could be applied for the first time to the resolution of real engineering problems. During the 80s and the 90s CFD technologies advanced further on and achieved the resolution of the complete Navier-Stokes equations, helping CFD methods to populate other engineering fields such as the automotive, biomedical or civil engineering.

In addition, during the early 80s commercial CFD codes came into the open market to become one of the daily used Computer-Aided Engineering (CAE) tools. Nowadays, commercial and open-source CFD codes provide technicians, designers and researchers with a virtual wind tunnel toolbox on their own laptops. Thus, CFD have become in an indispensable part of any design process involving fluids, making possible the development of technologies that have helped to achieve the all-times most rapid mankind evolution.

1.2 The fluid–structure interaction problem

From the previous section, it can be concluded that the comprehension of the interaction mechanisms that appear after the immersion of a body into a fluid has been one of the most important scientific challenges of the past centuries. In a wide range of applications, such as the hydrostatic analysis of massive civil engineering structures or the laminar flow analysis in open water channels, the fluid dynamics analysis can be efficiently achieved without considering the deformation of the “wet” bodies. However, there are lots of engineering and natural phenomena in which such simplification is no longer valid, leading to the so called Fluid–Structure Interaction (FSI) problems.

At this point, it is important to remark that the resolution of a fluid dynamics problem with and without considering such fluid–structure interaction mechanisms can lead to a completely different response. Hence, the FSI problem can be somehow considered as the union of a fluid dynamics and a solid mechanics problems such that the solution of each one of them is linked to the one of the other. This mutual dependency comes from the fact that the structure is deformable under the action of the fluid load. At the same time, the structure deformation implies a modification in the fluid problem geometry, which turns into a change in the load distribution over the structure.

Nowadays this might seem kind of obvious, but uncontrolled FSI is behind some of the most famed engineering failures of the modern era. For instance, fluttering, which is a dynamic instability that causes self-oscillation of fluid-surrounded slender bodies, is understood to be the main cause of the Tacoma Narrows bridge collapse in 1940 as well as the Braniff International Airways flight 542 crash in 1959.

Furthermore, there also exists lots of natural phenomena, such as the flapping of animal wings or the hydrodynamics of fish fins that are indeed FSI problems [6,7,8]. The study of such animal motion mechanisms from a computational perspective [9,10] is of great value for the development of novel technologies such as newfangled unmanned vehicles. Such natural phenomena also include many of the human body processes, mainly those related to the circulatory and respiratory systems. The numerical simulation of these has helped the clinicians in the decision making in front of cardiovascular [11,12,13] and airways [14] diseases, in the development of new prosthetic devices such as artificial heart valves [15,16,17], and more recently, in the study of the mechanics of blood cell motion [18].

More specifically, there are also lots of civil engineering structures that are affected by the FSI problem. These can be basically summarized in thin-walled structures [19], namely membranes and shells, and slender structures such as bridge decks [20], chimneys, elevated tanks that suffer from sloshing during earthquakes [21,22] or the blades of wind energy harvesting mills [23].

In some cases, the FSI effects can be efficiently reproduced using reduced order models [24] (e.g. water hammer in fluid ducts). However, such simplifications are no longer possible in the vast majority of applications since the interaction mechanisms involve a three-dimensional behavior. This thesis focuses on the latter, particularly in those problems that involve lightweight thin-walled structures.

1.3 Lightweight thin–walled structures

Since their appearance in the mid past century, the interest in the use of highly flexible thin–walled structural solutions has been constantly growing. In this regard, the German architect and engineer Frei Otto (1925-2015) must be acknowledged as the pioneer in the design of lightweight structures. Under the idea of creating more sustainable structures in a world with resources shortages, Frei Otto abandoned the traditional structural design in favour of more efficient and creative solutions based on the use of metal reinforced tensile membranes and cables [25]. Although this new technology was initially used for basic structures, such as the Musikpavillon in the 1955 Kassel Federal Garden exhibition (Fig. 1.1a), their advantages rapidly became apparent and encouraged their use as permanent solutions for the roof covers of the 1967 Montreal World Expo German pavilion and the 1972 Munich Olympic Stadium (Olympiastadion München) (Fig. 1.1b).

(a) Kassel Federal Garden exhibition Musikpavillon (1955).	(b) Olympiastadion München interior view (1972).
Figura 1.1: Evolution of lightweight structures by Frei Otto.

The pioneer understanding of how structures must be conceived that Frei Otto had impinged the current structural design criteria, which have evolved to use of much more efficient materials such as reinforced textile composites. Even though such extremely lightweight structures can be also employed as roof covers or shadowing structures (Fig. 1.2), their versatility and resource-saving manufacturing opens a wide variety of novel applications.

Figura 1.2: Textile roof cover in the airport of München.

For instance, these can be used for the creation of temporary rapid deployment inflatable structures that can serve as emergency shelters, exposition pavilions or hangars (Fig. 1.3) [26]. Furthermore, they can also be used for the creation of foldable or moving structures capable of adapting their shape according to the environmental conditions (Fig. 1.4). Last but not least, there are other problems out of the civil engineering scope, such as the optimization of sail navigation or the design of parachutes, that could directly benefit from the research on this field.

(a) Exterior view.	(b) Interior view.
Figura 1.3: Inflatable 75m span hangar in Saudi Arabia (courtesy of Buildair).

(a) Shading structure closed.	(b) Shading structure opened.
Figura 1.4: Convertible shading roof in the Prophet's Mosque in Medina. In summer the umbrellas open to provide shade while in winter they do so to prevent the day warmth from escaping (courtesy of SL Rash GmbH).

It is clear that the structural safety analysis of this family of structures depends hugely on the characterization of the mechanical properties of the textile materials they are composed of. Their response is also strongly influenced by the wind load they are subjected to.

Normally, the structural integrity under the action of wind loads is ensured by the accomplishment of the design codes requirements. Unfortunately, these newfangled structural solutions are commonly out of the scope of the standards, and thus require alternative ways to guarantee safety and reliability. Moreover, they are prone to suffer from aeroelastic effects, which should definitively be considered in the response assessment, owing to their extremely flexible nature.

Even though wind tunnel facilities have been of great benefit to characterize the flow around flexible objects, they have some inherent drawbacks. On the one hand, wind tunnels involve enormous construction and operation costs that compromise their economic viability. This turns into an increased cost per experiment, which restricts their use to reference research centers or big private companies. On the other hand, depending on the dimensions of the problem and of the ones of the wind tunnel, the experiment may require the use of a reduced scale model. In some applications this might compromise the quality of the results as it is not always possible to properly reproduce the scale effects in the reduced experiment in such a way that the aeroelastic behavior of the real structure is represented. In this context, numerical simulation arises as a feasible alternative to overcome these limitations, which surely appear in the FSI wind analysis of lightweight structures.

1.4 Objectives

The main objective of this thesis is the development of a numerical tool, henceforth the Virtual Wind Tunnel (VWT), for the reliable resolution of FSI numerical experiments. It has to be emphasized that the VWT must be capable of dealing with any type of structure, namely volumetric and membrane-like bodies.

Considering that the target application is the resolution of extremely flexible structures under the action of severe wind loads, the VWT must be robust enough to cope with extremely large displacements and/or rotations. On top of this, there might be cases in which the analysed bodies can suffer from self-contact and wrinkling.

In this context, the traditional Arbitrary Lagrangian-Eulerian (ALE) approaches that have been traditionally used for similar purposes are prone to suffer from element deterioration. Hence, rather than focusing on the Computational Solid Mechanics (CSM) problem, this work focuses on the developing of embedded Finite Element (FE) methods for the robust resolution of CFD and FSI problems involving thin-walled structures. It is pointed out that efficiency is crucial for the effective technology transfer of the VWT to real engineering applications, in which the computational times must be in accordance to the project development ones. In consequence, all the techniques in this work take this as an essential requirement.

At the end of this thesis, the VWT tool will help to

• extend the use of CFD and FSI simulations to the civil engineering field as it already occurs in other engineering disciplines (e.g. aerospace engineering)
• study the aeroelasticity phenomena that occur in extremely flexible structures that cannot be tested in real wind tunnel facilities. The better understanding of these during the design phase would allow optimizing the geometry and materials of the structure, thus reducing the construction overall cost and material consumption.
• optimize and simplify the design process by reducing the amount of required wind tunnel experiments, something that should turn into a more efficient project development.

Last but not least, it is due mentioning that all the implementations of this work must remain as open source software. In order to fulfil this requirement the VWT will be developed within the Kratos Multiphysics open source framework [27,28].

1.5 Contents

First of all, it has to be clearly stated that this document is conceived as a compendium of publications. Therefore, three of the chapters are made as a brief introduction to the corresponding paper and the postprint submitted to the journal.

Chapter 2 presents the state of the art on fixed mesh CFD methods and FSI coupling. In Chapter 3 the quasi–incompressible Navier–Stokes stabilized formulation used throughout the thesis is derived. Chapter 4 is the first paper of the compendium and presents a discontinuous embedded formulation for the CFD analysis of slip thin-walled bodies. This formulation is enhanced in Chapter 5 to deal with both stick and slip thin-walled bodies. Chapter 6 is the second paper of the compendium and validates the formulation described in Chapter 5 with data from in vitro biomedical experiments. The extension of the embedded CFD solver to FSI problems is studied in the third paper of the compendium, which is collected in Chapter 7. In Chapter 8 embedded CFD methods are exploited to efficiently deal with extremely large wind engineering problems. Finally, the conclusions and further work lines are discussed in Chapter 9.

2 State of the art

This chapter reviews the state of the art on the use of fixed mesh techniques in CFD and FSI problems. The first section describes how this family of numerical methods can help to efficiently solve the VWT target applications. Secondly, the implicit representation of geometries is briefly discussed, particularly focusing on the discontinuous level set representation of membrane-like bodies. The third section reviews the characteristic features of each one of the literature fixed mesh approaches in order to choose the suitable one for the problems at hand. The forth section reviews the literature on the chosen alternative. The fifth section presents the common issues that appear when fixed mesh methods are applied to moving interface problems. Finally, the last section reviews the coupling techniques to approach the FSI problem, specially that one applied in this thesis.

2.1 Fixed mesh methods. Motivation

This thesis aims at creating a framework for the FSI analysis of any type of bodies, with a particular emphasis on the lightweight thin-walled structures. Taking into account that these are highly deformable structures, the traditional Arbitrary Lagrangian-Eulerian (ALE) methods [29] might lack robustness since the mesh motion problem tends to yield degenerated elements, namely highly distorted or even inverted, when large displacements and rotations of the boundaries appear. Moreover, problems involving thin-walled bodies might also fall into topology changes if self-contact or wrinkling appears.

Having said this, it can be easily guessed that the selection of the fluid-structure interface tracking technique becomes crucial for the robust and efficient resolution of the problems at hand. Taking this into account, mesh-based approaches can be classified according to how the analysed bodies are represented in body conforming (also known as body fitted) and non-conforming approaches.

Even though body conforming purely Lagrangian methods such as the Particle Finite Element Method (PFEM) [30,31] are a robust alternative for the interface tracking in moving boundary problems, they rapidly become prohibitive as remeshing is required at each time step. Alternatively, the CFD domain can be meshed regardless the analysed objects to decouple the computational mesh (henceforth denoted as background mesh) from the immersed geometries, which are no longer represented by the boundaries of the mesh. Instead, the analysed bodies are implicitly characterized on the background mesh from their skin representation by using techniques such as the level set method [32].

Although this introduces an extra complexity to the problem, the implicit representation of boundaries becomes advantageous for the interface tracking. Hence, all the potential mesh updating issues that may appear under the presence of large boundary displacements and rotations are bypassed since the ALE mesh motion problem is substituted by a simple update of the implicit representation of the bodies. On top of that, topology changes are naturally dealt with without requiring any operation in the background fixed mesh.

Furthermore, the use of an implicit representation also optimizes the simulation pipeline since the preprocess becomes much more efficient. On the one hand, it avoids all the manual repairing operations such as the removal of holes, overlapped surfaces or duplicated entities that are required to create a volume mesh from a “dirty” input geometry. These operations, which can amount up to the 50% of the total simulation time for real engineering cases [33], can be completely bypassed as the geometrical defects are filtered by the level set calculation algorithm, which is a more robust and effective operation [34,35]. On the other hand, it is also known that the creation of a volume mesh from a thin-walled geometry is a rather troublesome operation, specially when complex curved geometries are involved. This is because the original membrane geometry needs to be duplicated in order to generate the volume mesh. These two geometrical entities tend to intersect each other when working in the machine precision limit, thus leading to a failure of the meshing operation. Using an implicit representation of the volumeless immersed bodies automatically avoids such limitation.

2.2 Level set representation of immersed geometries

At this point, it is clear that using a level set representation is of great utility in some applications. This section briefly describes the basics of implicit representations and how they can be applied depending on the volumetric nature of the analysed geometries.

The main idea behind any level set method relies on the calculation of a signed distance function such that its zero isosurface matches, as closely as possible, the skin of the analysed bodies. Even though all the level set algorithms are based on the intersection between the background and overlapped skin meshes, they can be classified according to their treatment of the volumetric nature of the bodies. Hence, those bodies that feature a well defined internal volume are typically represented using a continuous level set algorithm (Fig. 2.1). In this case the distance is computed node-by-node and the inside/outside concept is treated using computer graphics techniques such as the raycasting [36]. These methods always return a smooth and continuous representation of the immersed bodies.

Figura 2.1: Continuous distance function. Body with a well-defined internal volume (left) and continuous distance representation (right). Dashed lines represent the continuous distance representation.

Unfortunately, these methods are not valid for membrane-like bodies as the inside/outside concept becomes meaningless in this case. In consequence, volumeless bodies cannot be represented in terms of a nodal-based distance function. Nevertheless, it is possible to describe these in terms of an element-based level set. As it can be observed in Fig. 2.2, such capability comes however at the price of having a potentially discontinuous distance function, meaning that the same node can have different distance sign and value depending on the element considered. This is exploited in the subsequent chapters, as well as in [37,38,39,40], to model a wide variety of thin-walled structures such as boat sails, lightweight membranes or thin biological tissues.

Figura 2.2: Discontinuous distance function. Body without internal volume (left) and discontinuous distance representation (right). Red and green portions of the cut elements indicate the positive and negative discontinuous distance regions. Light green denotes the non-intersected elements.

2.3 Level set methods for CFD and FSI problems

Prior to any further discussion, it is interesting to mention that the first FSI numerical approach was based on a non-conforming mesh method. It was presented in the early 70s by Prof. Peskin and relies on the addition of an artificial body force to consider the immersed structure motion into the fluid problem [41].

Introducing the use of a level set-based representation of the immersed geometries opens the possibility to use a wide variety of methods such as the eXtended Finite Element Method (xFEM) [42,43,44,45], the Immersed Boundary Method (IBM) [46,47,48,34,35,49], the Embedded Boundary Method (EBM) (also known as Cut-FEM) [50] or the Shifted Boundary Method (SBM) [51,52]. Regardless this thesis only considers traditional FE discretizations, the novel Immersogeometric approaches [16,53] that combine exact geometries and embedded formulations also deserve to be mentioned.

Table 2.1 compares the aforementioned methods according to different criteria. The first one is the computational effort in terms of locality and matrix graph preservation, according to which xFEM is the most expensive technique as it requires the introduction of extra Degrees Of Freedom (DOFs) associated with the FE enrichment unknowns. Aside of introducing the need for blending elements, which turns into a more expensive distributed memory assembly because of the communication overhead, this also requires the system matrix graph to be reconstructed each time the level set is updated.

The second criterion compares the BC imposition. In this regard, the IBM is possibly the most straightforward approach, as it allows direct imposition of the BC. However, this is not done over the interface but on its closest dry nodes, something that makes the method prone to suffer from mass conservation issues. The SBM method aims at overcoming this limitation by imposing a modified boundary value that takes into account the distance between the interface and the closest dry nodes where the BC is enforced. Contrariwise, the xFEM and the EBM rely on a weak BC imposition. This is normally done using variational techniques such as the penalty method [54,37], the Nitsche method [55,56,57,58,59,60] or the Lagrange multipliers method [61]. Although the weak imposition entails an extra complexity, that might turn into accuracy and stability issues, it allows to consistently enforce the BC over the level set intersections, which are the best non-conforming representation of the immersed geometries that can be implicitly achieved.

Last but not least, the table compares the capability of the methods to be applied to membrane-like bodies. As it can be observed, any method that requires an interior (i.e. dry node) for the BC imposition becomes no longer valid for the simulation of thin-walled geometries, meaning that only the xFEM and the EBM can be applied for such purpose.

Taking into account that industrial CFD and FSI problems likely involve a huge number of DOFs, the chosen alternative must be as computationally efficient as possible. Considering that only the xFEM and the EBM fulfill the thin-walled bodies requirement, the more suitable method for the target applications is the EBM owing to the absence of blending elements and the ability to keep the initial system matrix graph during the entire simulation.

2.4 The Embedded Boundary Method. An overview

After having selected the EBM as reference fixed mesh approach, this section briefly reviews its basis as well as the different alternatives that can be employed for the weak imposition of the level set BC. To ease the discussion, from now on the fluid computational domain is denoted as ${\textstyle \Omega _{f}}$ while the immersed objects skin representation is denoted as ${\textstyle \Gamma _{f}}$.

The first distinguishable feature of the EBM is that the computational domain is not conformed by the element faces but by the intersections of the level set function with the background mesh elements. Aside of being a more accurate representation of the analysed bodies, this enables the application of the method to thin-walled bodies. This is schematically depicted in Fig. 2.3, which recovers the previously presented volumetric (Fig. 2.1) and volumeless (Fig. 2.2) toy examples to show their corresponding embedded computational domain.

(a) Continuous level set case.	(b) Discontinuous level set case.
Figura 2.3: EBM sketch for the examples in Figs. 2.1 and 2.2. Light orange denotes the active computational domain (${\displaystyle \Omega _{f}}$). Red lines denote the level set intersections (${\displaystyle \Gamma _{f}}$) that represent the embedded skin.

From the previous images it can be deduced that an auxiliary splitting operation is required in each intersected element to perform the integration properly. As it is depicted in Fig. 2.4, the division of the parent geometry in a set of smaller ones allows relocating the integration points such that only the “wet” part of the cut elements is integrated. The same procedure is equivalently applied to thin-walled bodies (Fig. 2.4b) to consistently integrate both (“wet”) sides of the membrane bodies.

Once the volume integration has been performed, the fluid BC needs to be imposed over the either continuous or discontinuous level set intersections that conform ${\textstyle \Gamma _{f}}$ (highlighted in red in Figs. 2.3 and 2.4). Although it is possible to create new nodes from the splitting pattern to strongly impose the BC [34], this operation is commonly done in a weak (i.e. variational) sense. There is a wide variety of methods that can be used for such purpose, leading to different flavors of the EBM. In the following a brief review of all the approaches that have been successfully applied to CFD problems is presented.

(a) Continuous level set case.	(b) Discontinuous level set case.
Figura 2.4: EBM single element sketch. Light orange denotes the active part of the element belonging to ${\displaystyle \Omega _{f}}$. Red lines denote the level set intersections ${\displaystyle \Gamma _{f}}$ that represent the embedded skin. Black dashed lines represent the auxiliary divisions for the integration. The black cross markers represent the outside nodes. Light green cross markers represent the fluid nodes. Light green round markers represent the integration points.

The most extended approach to weakly enforce the BCs in the context of EBM methods is possibly the Nitsche method [55,62], which has been applied for the imposition of the stick (no-slip) condition in both Stokes [58] and Navier-Stokes [16,59] problems. The modified Nitsche method for no-slip conditions in [56] should be also acknowledged because of its simplicity and remarkable good performance. Regarding the slip wall behavior, the Nitsche method is used in [57] to enforce a pure slip condition. In [60] this is extended to the Navier-slip condition.

The Lagrange multipliers method can be alternatively used to impose Dirichlet BCs [63,64]. In [57] a Lagrange multiplier approach for the imposition of the pure slip condition in Stokes flows is presented. Furthermore, the Lagrange multipliers technique in [61], which is based on a static condensation, deserves to be highlighted after its ability to keep constant the number of DOFs.

Unfortunately, all the aforementioned EBM approaches require a well-defined internal volume. In consequence, none of the methods above can be directly applied for the analysis of thin-walled bodies, thus precluding the simulation of membrane and shell structures. This could be achieved in an IBM fashion by adding an artificial volume force body [41,65,66,67]. However this kind of methods might lack accuracy and suffer from mass conservation issues in some applications. In [9] an EMB ghost cell method able to work with any type of bodies (volumetric and volumeless) is presented. Nonetheless, such capability comes at the price of having a no longer purely local method as the BC imposition requires the element neighbors. Alternatively, the approach in [54] addresses the problem in a EBM-like manner by introducing a modified FE space in combination with a penalty imposition of the no penetration BC. An important advantage of this method is that the graph of the discrete stiffness matrix is preserved, something that becomes crucial when the level set function needs to constantly updated (i.e. FSI problems).

After having reviewed the literature on Cut-FEM approaches it can be noted that the volumetric problem is pretty well resolved. There is a wide variety of approaches that can efficiently and robustly deal with the problem at hand. Among these, the approaches in [56] and [60] are used in this work for the volumetric no-slip and slip (in a general Navier-slip sense) wall behavior modeling.

However, there is still room for the research on the application of the EBM to the resolution of CFD problems involving thin-walled bodies. More specifically, there is no purely local approach in the literature capable of representing the discontinuities arising from the immersion of membrane-like structures. The chapters 4 and 5 of this thesis aim at accomplishing such objective.

2.5 Fixed mesh methods and moving bodies

The enhanced robustness that fixed mesh approach have to deal with arbitrary large displacements and rotations comes at the price of some particularities that require to be taken care of. These are the small cut instability and the historical data initialization.

2.5.1 Small cut instability

As it can be easily guessed from its name, this problem is an instability that pops up when there is an ill-conditioned intersection between the level set and the background mesh. Even though the concept of ill-conditioned intersection varies depending on the fixed mesh approach used, it can be generalized to any intersection of the embedded skin that occurs extremely close to one node of the background mesh, which is equivalent to say that close to zero values of the level set are potentially ill-conditioned. It is important to note that no distinction is done yet between continuous and discontinuous level set functions, meaning that zero values of the distance function are always ill-defined. For those cases involving volumetric bodies represented by continuous distance functions, the previous definition can be restricted to those level set values close to zero that yield an almost empty element.

The small cut instability commonly appears when the interface moves across the background mesh. Nonetheless, it may also happen in problems involving a steady interface, specially in those ones involving poorly defined input geometries. In consequence, an auxiliary strategy to prevent this undesired situation is always demanded to robustly apply any fixed mesh method.

There are several approaches to prevent the convergence of the problem to be ruined by the ill-conditioned cuts. Some of these are the ghost penalty method [50] or the ghost cell method [68]. Although the stability and robustness of these methods is proved, they require the neighbors of the ill-conditioned elements to be applied. This adds a search operation that turns into an extra computational overhead, which might be particularly relevant when distributed memory platforms are used.

Alternatively, the ill-conditioned cuts can be avoided rather than controlled by using a level set quality check and correction algorithm. This approach starts from the idea that the level set is indeed a representation of the real geometry, whose quality is directly related to the background mesh element size. Taking this into consideration, the distance can be slightly modified by a small quantity depending on the element size without incurring in a noticeable error. For the continuous level set case, this means to avoid the “almost empty” intersections (highlighted in light orange in Fig. 2.5) by slightly moving the interface towards the positive distance side. Likewise, for the discontinuous level set case ill-conditioned cuts can be avoided by slightly modifying the close to zero elemental level set values (Fig. 2.6). This approach arises as a more computationally efficient (almost) purely local alternative that only requires one synchronization operation. The capability of this apporach to robustly deal with the problems at hand is evinced in Chapters 6, 7 and 8, where it is applied to complex real-life geometries.

Figura 2.5: Continuous distance check and correction. Grey regions represent the computational domain ${\displaystyle \Omega _{f}}$. Red dashed lines represent the level set intersections ${\displaystyle \Gamma _{f}}$. Before (left) and after (right).

Figura 2.6: Discontinuous distance check and correction. Grey regions represent the computational domain ${\displaystyle \Omega _{f}}$. Red dashed lines represent the level set intersections ${\displaystyle \Gamma _{f}}$. Before (left) and after (right).

2.5.2 Consistent initialization of values. The FM-ALE method

This issue appears when the level set function, regardless its continuous/discontinuous nature, evolves in time. As it is graphically described in Fig. 2.7, when the distance function modifies its spatial configuration the background mesh nodes that used to lie in one side of the level set might end up in the other one. In consequence, the historical data of these nodes is no longer valid, thus leading to a wrong approximation of all the inertial terms involved in the problems.

It is due mentioning that no distinction is done between continuous and discontinuous distance functions, as the historical values need to be consistently updated in both cases. While in the continuous level set case (volumetric bodies) it is needed to consistently initialize the historical values in those nodes that change their status from “dry” to “wet”, in the discontinuous level set case (volumeless bodies) it is required to do so as the historical data becomes meaningless as soon as the nodes move from one side to the other of the immersed thin-walled body. Having said this, it can be easily guessed that the proper treatment of this issue is crucial for the application of fixed mesh methods to any moving boundaries problem (i.e. FSI).

Figura 2.7: Embedded nodal initialization example. The node highlighted in green changes its position from one side to the other of the level set (red dashed line), which moves as indicated by the black arrow.

The consistent historical data initialization can be achieved in an ALE-fashion by using the Fixed Mesh–Arbitrary Lagrangian Eulerian (FM–ALE) method. The FM–ALE was firstly described in [69] to address this issue in moving boundaries embedded CFD problems. In [70] the authors extend the method to solid mechanics and FSI problems. Similarly, they do so for free-surface problems involving floating objects in [71,72].

The main idea behind the FM-ALE method is to retrieve the consistent historical values from the so called virtual mesh, which is created as a copy of the background mesh used in the resolution of the problem. Each time the level set changes its spatial configuration, the virtual mesh is accordingly deformed by solving a mesh motion problem. Then all the historical values, as well as the mesh velocity, can be projected from the virtual mesh to the background one. Once the projection is completed, the deformation of the virtual mesh is reverted back to the initial position.

It is due mentioning that this last step is the key of success of the method when it is used in presence of arbitrary large movements of the embedded interfaces. Hence, bringing back the virtual mesh to its original position each time the FM–ALE algorithm is applied implies that the mesh motion problem most likely involves small displacements and rotations, thus avoiding the common element distortion issues associated to body conforming ALE approaches. Besides that, reverting the mesh deformation also makes the problem presumably cheaper in terms of the computational effort.

2.6 Fluid-Structure Interaction

This section reviews the ingredients that are required to extend the embedded CFD solver to FSI problems. Hence, the different alternatives that can be selected to achieve the coupling are reviewed. The last part of the section focuses in detail on the particular choice of this work. To ease the discussion, the fluid and structure domains are denoted as ${\textstyle \Omega _{f}}$ and ${\textstyle \Omega _{s}}$ from now on. Similarly, their corresponding FSI interfaces are denoted as ${\textstyle \Gamma _{f}}$ and ${\textstyle \Gamma _{s}}$.

2.6.1 Coupling schemes

For a coupled problem of any type it is required to satisfy the interface transmission conditions to guarantee that the coupling is effectively achieved. Such transmission conditions are of course problem dependent and thus need to be particularized depending on the physics to be coupled. For the FSI specific case these are the interface force equilibrium as well as the mass continuity, which demands that ${\textstyle \Gamma _{f}}$ and ${\textstyle \Gamma _{s}}$ spatial configurations match.

To ensure that the transmission conditions of the problem are accomplished, a coupling scheme is required. These can be classified according to how the transmission conditions are imposed in Dirichlet-Neumann (DN), Neumann-Neumann (NN) and Robin-type schemes (i.e. Robin-Robin (RR) or Neumann-Robin (NR)).

Owing to its proved well performance in the FSI field [73,74,75,54], the DN coupling scheme is selected as reference coupling scheme in this thesis. In short, the DN scheme consists in enforcing a Dirichlet BC in one subdomain and a Neumann BC in the other one. The Dirichlet condition is commonly applied in the lower density subdomain, which in the FSI case is likely ${\textstyle \Omega _{f}}$. In consequence, the Neumann BC is applied over ${\textstyle \Gamma _{s}}$ as a surface load whose value comes from the stress in the fluid boundary. For a thorough review on the stability and performance of the DN scheme, as well as on its interaction with several time discretization techniques the reader is referred to [76].

Furthermore, the implementation of the DN scheme is reasonably simple and does not require to modify the subdomain formulations. This last feature becomes crucial for the efficient application of non-intrusive coupling algorithms like the ones that are to be exploited in chapter 7 of this thesis.

Nevertheless, it has to be mentioned that there exist cases in which the DN may lack convergence or suffer from stability issues [77]. In this context, NN and Robin-type schemes arise as a feasible alternative to overcome these limitations. Although these applications are out of the scope of this thesis, the reader is referred to [78,79] and [80,81] for more information regarding RR and NR schemes.

2.6.2 Implementation of coupling schemes

The previously defined transmission conditions require a coupling strategy to ensure that they are satisfied in all the coupling interfaces. Such strategies can be classified according to how the coupled problem is solved in monolithic and partitioned (also known as staggered) strategies.

The idea behind any monolithic strategy is to put together all the subdomain contributions, together with the corresponding transmission conditions, in a unique large system of equations. This approach has been successfully applied in the resolution of body conforming FSI problems in [82,30]. More related to this thesis is the method in [70], where the authors present a monolithic FSI technique based on a non-conforming discretization in both the fluid and the structure subdomains.

Although monolithic strategies are possibly the most robust approach for the resolution of strongly coupled FSI problems, they tend to yield poorly conditioned systems of equations, which might eventually impede the use of fast iterative solvers. Besides this, the implementation of a monolithic strategy is rather intrusive and likely requires the development of a new solver for the specific purpose of the coupling.

Partitioned strategies avoid these disadvantages by keeping a separated solver for each one of the subdomain problems, meaning that the fulfillment of the transmission conditions is achieved by exchanging information between the coupled subdomains. Although the convergence rate of partitioned strategies might be suboptimal, they allow reusing already existing, possibly robust and thoroughly validated, specific solvers. This idea brings up the concept of black-box coupling [83,84,85], in which the coupling strategy is also requested to not interfere with the subdomain solvers more than to get and set information. Hence, the subdomain solvers are considered as a sort of black boxes that return a solution for the input data, that must be provided according to their Application Programming Interface (API).

In the framework of this thesis, this makes possible to directly use the already existent structural mechanics module of Kratos Multiphysics in order to focus the implementation effort on the development of the embedded CFD technologies and its extension to FSI problems.

Partitioned strategies can be further classified according to how many times the information is exchanged per time step. On the one hand, there are the loosely-coupled (or explicit) strategies in which the information is exchanged once per time step [86]. This is obviously the cheapest approach in terms of computational and implementation effort. However, they might lack convergence when there is strong interaction between subdomains or when the added mass effects are important [87,88,89].

On the other hand, there are the strongly-coupled (or implicit) strategies in which the information exchange occurs in an iterative manner until a certain convergence criterion is reached. Depending on how the information is used within each iteration, it is possible to distinguish between Jacobi and Gauss–Seidel strongly-coupled strategies. While in Jacobi-type iterations the information exchange occurs at the same time, allowing the parallel resolution of the subdomain problems, in Gauss-Seidel iterations this is done in a sequential manner, enabling the use of the latest information obtained in one subdomain to solve the other one within the same iteration.

The iterative nature of Gauss-Seidel approaches makes them more suitable for the target applications of this thesis. Although their convergence rate might be not as good as that of monolithic approaches, it can be astonishingly improved, specially when facing strongly coupled problems, by using them in combination with convergence acceleration algorithms.

Finally, it is interesting to point out that there exist some semi-implicit strategies that are somehow between the explicit and implicit ones. A very detailed review about several semi-implicit approaches can be found in [81].

2.6.3 Black-box interface residual minimization techniques

The choice of a Gauss–Seidel iterative strategy requires the definition of an interface residual function such that its minimization ensures the fulfilment of the problem transmission conditions. As the interface residual is completely problem dependent, it is assumed to be already defined for the discussion at hand. Hence, the literature review in this subsection focuses on the residual minimization rather than on its definition, which will be addressed in chapter 7.

As it is commented above, the use of a black-box coupling strategy implies restricted access to the subdomain solvers. The most straightforward approach that accomplishes with such requirement is a fixed-point iteration with relaxation. Although convergence can be achieved using a constant relaxation parameter, these methods are commonly used in combination with dynamic relaxation formulas. Among all of them, the second order Aitken scheme is possibly the most extended technique in the FSI community [73,75] owing to its reasonably good performance and relatively easy implementation.

Alternatively, the interface residual problem can be considered as a vector of unknowns to be minimized. Hence, the problem can be posed as a Newton–Raphson (N–R) iterative procedure, leading to the so called Jacobian-based residual minimization techniques. Ideally, the exact Jacobian that minimizes the vector interface residual could be computed by accessing the subdomain solvers. However, this intentionally precluded by the selection of a black-box coupling strategy. Nevertheless, the exact Jacobian calculation can be avoided by using Jacobian-free Newton-Krylov (JFNK) [90,14] or Quasi-Newton (QN) [84] methods.

Rather than approximating the complete interface Jacobian, JFNK methods are based on approximating its projection onto the iteration update vector by using a finite differences formula. Owing to the inherent non-linearity of the FSI problem, such projection needs to be linearized by introducing a small perturbation, which needs to be selected by the user. Although this issue can be addressed by using a formula to automatically do that [91], the optimal value for the small perturbation is still completely problem dependent.

In this context, QN methods arise as an alternative to avoid the selection of the perturbation constant. Unlike JFNK methods, QN approaches approximate the complete interface Jacobian by means of a linearized formula that takes the information from the previous FSI iterations. The different QN methods can be characterized according to how such Jacobian approximation formula is conceived [92].

Among all of them, the Broyden's iteration scheme should be acknowledged as the first QN method applied to FSI problems [91]. A more popular approach is the Interface Quasi-Newton with Inverse Jacobian from Least Squares model (IQN-ILS) [93,92], as well as its block equations version (IBQN-ILS) [92]. In [92] the authors prove that the performance of the IQN-ILS overcomes that of Aitken and of Interface-GMRES schemes for all the reported examples. A very similar approach to the IQN-ILS is the MultiVector Quasi-Newton (MVQN) algorithm presented in [84]. Although the authors report that it converges slightly faster than the IQN-ILS, this comes at the price of requiring a matrix inversion, whose size equals the interface residual one, each time the Jacobian is approximated. This disadvantage has been recently overcome in [94].

Considering that QN approaches involve a rather large number of matrix operations, each QN iteration is possibly more computationally expensive than that one of a JFNK technique. Notwithstanding this, they are reported to be the most efficient Jacobian-free approach in terms of total number of iterations. Taking into account that the computational bottleneck in partitioned FSI problems is likely the resolution of the CFD problem, the reduction of the total coupling iterations is what effectively turns into a performance improvement. Taking this into account, the embedded FSI coupling Gauss-Seidel technique presented in chapter 7 will be based on the use of QN algorithms. More specifically, the MVQN is chosen as reference black-box convergence accelerator owing to its reported good performance.

input-ch2-StateArt input-ch3-NavierStokes input-ch4-CMAME_emb_1 input-ch5-AusasNavierSlip input-ch6-BIO input-ch7-CMAME_fsi

3 Quasi-incompressible Navier–Stokes stabilized formulation

This section presents the quasi–incompressible Navier-Stokes formulation that serves as basis for all the embedded CFD developments of this thesis. Firstly, the viscous incompressible Navier-Stokes governing equations as well as the inclusion an extra pseudo–compressibility term are discussed. Secondly, the FE discretization and the subscales stabilization technique are introduced. In the third section the variational form of the problem is derived together with some notes on its symbolic implementation. Finally, the correctness of the implementation is proven by the method of manufactured solutions.

3.1 The Quasi-incompressible Navier–Stokes equations

This thesis focuses on the resolution of CFD and FSI problems involving viscous fluids. Considering that the Mach number (Ma) of the target applications is below 0.3, the flow can be considered as incompressible and thus modeled with the well-known viscous incompressible Navier–Stokes (N–S) equations [29,95].

In this context, the Cauchy stress tensor ${\textstyle {\boldsymbol {\sigma }}}$ can be defined as

where ${\textstyle \mathbf {u} }$ denotes the velocity, ${\textstyle p}$ the pressure, ${\textstyle {\boldsymbol {\nabla }}^{s}}$ the symmetric gradient operator, ${\textstyle \mathbf {I} }$ the 2nd order identity tensor and ${\textstyle \mathbb {C} }$ the 4th order viscous constitutive tensor.

Although any viscous constitutive relation could be used, only isotropic Newtonian fluids are considered in this thesis. Hence, the viscous constitutive tensor reads

where ${\textstyle \otimes }$ denotes the tensor product and ${\textstyle \mathbb {I} }$ the 4th order symmetric identity tensor. ${\textstyle \mu }$ and ${\textstyle \lambda }$ are the viscosity constants. The first one represents the dynamic viscosity while the second one can be computed as ${\textstyle \lambda =-2\mu /3}$ after taking the Stokes' assumption, which is valid in those cases in which the compressibility of the medium is negligible.

Inserting the previous definition of ${\textstyle {\boldsymbol {\sigma }}}$ into the balance of linear momentum and mass conservation equations yields the viscous incompressible N–S equations

being ${\textstyle \rho }$ the density and ${\textstyle \mathbf {b} }$ the vector of volumetric (body) forces. ${\textstyle \partial (\bullet )/\partial t}$ and ${\textstyle {\boldsymbol {\nabla }}}$ denote the partial time derivative and the gradient operator. ${\textstyle \mathbf {a} }$ is the convective velocity, which is defined from the mesh velocity ${\textstyle \mathbf {u} _{m}}$ as ${\textstyle \mathbf {a} =\mathbf {u} -\mathbf {u} _{m}}$.

At this point, it is important to provide some remarks about the convective term in Eq. 3.3.a since it strongly affects the applicability and convergence of the method. The first one concerns the inclusion of ${\textstyle \mathbf {u} _{m}}$, which equals ${\textstyle \mathbf {0} }$ in a fixed mesh Eulerian framework. Although this is a valid simplification, it would impede the application of the FM–ALE algorithm in chapter 7. In consequence, it is convenient to keep such ALE definition of the convective velocity ${\textstyle \mathbf {a} }$.

The second one concerns the non–linearity. As it is evinced in [96], it is convenient to do a Picard linearization of the convective term for the sake of the non-linear iteration robustness. Hence, the convective velocity is linearized as ${\textstyle \mathbf {a} =\mathbf {u} ^{k-1}-\mathbf {u} _{m}}$, where ${\textstyle k}$ superscript denote the current iteration. In consequence, the convective term turns into ${\textstyle \left(\mathbf {u} ^{k-1}-\mathbf {u} _{m}\right)\cdot {\boldsymbol {\nabla }}\mathbf {u} ^{k}}$, which is the correct Picard linearization that ensures stability and proper convection of information [96].

With regard to the mass conservation equation, when dealing with purely incompressible fluids it is customary assumed that ${\textstyle {\frac {D\rho }{Dt}}=0}$, meaning that that the pressure field is defined up to a constant. In the standard body fitted case, this issue is immediately fixed by imposing a Neumann BC in any part of the domain. Unfortunately, this solution might be no longer valid when working with a level set-based fixed mesh approach. This is due to the fact that, despite the level set is a robust and always defined operation, there are situations in which isolated closed domains of fluid (i.e. “bubbles”) can appear. Taking into account that such undesired fluid cavities are likely related to poorly defined input files (e.g. stl) or to moving boundary problems involving complex geometries, it is impossible to a priori known their location. Trying to locate them in order to apply a Neumann BC is not a feasible solution since this is an extremely expensive, and possibly non-deterministic operation.

Alternatively, the unbounded pressure issue can be readily fixed if a slight compressibility is included in the formulation. This can be achieved by rewriting the density time derivative in Eq. 3.3.b such that a state equation linking the density to the pressure field can be introduced. For almost incompressible fluids, it is possible to take the simplified speed of sound equation of state

Hence, developing the material time derivative of the density in Eq. 3.3.b as

which is equivalent to

allows to get rid of the density partial time derivative by inserting Eq. 3.4 as

By further assuming that ${\textstyle {\boldsymbol {\nabla }}\rho \approx \mathbf {0} }$, which physically means that the density space fluctuations are negligible, the term ${\textstyle \mathbf {a} \cdot {\boldsymbol {\nabla }}\rho }$ can be neglected to give

Combining this last equation with the momentum conservation law in Eq. 3.3.a yields the final form of the quasi-incompressible N–S equations

It is important to mention that the fully-incompressible form ${\textstyle \nabla \cdot \mathbf {u} =0}$ is recovered as the sound speed ${\textstyle c\rightarrow \infty }$. This makes possible to effectively deactivate the pseudo-compressible term in those cases in which isolated fluid cavities are not expected to appear. In chapter 5 a test example that proves the capability of this approach to deal with such undesired scenarios is presented.

3.2 Discrete form and stabilization

The FE discretization of mixed formulations, such as the Stokes (or N–S) equations, is commonly associated to saddle point problems [97]. These give rise to numerical artifacts and spurious oscillations, which most probably ruin the solution, when the discretizations of the mixed variables are not compatible between them. In this context, the Ladyzhenskaya–Babuska–Brezzi (LBB) condition, also known as inf–sup condition, provides a criterion to distinguish between stable and unstable discretizations. For more details on inf–sup stable and unstable FE discretizations, the reader is referred to [98,95].

Aiming to keep the embedded elemental splitting as simple as possible, this thesis only considers simplex elements with same velocity and pressure interpolation order (i.e. linear triangles and tetrahedra [99]). Furthermore, simplex elements are the most straightforward approach to generate unstructured grids around complex geometries. However, they are inf–sup unstable and therefore require the use of a stabilization technique to avoid spurious oscillations. Note that using an inf–sup stable discretization would require the use of higher order approximations, which in the context of EBM would turn into an extremely cumbersome implementation of the splitting algorithms.

There is a wide variety of stabilization techniques that have been proved to be effective for similar N–S problems. Some of these are the Streamline Upwind/Petrov–Galerkin (SUPG) [100,101], the Galerkin Least-Squares (GLS) [102], the Finite Increment Calculus (FIC) [103] or the Variational MultiScales (VMS) [104,105], which is the adopted one.

The VMS method relies on the separation of the solution fields ${\textstyle \mathbf {u} }$ and ${\textstyle p}$ in two scales as

${\textstyle \mathbf {u} _{h}}$ and ${\textstyle p_{h}}$ being the velocity and pressure FE resolvable scales. Contrariwise, ${\textstyle \mathbf {u} _{s}}$ and ${\textstyle p_{s}}$ are the subscales representing the velocity and pressure fluctuations that cannot be captured by the FE solution and thus need to be modeled.

The previous separation of scales can be inserted into Eq. 3.6 to yield

At this point it is required to define a model for the subscales ${\textstyle \mathbf {u} _{s}}$ and ${\textstyle p_{s}}$, which are assumed to have a bubble shape that implies zero value of all their boundary integrals [104]. In general terms, it can be said that the subscales are built as a projection of the FE residuals onto the space of the small scales. In consequence, the nature of such projection defines the subscale model to end up with. Hence, assuming an identity projection operator leads to the Algebraic Sub-Grid Scales (ASGS) [106] while considering an orthogonal one, which implies the assumption that the subscales space is orthogonal to the FE one, leads to the Orthogonal SubScales(OSS) [107]. Owing to its remarkably good performance and simplicity, which turns into a more efficient implementation, the ASGS is the option selected in this work.

Hence, the ASGS velocity and pressure subscales can be calculated from the momentum and mass conservation discrete residuals ${\textstyle \mathbf {R^{M}} }$ and ${\textstyle R^{C}}$ as

being ${\textstyle \tau _{1}}$ and ${\textstyle \tau _{2}}$ the stabilization constants, which are taken from [106] and computed as

where ${\textstyle c_{1}=4.0}$ and ${\textstyle c_{2}=2.0}$ , ${\textstyle \lVert \mathbf {a} \rVert }$ is the convective velocity norm and ${\textstyle h}$ is the element size, which is computed as the average of the heights associated to each node of the element. ${\textstyle \tau _{dyn}}$ is a parameter bounded between 0 and 1 that allows considering the time step value in the calculation of ${\textstyle \tau _{1}}$.

Note that the previous subscales definition (Eq. 3.9) introduces a dependency on the FE solution through the algebraic momentum and mass conservation residuals, which are defined as

Besides this, Eq. 3.8 also shows the time dependent nature of the subscales. Considering such time dependency leads to a so called dynamic subscales formulation. Conversely, neglecting it by assuming that ${\textstyle \partial \mathbf {u} _{s}/\partial t\approx \mathbf {0} }$ as well as ${\textstyle \partial p_{s}/\partial t\approx 0}$ leads to a quasi-static subscales formulation. Even though the dynamic approach has superior characteristics [108], these come at the cost of increasing the complexity of the formulation. Moreover, they also involve a huge increase in the memory consumption as they introduce the need of storing the subscales historical data for each integration point. Taking these considerations into account, the quasi-static approach is the preferred option in this thesis.

3.3 Variational form and automatic differentiation

After having defined the stabilized governing equations, the variational (weak) form of the problem needs to be derived. To do that, the velocity and pressure test functions ${\textstyle \mathbf {w} }$ and ${\textstyle q}$ are defined. For the sake of simplicity, it is also convenient to use the standard notation ${\textstyle (\bullet ,\bullet )_{\Omega }}$ for the inner product volume integral. Likewise, the inner product boundary integrals can be denoted as ${\textstyle \left\langle \bullet ,\bullet \right\rangle _{\Gamma }}$.

Hence, the Galerkin FE functional to be minimized is built as the inner product volume integral of the test functions times the FE residuals as

Note that Eq. 3.12 involves the complete FE residuals, which include not only the FE solution but also the subscales. After introducing the scales separation and the quasi-static subscales assumption, the complete FE residuals ${\textstyle \mathbf {R^{M}(\mathbf {u} ,p)} }$ and ${\textstyle R^{C}(\mathbf {u} ,p)}$ read as

Inserting these into the Galerkin functional above yields the stabilized quasi–incompressible Navier–Stokes Galerkin functional

Taking into account that this formulation is to be discretized using simplicial elements, integration by parts is demanded for getting rid of as higher order derivatives as possible. After doing so and dropping all the remaining higher order terms, the final Galerkin functional to be implemented reads as

being ${\textstyle \mathbf {t} }$ the Cauchy traction vector, computed as ${\textstyle \mathbf {t} =\left(\mathbb {C} :{\boldsymbol {\nabla }}^{s}\mathbf {u} _{h}-p_{h}\mathbf {I} \right)\cdot \mathbf {n} }$, and ${\textstyle \mathbf {n} }$ the outwards unit normal vector. Note that this equation already takes into account the subscales null boundary value assumption.

The algebraic form of the discrete functional in Eq. 3.15 can be automatically obtained by using a Computer Algebra System (CAS). Among all the available alternatives, in this work the Python library Sympy [109] is used. Aside of being a much more efficient alternative to the traditional implementation in terms of the human effort, following a symbolic implementation ensures that the differentiation of the functional is always correct since this is automatically performed by the CAS library. This limits the possible error sources to the definition of the Galerkin functional, and thus minimizes the potential implementation mistakes.

Therefore, the elemental Left Hand Side matrix (${\textstyle \mathbf {LHS} }$) and Right Hand Side vector (${\textstyle \mathbf {RHS} }$) can be automatically obtained by expressing the symbolic functional in terms of the nodal test functions (${\textstyle \mathbf {w} _{I}}$ and ${\textstyle q_{I}}$) and the nodal unknowns (${\textstyle \mathbf {u} _{I}}$ and ${\textstyle p_{I}}$). Then, assuming a symbolic description of the shape functions and of their derivatives, allows obtaining the elemental ${\textstyle \mathbf {RHS} }$ by automatic differentiation as

while the elemental ${\textstyle \mathbf {LHS} }$ can be similarly derived as

The symbolic implementation of the quasi–incompressible Navier–Stokes stabilized formulation presented in this chapter is available in the Kratos Multiphysics repository.

Last but not least, it is due mentioning that this section assumes previous knowledge on the FE method. In case more details about the basis of the method are required, the reader is referred to [110,99]. Concerning the specific application of the FE method to CFD problems, the reader is referred to [111,29,95].

3.4 Validation

This section validates the FE implementation of the quasi–incompressible Navier–Stokes formulation previously described. To that purpose the Method of Manufactured Solutions (MMS) is used [112,113]. The problem geometry consists in a 1x1m square while the material properties ${\textstyle \rho }$ and ${\textstyle \mu }$ are equal to 10³ and 10^-4. The speed of sound velocity ${\textstyle c}$ is set to 10⁴ in order to introduce the effect of the pseudo–compressible term.

The problem is solved for two different source terms that correspond to two different transient solutions. First the source term

that corresponds to the sinusoidal transient field

is set.

Secondly, the experiment is repeated for the source term

which corresponds to the non-linear transient field

Fig. 3.1 collects the weighted average error for both the velocity and pressure fields after 1s of simulation. In both cases it can be observed that the velocity and pressure convergence rates match the expected ${\textstyle h^{2}}$ and ${\textstyle h}$ ones.

(a) Sinusoidal transient field.	(b) Non-linear transient field.
Figura 3.1: Manufactured solution experiments results.

4 A modified Finite Element formulation for the imposition of the slip boundary condition over embedded volumeless geometries

4.1 Article data

Title: A modified Finite Element formulation for the imposition of the slip boundary condition over embedded volumeless geometries

Authors: R. Zorrilla, A. Larese and R. Rossi

Journal: Computer Methods in Applied Mechanics and Engineering 353 (2019) 123–157

Received: 2 May 2018 / Accepted: 10 May 2019 / Available online: 15 May 2019

DOI: 10.1016/j.cma.2019.05.007

4.2 Scientific contribution

This chapter collects the first paper of the compendium, which presents a novel embedded mesh method for the resolution of Navier–Stokes problems. Even though the proposed formulation is capable of dealing with any type of body, it is intended to be used in combination with discontinuous level set functions to allow the resolution of CFD problems involving embedded thin-walled bodies.

The proposal is based on the substitution of the standard FE space by an alternative one in those elements that are intersected by the level set. Such alternative space, which was firstly proposed in [114], makes possible the representation of both the velocity and pressure discontinuities arising from the immersion of any body. In addition to this, a no penetration condition is weakly imposed by doing a convenient integration by parts of the mass conservation equation and the addition of an extra term that penalizes the normal projection of the embedded wall velocity. Taking into account that no constraint is enforced in the tangential direction, the imposition is equivalent to a slip BC.

The paper also discusses alternative approaches to impose the slip condition, namely a MultiFreedom Constraint (MFC) approach for body fitted discretizations and a Nitsche-based imposition of the Navier-slip BC for volumetric embedded bodies. These alternative approaches are also applied to validate the discontinuous embedded proposal. The accuracy of the method is proved as well as its ability to keep the ${\textstyle h^{3/2}}$ convergence order of the modified FE space. Furthermore, the chapter also showcases a potential real application that the proposed technique could have in the context of membrane structures CFD by solving the flow around two boat sails.

The results obtained in this chapter are a very first step towards the final VWT objective. More specifically, this chapter settles the basis of the discontinuous embedded Navier–Stokes formulations that are to be exploited in the subsequent chapters for the resolution of more complex CFD and FSI problems.

5 Discontinuous Nitsche-based formulation for the imposition of the Navier-slip condition

5.1 Introduction

This chapter advances further in the investigation of embedded CFD methods for thin-walled bodies. Hence, taking as starting point the purely-slip formulation presented in Chapter 4, this chapter introduces some enhancements in order to extend the applicability range of the method to any wall behavior from the slip to the stick (no-slip) limits. This is achieved by substituting the slip BC by a general Navier-slip BC, which can be somehow understood as a linear wall law. In addition, the BC enforcement is also improved by substituting the penalty-based method by a Nitsche imposition.

The chapter is organized as follows. First of all the discontinuous Nitsche-based method for the Navier-slip BC imposition is presented. Secondly, the new proposal is validated by solving several examples, which involve different wall behaviors ranging from the slip to the no-slip limits. Moreover, a potential industrial application is also presented. The last section collects the conclusions and further enhancements.

The formulation in this chapter has been submitted for its publication in Computer Methods in Applied Mechanics and Engineering (CMAME). The current status of the paper is under review.

5.2 Embedded Navier-slip boundary condition

This section presents a novel formulation for the imposition of the Navier-slip BC over embedded thin-walled bodies. Hence, the substitution of the standard FE space by the Ausas discontinuous one is also introduced in order to represent the discontinuities arising from the immersion of a (possibly thin-walled) body. However, the space substitution is combined in this case with the Navier-slip Nitsche method presented in [60]. Although this might seem a simple mixture of the ideas in Chapter 4 with the ones in [60], it is not obvious that the combination of the Ausas FE space with a stick wall behavior works out of the box owing to the worse interpolation properties featured by the Ausas FE space. In the following, the properties of the Ausas FE space are briefly recalled to ease the comprehension of the rest of the chapter. The next subsections present the Nitsche method for the imposition of the Navier-slip BC and a more accurate alternative to calculate the drag force when the Ausas FE space is involved.

5.2.1 The Ausas Finite Element space

The Ausas FE space was initially conceived to represent the discontinuities arising from the resolution of two-phase flow problems [114,115]. In this thesis, the ability of the Ausas FE space to disconnect the solution fields in the two sides of the level set intersection is exploited to directly embed thin-walled bodies without requiring any interface operation.

A particularly interesting feature of the Ausas FE space is that it is conforming with the standard FE one. Hence, no special treatment of the blending elements (the elements neighbouring the intersected ones) is required, thus leading to a purely local formulation. This becomes in a great advantage for the extension of the formulation to distributed memory platforms.

Aiming to describe the main geometrical features of the Ausas FE space, Fig. 5.1 depicts the same intersected element example that the original authors used in their proposal. By inspecting the Ausas shape functions that are obtained from the corresponding division in subelements, it can be observed that:

• the shape function values along an intersected edge are constant. Note that in the green side the shape function at the intersection points Q and P takes the same value as in A (Fig. 5.1b). On the other hand, in the red side the value at the intersection point P takes the same value as in B (Fig. 5.1c). The same can be observed for Q and C (Fig. 5.1d).
• in consequence, the shape functions values are constant in one of the split sides (the green one in Fig. 5.1a), and thus lead to an always constant interpolation with no gradient. Hence, only piecewise constant functions can be exactly reproduced.
• in the non-constant side, the shape function gradients are approximatively null in the intersection normal direction.

(a) Partition of ABC following the interface PQ.	(b) ${\displaystyle N_{A}}$
(c) ${\displaystyle N_{B}}$	(d) ${\displaystyle N_{C}}$
Figura 5.1: Ausas shape functions for the intersected finite element ABC (source [114]).

5.2.2 Nitsche imposition of the Navier-slip Boundary Condition

This subsection presents the general Navier-slip BC as well as the Nitsche method that is employed for its imposition. As it has been already mentioned, this technique is based on the work by Winter et. al. in [60], which is adapted to be used in combination with the aforementioned Ausas FE space.

In short, the Navier-slip BC can be understood as a sort of linear wall law whose behaviour is regulated by the slip length parameter ${\textstyle \varepsilon }$ as

where ${\textstyle \mathbf {u} ^{t}}$ and ${\textstyle x^{n}}$ denote the tangential wall velocity and the normal to the wall direction. As it can be observed, the Navier-slip BC approaches the completely stick condition as ${\textstyle \varepsilon \to 0}$ while it tends to the slip limit as ${\textstyle \varepsilon \to \infty }$.

From a numerical perspective, the general Navier-slip BC is defined as the conjunction of a no penetration constraint in the normal direction and a shear force imposition in the tangential one. While the first one is enforced on the normal projection of the velocity, the second is a boundary value imposition of the tangential projection of the viscous traction, meaning that the general Navier-slip condition is indeed a Robin-type BC. Accordingly, the Navier-slip BC can be split in a normal and a tangential contribution as

in which ${\textstyle \mathbf {g} }$ and ${\textstyle \mathbf {h} }$ are the velocity and the tangential traction to be imposed over the embedded boundary. ${\textstyle \mathbf {P} ^{n}}$ and ${\textstyle \mathbf {P} ^{t}}$ are the normal and tangential projection operators, which can be computed as ${\textstyle \mathbf {P} ^{n}=\mathbf {n} \otimes \mathbf {n} }$ and ${\textstyle \mathbf {P} ^{t}=\mathbf {I} -\mathbf {n} \otimes \mathbf {n} }$.

The imposition of the BC in Eq. 5.2 is achieved by using a stabilized Nitsche method in both the normal and tangential directions. The Nitsche imposition of the normal component reads

while the tangential one is

${\textstyle \gamma }$ is the penalty constant while ${\textstyle \zeta \in \{-1,1\}}$. The choice ${\textstyle \zeta =-1}$ leads to the adjoint inconsistent Nitsche formulation, which enjoys improved inf-sup stability for any value of ${\textstyle \gamma }$ [60]. Even though optimal convergence is not guaranteed for the velocity L²-error in this case, the adjoint inconsistent formulation is taken owing to its reported better stability properties. ${\textstyle \phi _{u}}$ is a stabilization constant defined as

It is important to remark that these two constraints need to be integrated in both sides of the level set interface after the ability of the Ausas FE space to completely disconnect the positive and negative sides of the intersection. Hence, the problem somehow translates to impose the same constraint in two different but overlapping interfaces.

In addition, it is also important to point out that the Ausas FE space reduces the approximation properties within the cut elements, so exactly where the Nitsche terms are active. This effect is particularly relevant since, depending on the splitting pattern, the gradient of the velocity and of the test functions is kinematically forced to be zero on one of the sides of the cut. In this context, the verification of the correct imposition of the Navier-slip BC is fundamental.

5.2.3 Drag calculation

Owing to the worse interpolation properties featured by the Ausas FE space, the boundary traction calculation might be inaccurate, specially under the presence of important wall shear effects. This section aims to address such undesired situation, which might limit the extension of the method to FSI problems.

In general terms, the computation of the embedded boundary traction can be performed by integrating the formula

over the level set intersections. Alternatively, it can be assumed that the Navier-Slip condition requires ${\textstyle t_{\parallel }=-{\frac {\mu }{\varepsilon }}\mathbf {P} ^{t}\left(\mathbf {u} -\mathbf {g} \right)}$, or equivalently that ${\textstyle \mathbf {P} ^{t}\left(\left(\mathbb {C} :{\boldsymbol {\nabla }}^{s}\mathbf {u} \right)\cdot \mathbf {n} \right)\approx -{\frac {\mu }{\varepsilon }}\mathbf {P} ^{t}\left(\mathbf {u} -\mathbf {g} \right)}$, and rewrite Eq. 5.6 in the form

These two expressions are equivalent in the continuum and when a “standard” FE formulation is employed. Not so however in the case of using Ausas shape functions, since in presence of zero velocity gradients would predict a zero skin friction in the first while the second would correctly estimate a skin friction provided that the tangential velocity is non-zero on the cut boundary (even when the gradient is zero in the immediate vicinity of the intersection).

In order to explain such inconsistency, the use of Ausas shape functions can be re-interpreted as the use of the standard FE space subjected to the additional constraint that the solution needs to be piecewise constant along each cut edge. The action of such geometric constraint is fundamental in the way equilibrium is established in the vicinity of the immersed objects. Indeed, the constraint allows forces to be transmitted from the cut boundary to the surrounding nodes without the need of producing a gradient in the solution, and is the reason for which the effective tangential stress can be non-zero in the vicinity of the embedded boundaries.

5.3 Validation

This section collects the benchmarking of the proposed formulation. To ensure that all the examples can be reproduced, the geometry, BCs and simulation settings are detailed prior to the results assessment.

In the following all the examples are listed with a brief description of their motivation

• 2D straight channel: this is a simple test case to check the implementation and convergence rates of the formulation considering three different wall behaviors (slip, intermediate and no-slip).
• 2D flow inside a ring: to obtain the convergence rates of the formulation when dealing with faceted approximations of curved geometries. The influence of the Nitsche constant ${\textstyle \gamma }$ is also evaluated with this example.
• 2D flow around pressurized cylinder: to prove the utility of the pseudo-compressible term to keep the pressure bounded when isolated fluid cavities appear.
• 2D elbow with internal wall: to test the performance of the formulation with a more complex geometry.
• 2D flow around thin-walled arch profile: to assess the performance of the proposed technique to solve high Re flows. This example is conceived as a 2D idealization of the final 3D example.
• 3D flow around a sailboat: this example presents a potential industrial application that consists in the CFD analysis of a sailboat.

5.3.1 2D straight channel

The objective of this very first example is to preliminary assess the correctness and convergence rates of the proposed formulation. Thus, the geometry and boundary conditions have been selected to be as simple as possible.

The geometry consists in a 2D 1.1x2m rectangular channel. The discontinuous level set function is set such that the effective width of the channel is 1m. Regarding the wall behaviour, the bottom edge is considered as a no-slip boundary while the embedded top one is tested with both slip and no-slip BCs. Moreover, an extra intermediate case in which the embedded skin ${\textstyle \varepsilon }$ equals 10^-1m is also solved.

The material properties are selected such that the Re number is 100. This is achieved by setting the density to 1kg/m³ and the viscosity to and 10^-2kg/ms.

A second order Backward Differentiation Formula (BDF2) scheme is used for the time discretization. To ensure that the steady state solution is reached, 20 time steps of 100s are solved.

The problem is solved for a set of structured meshes, whose characteristic element size and edge subdivisions are collected in Table 5.1.

5.3.1.1 Slip interface

In this case the flow is induced by applying a volume force of 0.5m/s² in the horizontal direction. The problem reaches a steady state that can be characterized by the rigid body movement

Table 5.2 collects the ${\textstyle L^{2}\left(\Omega \right)}$-norm of both the velocity and pressure errors. Fig. 5.2a and Fig. 5.2b depict the previous results. As it can be observed, the convergence order is ${\textstyle h^{2}}$ for the velocity field while it is superlinear for the pressure one, meaning that the obtained convergence rates match the linear FE ones even though the Ausas FE space is limited to a ${\textstyle h^{3/2}}$ convergence rate [114]. Such improved convergence is a particularity of the current test case, which is not affected by the null gradient limitation of the Ausas shape functions since the analytical solution gradient is zero in the top embedded boundary.

Finally, Fig. 5.3 presents the velocity field for three different mesh refinement levels (meshes 1, 3 and 5 in Table 5.1). The obtained solution matches the parabolic rigid body movement in Eq. 5.8.a in all cases. It can be also observed that, as expected, the gradient tends to zero in the region close to the embedded boundary as the mesh is refined.

(a) Velocity convergence.	(b) Pressure convergence.
Figura 5.2: 2D straight channel (slip interface). Solid lines represent the obtained results. Dashed lines represent ${\displaystyle h}$ and ${\displaystyle h^{2}}$ convergence rates.

(a) Mesh 1.	(b) Mesh 3.
(c) Mesh 5.
Figura 5.3: 2D straight channel (slip interface). Coarse, medium and fine meshes ${\displaystyle \lVert \mathbf {u} \rVert }$ field.

5.3.1.2 No-slip interface

Instead of using a body force to induce the flow, in this case a pressure gradient between the inlet and the outlet is applied. Hence, the pressure is fixed to 1Pa in the left edge and to 0Pa in the right one. This, in combination with the no-slip behaviour in both top and bottom boundaries, leads to the well known pressure driven flow between two steady plates benchmark, also known as Poiseuille flow. The steady state analytical solution of the Poiseuille flow can be particularised for the above mentioned material properties as

Table 5.3 collects the ${\textstyle L^{2}\left(\Omega \right)}$-norm of the velocity and pressure errors. The results exhibit a convergence rate around ${\textstyle h}$ for both velocity and pressure fields (Figs. 5.4a and 5.4b).

Fig. 5.5 depicts the obtained solution for three different mesh refinement levels (meshes 1, 3 and 5 in Table 5.1). The expected parabolic velocity profile can be noted in all cases. However, the peak value in the channel mid axis is far from the analytical solution value (6.25m/s) in the coarsest mesh. Regarding the pressure, the horizontal linear gradient can be observed in all cases.

Furthermore, it is worth to comment about the perturbations in the solution close to the embedded boundary, which can be clearly noted in the coarsest mesh pressure field (Fig. 5.4b). These are the direct consequence of the inability of the Ausas FE space to represent a gradient in one side of the intersected elements, which forces the solution in the embedded boundary vicinity to accommodate according to such extra constraint.

(a) Velocity convergence.	(b) Pressure convergence.
Figura 5.4: 2D straight channel (no-slip interface). Solid lines represent the obtained results. Dashed lines represent ${\displaystyle h}$ and ${\displaystyle h^{3/2}}$ convergence rates.

(a) Mesh 1 velocity norm.	(b) Mesh 1 pressure.
(c) Mesh 3 velocity norm.	(d) Mesh 3 pressure.
(e) Mesh 5 velocity norm.	(f) Mesh 5 pressure.
Figura 5.5: 2D straight channel (no-slip interface). Coarse, medium and fine meshes solution.

5.3.1.3 Slip length 10^-1 interface

This last variant of the straight channel test checks the performance of the proposed technique when it is applied in a wall law fashion. This means to introduce a controlled wall friction such that the embedded skin behavior is between the slip and no-slip limits. Hence, the slip length is set ${\textstyle \varepsilon }$ to 10^-1m in the embedded boundary. The inlet profile ${\textstyle v_{x}=1-(1/(1+\varepsilon ))y}$ is accordingly imposed in the left edge while the bottom edge velocity is fixed to [1,0]m/s.

These settings yield the analytical velocity and pressure steady state solution

which exerts a 0.0182N total horizontal drag force over the embedded boundary.

Table 5.4 collects the ${\textstyle L^{2}\left(\Omega \right)}$-norm of both the velocity and pressure errors as well as the ${\textstyle L^{2}\left(\Gamma \right)}$-norm of the horizontal drag force error. These results are presented for three different values of the penalty constant ${\textstyle \gamma }$. As it can be observed in Fig. 5.6, the convergence rate of the three analysed magnitudes is in all cases around ${\textstyle h}$, which is in line with the results previously obtained for the pure no-slip case. Besides that, it can be also noted that the convergence rate of the formulation is not affected by the ${\textstyle \gamma }$ value. As expected, the error constant minimizes for smaller ${\textstyle \gamma }$ values.

Complementary, Fig. 5.7 shows the obtained velocity fields for three mesh refinement levels (meshes 1, 3 and 5 in Table 5.1). The expected linear velocity profile can be clearly identified in all cases.

Finally, it has to be highlighted that, even though this test involves a very simple geometry, the obtained results are a very first, but essential, proof of convergence required prior to the application of the proposed technique to more realistic scenarios.

(a) Velocity convergence.	(b) Pressure convergence.
(c) Drag convergence.
Figura 5.6: 2D straight channel (Slip length 10-1 interface). Solid lines represent the obtained results. Dashed lines represent ${\displaystyle h}$ and ${\displaystyle h^{3/2}}$ convergence rates.

(a) Mesh 1.	(b) Mesh 3.
(c) Mesh 5.
Figura 5.7: 2D straight channel (Slip length 10-1 interface). Coarse, medium and fine meshes velocity field.

5.3.2 2D flow inside a ring

The aim of this test is to obtain the convergence rates of the presented formulation, in both the slip and no-slip limits, when dealing with polygonal approximations of curved boundaries. It is known that the imposition of a slip BC over faceted approximations of curved boundaries can be problematic [57]. This is explained by the fact that the slip BC requires the velocity to be parallel to each face of the discrete polygonal approximation. In consequence, the unique solution that fulfills such requirement in the vertices that are shared between faces is the null velocity. As it is observed in [57,37], this can lead to a wrong velocity approximation in the slip wall, specially in those cases that present an insufficient discretization of the curved boundary.

The geometry of the problem consists in two concentric cylinders separated by an interior fluid cavity. The inner cylinder radius is 0.5m while the outer one is 1.0m, meaning that the fluid domain ${\textstyle \Omega }$ can be described as ${\textstyle \left(x,y\right)\in \mathbb {R} ^{2}:0.25\leq x^{2}+y^{2}\leq 1}$. The fluid is initially at rest and a unit angular velocity ${\textstyle {\bar {v}}_{\Theta }}$ is imposed in the outer cylinder to induce the flow. The pressure is fixed to zero along the outer cylinder boundary to obtain a radially symmetric solution.

With regard to the material properties, the dynamic viscosity ${\textstyle \mu }$ and density ${\textstyle \rho }$ are 1.0kg/m³ and 10^-3kg/ms. A second order BDF2 scheme is used for the time discretization. To ensure that the steady state solution is reached, 10 time steps of 200s are solved.

A centered structured mesh is used in all the refinement levels of the study. Note that the outer cylinder is completely meshed disregarding the inner one, which is represented by a radial discontinuous distance function. The characteristic element size together with the number of radial and perimeter subdivisions are collected in Table 5.5. Fig. 5.8 shows the Mesh 1 refinement level in Table 5.5. More specifically, Fig. 5.8a depicts the centered structured pattern while Fig. 5.8b details the discontinuous level set function intersection pattern.

In the following, the convergence of the obtained solutions when a slip and no-slip BC is imposed over the inner boundary is studied by means of the analytical solutions, which can be easily obtained from the mass and momentum conservation equations written in polar coordinates.

(a) Mesh capture.	(b) Level set intersection pattern.
Figura 5.8: 2D flow inside a ring test case geometry (mesh 1).

5.3.2.1 Slip interface

In order to ensure a slip-like behavior, the slip length ${\textstyle \varepsilon }$ is set to 10⁸ in all the cases presented in this subsection. The Nitsche penalty constant ${\textstyle \gamma }$ is 0.1 in all the examples. The analytical solution can be expressed in polar coordinates as

Table 5.6 collects the ${\textstyle L^{2}\left(\Omega \right)}$-norm of the velocity and pressure errors. For both velocity and pressure fields, the results display a superlinear convergence rate close to ${\textstyle h^{3/2}}$ (Figs. 5.9a and 5.9b), which is in line with slip straight channel observations as well as with the results in [114,37].

The obtained velocity and pressure fields are shown for three of the mesh refinement levels (mesh 0, 3 and 6) in Fig. 5.11. At first sight, it can be easily observed that the coarsest mesh solution if far from the expected slip behavior. However, the method converges to the slip solution as the mesh is refined. These observations are in agreement with the ones in [57] about the resolution of a similar Stokes problem using a Nitsche imposition of the slip BC. The same issues are also reported in Chapter 4 when solving the same test but using the Navier-Stokes embedded discontinuous purely slip formulation. Having said this, the deterioration of the convergence rates to a superlinear rate below the expected ${\textstyle h^{3/2}}$ one is presumably associated to the polygonal approximation of curved boundaries rather than to the Nitsche imposition of the slip BC.

(a) Velocity convergence.	(b) Pressure convergence.
Figura 5.9: 2D flow inside a ring (slip interface). Solid lines represent the obtained results. Dashed lines represent ${\displaystyle h}$ and ${\displaystyle h^{3/2}}$ convergence rates.

Complementary, the sensitivity of the formulation to the value of the penalty constant ${\textstyle \gamma }$ when ${\textstyle \varepsilon \rightarrow \infty }$ is also assessed. To that purpose, the same example is solved but using the finest mesh and a set of ${\textstyle \gamma }$ values ranging from from 10^-4 to 10². Table 5.7 collects the velocity and pressure error ${\textstyle L^{2}\left(\Omega \right)}$-norms for all the ${\textstyle \gamma }$ values. These results are also depicted in Fig. 5.10. It can be observed that both error norms minimize for ${\textstyle \gamma }$ values between 0.1 and 10, being the minimum values around the unit ${\textstyle \gamma }$ value.

	Pressure error norm.
(a) Velocity error norm.	(b) Pressure error norm.
Figura 5.10: 2D flow inside a ring (slip interface). Error norms for different ${\displaystyle \gamma }$ values.

(a) Mesh 0 velocity.	(b) Mesh 0 pressure.
(c) Mesh 3 velocity.	(d) Mesh 3 pressure.
(e) Mesh 6 velocity.	(f) Mesh 6 pressure.
Figura 5.11: 2D flow inside a ring (slip interface). Coarse, medium and fine meshes solutions.

5.3.2.2 No-slip interface

In this case the slip length ${\textstyle \varepsilon }$ is set to 0 to enforce a no-slip wall behavior in the embedded interface. The Nitsche penalty constant ${\textstyle \gamma }$ is kept as 0.1. As in the slip case, the results and convergence assessment are based on the analytical solution, which in in polar coordinates reads as

Table 5.8 collects the velocity and pressure error ${\textstyle L^{2}\left(\Omega \right)}$-norms. These values are also depicted in Figs. 5.9a and 5.9b. It can be observed that both velocity and pressure fields exhibit a linear convergence rate. Such deteriorated convergence rate can be again associated to the inability of the Ausas FE space to properly capture the shear effects within the intersected elements.

Fig. 5.11 depicts the solution for three of the mesh refinement levels (meshes 0, 3 and 6). It can be observed that the solution tends to the analytical one as the mesh is refined.