UniversitÃ degli Studi di Roma "Tor Vergata"

Università degli Studi di Roma "Tor Vergata"

Mechanical Engineering Department Doctor of Philosophy Thesis in Mechanical Engineering Advanced FSI analysis within ANSYS FLUENT by means of a UDF implemented explicit large displacements FEM solver

Tutor:

Ph.D. Candidate:

Prof. Fabio Gori

Emiliano Costa

Co-Tutor: Prof. Marco Evangelos Biancolini

Academic Year 2011/2012 Ph.D. Program in Environment and Energy Cycle XXIV

This page has been left intentionally blank

To Valentina, who helps me in opening and closing the parentheses of life

Questo secolo oramai alla ne saturo di parassiti senza dignità mi spinge solo ad essere migliore con più volontà

F.B.

Abstract In the present thesis the development of an eective and ecient methodology to handle three-dimensional FSI studies involving solid isotropic thin-walled structures by means of the commercial CFD software ANSYS FLUENT is described. With the aim to extend the standard analysis capabilities of the computational tool accounting for the geometrically non-linear behaviour of deformable components during simulation, a C-based parallel algorithm founded on FEM explicit co-rotational formulation approach has been implemented through the UDF feature. Moreover, in order to fruitfully support the user throughout the set-up of the FSI study, an intuitive and easy-to-use GUI within the same workbench user environment of the CFD solver has been developed. Several test cases have been performed to validate the correctness and reliability of numerical outputs obtained using the implemented FEM solver. These results have been compared with either known solutions or with those gained through a consolidated commercial code. Finally, the employment of the proposed methodology for the numerical reproduction of FSI in real applications is presented, so as to demonstrate its actual potentialities, and candidate it as a mature and robust computational means to cope with ever-growing engineering design requirements.

Keywords: FSI, Fluid-Structure Interaction, FEM, Finite Element Method,

i

explicit method, co-rotational, large displacements, small strains, partitioned approach, CFD, Computational Fluid Dynamics, FVM, Finite Volume Method.

Correspondence to E-mail:

[email protected]

ii

Acknowledgements This doctoral thesis encloses the eorts not only mine, but also from people who enjoyed, at least I hope, sharing this not easy trip with me. First and foremost I want to thank Prof. Gino Bella for transferring me his a

priori optimism and for encouraging me to undergo this exciting experience which has made me renew my private engineering passion. I wish to express my particular gratitude and estimate to my co-tutor, Dr. Marco Evangelos Biancolini, whose expertise, open-minded and sparkling approach have considerably enhanced my Ph.D. experience. My sincere appreciation to Mr. Fabrizio Lagasco for his support since the beginning. Not the least of my thanks is to Ignazio Maria Viola for his availability to let me use the model for the second industrial test case, Tanya Scalia and Margaret Teale for reviewing this manuscript (actually I have not realized yet who was the murder too), Carlos de la Cuesta de Bedoya and Enrico Fiore for their translations of the abstract, and Ernesto Mottola for his really helpful suggestions (don't worry, any "Maccheronic linguistic crimes" are just mine). Emiliano Costa

...the wise old man living in Rocca di Papa says: "'Ffonna zì, sinnò t'a leva..."

iii

Contents Abstract

i

Acknowledgements

iii

Notation

xi

Acronyms and abbreviations used I Introduction

xvi 1

I.1

Objective of the document and scope of the work . . . . . . .

3

I.2

Structure of the document . . . . . . . . . . . . . . . . . . . .

4

II FSI statement

6

II.1 Eulerian versus Lagrangian formulation . . . . . . . . . . . . .

6

II.2 FSI physical phenomenon . . . . . . . . . . . . . . . . . . . .

7

II.3 Strategies for simulating FSI problems . . . . . . . . . . . . .

9

II.4 FSI scientic publications review . . . . . . . . . . . . . . . . 12

III ELDS statement

15

III.1 FEM approach to CSD . . . . . . . . . . . . . . . . . . . . . . 15 III.1.1 Implicit versus explicit time integration . . . . . . . . . 17 III.2 Introduction to two-dimensional FEM calculation . . . . . . . 21 III.2.1 Basic concepts of plate bending behaviour . . . . . . . 21 III.2.2 Two-dimensional stresses and strains . . . . . . . . . . 22

iv

CONTENTS III.2.3 Finite Element stiness matrix calculation . . . . . . . 25 III.2.4 Fundamentals of discretisation . . . . . . . . . . . . . . 26 III.3 ELDS implementation . . . . . . . . . . . . . . . . . . . . . . 33 III.3.1 Co-rotational formulation . . . . . . . . . . . . . . . . 35 III.3.2 ELDS solution characterization . . . . . . . . . . . . . 38 III.3.3 ELDS element stiness matrix formulation . . . . . . . 40 III.3.4 Explicit method nodes management . . . . . . . . . . . 57 III.4 Overall ELDS operations . . . . . . . . . . . . . . . . . . . . . 60

IV Coupling of ELDS with FLUENT

64

IV.1 Introduction to FLUENT suite . . . . . . . . . . . . . . . . . 64 IV.2 Standard CFD analysis procedure . . . . . . . . . . . . . . . . 65 IV.2.1 Pre-processing stage . . . . . . . . . . . . . . . . . . . 67 IV.2.2 Computing stage . . . . . . . . . . . . . . . . . . . . . 72 IV.2.3 Post-processing stage . . . . . . . . . . . . . . . . . . . 73 IV.3 Customization of FLUENT computing . . . . . . . . . . . . . 74 IV.3.1 FLUENT solver . . . . . . . . . . . . . . . . . . . . . . 74 IV.3.2 Solution customization by means of the UDF feature . 79 IV.3.3 UDFs parallel implementation . . . . . . . . . . . . . . 83 IV.3.4 Customized GUI implementation . . . . . . . . . . . . 86 IV.4 Dynamic mesh and update methods in FLUENT . . . . . . . 87 IV.5 Outline of overall methodology process . . . . . . . . . . . . . 89

V ELDS validation

91

V.1 Preliminary validation tests . . . . . . . . . . . . . . . . . . . 91 V.1.1 Free falling test . . . . . . . . . . . . . . . . . . . . . . 92 V.1.2 Initial conguration inuence verication . . . . . . . . 95 V.2 Comparison with LS-DYNA code . . . . . . . . . . . . . . . . 98

VI Industrial case studies

102

VI.1 Case study 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

CONTENTS

v

CONTENTS VI.1.1 Case study 1: model set-up . . . . . . . . . . . . . . . 103 VI.1.2 Case study 1: computational outputs . . . . . . . . . . 107 VI.2 Case study 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 VI.2.1 Case study 2: model set-up . . . . . . . . . . . . . . . 111 VI.2.2 Case study 2: computational outputs . . . . . . . . . . 113

VIIConclusions

118

Bibliography

120

Index

126

Appendix 1

132

Appendix 3

133

Appendix 2

135

CONTENTS

vi

List of Figures II.1 Methodological approaches for the FSI numerical simulation . 11 II.2 Principal characteristics of the FSI methodological approaches

12

III.1 Characteristic time scale of typical transient dynamics analyses 20 III.2 Generic rectangular plate representation . . . . . . . . . . . . 22 III.3 Stresses within a dierential element of a plate . . . . . . . . . 24 III.4 Pascal's triangle for the complete polynomial characterization

27

III.5 Denition of area coordinates χ1 , χ2 , and χ3 . . . . . . . . . . 29 III.6 Triangular element in the local system and area coordinates . 30 III.7 Transformation from the reference to actual time element conguration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 III.8 Building-up of ELDS general at triangular element . . . . . . 41 III.9 Constant Strain Triangle (CST) . . . . . . . . . . . . . . . . . 42 III.10Kirchho's assumptions . . . . . . . . . . . . . . . . . . . . . 47 III.11Mesh manipulation to nodal masses conguration achievement 57 III.12A generic node normal calculation . . . . . . . . . . . . . . . . 59 III.13Nodal rotation convention . . . . . . . . . . . . . . . . . . . . 59 III.14ELDS initialization operations . . . . . . . . . . . . . . . . . . 61 III.15ELDS computing operations . . . . . . . . . . . . . . . . . . . 62 IV.1 Flow chart of the standard CFD unsteady analysis

. . . . . . 66

IV.2 Geometrical model of the diuser . . . . . . . . . . . . . . . . 68 IV.3 Geometrical model of the reed valves housing

. . . . . . . . . 68

IV.4 Cleaned-up model of the second industrial test case . . . . . . 70 vii

LIST OF FIGURES IV.5 Surface mesh detail of the rst industrial test case . . . . . . . 71 IV.6 Volume mesh detail of the rst industrial test case . . . . . . . 72 IV.7 Sails pressure distribution in steady ow condition . . . . . . . 73 IV.8 Stationary ow path lines around sails . . . . . . . . . . . . . 74 IV.9 Flow chart of the segregated approach . . . . . . . . . . . . . 81 IV.10Flow chart of the density-based approach . . . . . . . . . . . . 81 IV.11Flow chart of a transient FSI study through ELDS introduction 82 IV.12Scheme of data transfer in FLUENT parallel environment . . . 84 IV.13FSI panel to set up the FEM parameters . . . . . . . . . . . . 87 V.1 Mesh of the sheet model of the rst preliminary validation trial 92 V.2 Mesh of the parts constituting the free falling test case . . . . 94 V.3 The shape assumed by the sheet during simulation

. . . . . . 94

V.4 Comparison between ELDS and theoretical falling test prole V.5 Model of the second validation test case

95

. . . . . . . . . . . . 96

V.6 Congurations of the sheet during simulation . . . . . . . . . . 97 V.7 Vertical coordinate prole of the monitoring node . . . . . . . 97 V.8 Monitoring point displacement z-component prole . . . . . . 99 V.9 Comparison between serial and parallel computing results . . . 100 V.10 Numerical proles comparison . . . . . . . . . . . . . . . . . . 100 V.11 Final deformed conguration of the sheet . . . . . . . . . . . . 101 VI.1 Examples of commercial reed valve pack . . . . . . . . . . . . 103 VI.2 Generation of the model of the convergent cone of the carburettor104 VI.3 Solid model of the primary inlet duct of the crankcase . . . . . 105 VI.4 Detail of the surface mesh of the basement . . . . . . . . . . . 106 VI.5 Discretisation of the upper three petals . . . . . . . . . . . . . 106 VI.6 Flow path lines for 10 kP a gauge pressure case

. . . . . . . . 108

VI.7 Opening of petals for 10 kP a gauge pressure case . . . . . . . 109 VI.8 Monitoring nodes vertical displacement prole . . . . . . . . . 110 VI.9 Parts constituting the simulation volume . . . . . . . . . . . . 112 LIST OF FIGURES

viii

LIST OF FIGURES VI.10Clamped nodes of main sail (left) and spinnaker (right) . . . . 113 VI.11The sails blowing down sequence

. . . . . . . . . . . . . . . . 114

VI.12The sails blowing sequence . . . . . . . . . . . . . . . . . . . . 115 VI.13Proles of total force components acting on main sail . . . . . 116 VI.14Proles of total force components acting on spinnaker . . . . . 116

LIST OF FIGURES

ix

List of Tables IV.1 Characteristic parameters for the proposed methodology . . . 89 V.1 Physical parameters values for both preliminary validation trials 92 V.2 Geometrical dimensions of the sheet model . . . . . . . . . . . 93 V.3 Parameters of the computational model . . . . . . . . . . . . . 98 VI.1 Physical parameters values for the industrial test case 1 . . . . 107

x

Notation What follows is a list of the symbols used in the text arranged as they appear in dierent chapters.

Mathematical symbols {} { }T ∂ ∂a ∂b

{˙} {¨} , [ ] [ ]T [R ]−1 R dV PdA

CAP. III x, y , z î, ˆj, k ˆ [M] [C] [K] ¨ {d} ˙ {d}

column vector column vector tranpose partial partial derivative of generic parameter a with respect to generic parameter b time dierentiation of column vector double time dierentiation of column vector derivation rectangular or square matrix rectangular or square matrix transpose square matrix inverse integral over the volume integral over the area summation (very) strict inequality, is much less than

axes of the global Cartesian system unit vectors along x, y , and z axes system mass matrix system damping matrix system stiness matrix nodal acceleration vector nodal velocity vector

xi

{d} {Fe } {r} {r}tn {r}tn+1 t b, c u, v , w [E] E ν ρ {σ} σx , σy , σz τxy {γ} x , y , z γxy Ue [k] [B] [D] D {Fnodal } {Fel } {Fext } {Fgrav } {Ffluid } Mx , My Mxy {κ} κx , κy κxy ξ, η

Notation

nodal displacement vector system force vector nodal position vector nodal position vector at tn nodal position vector at tn+1 time, thickness of a plate largest dimensions of a plate components of displacement along x, y , and z axes material matrix modulus of elasticity of the material Poisson's ratio of the material density of the material stresses vector normal component of stress parallel to x, y , and z axes shearing stress component in rectangular coordinates strains vector unit elongation along x, y , and z directions shearing strain component in rectangular coordinates strain energy stiness matrix strain-displacement matrix modulus matrix exural rigidity of the plate nodal force vector nodal elastic force vector nodal external force vector nodal force vector due to gravity nodal ow force vector bending moment per unit length of sections of plate perpendicular to x and y axes twisting moment per unit length of sections of plate perpendicular to x axis curvature vector curvature of the middle-surface of the plane in xz and yz planes twist of the surface with respect to the x and y axes isoparametric coordinates xii

s lij lmin V am S tn tn+1 tref ∆t ∆tcrit ∆tELDS ∆tCF D nF EM [k]CST [k]DKT [k]DRILL [k]ELDS {a} ai {N} Ni Pn (x) Pn (x, y) P (x, y) (xi , yi ) A A 1 , A2 , A 3 χ1 , χ2 , χ3 L23 h φ(x) {φ} φi φe (x) {X} [A] [N] C ref

Notation

edge-tangent coordinate the length of side ij of a triangle smallest distance between two nodes of the surface mesh of the deformable components element volume speed of sound in the material time step safety factor time at n iteration time at n + 1 iteration, actual time reference time time step size critical time step size time step size adopted by ELDS time step size for CFD transient solution number of FEM subtimesteps of calculation CST element stiness matrix DKT element stiness matrix DRILL element stiness matrix ELDS element stiness matrix generalized coordinates vector ith generalized coordinate shape functions vector ith shape function complete one-dimensional polynomial of order n complete two-dimensional polynomial of order n generic point of the triangle coordinates of ith node of the element triangle area triangle subareas area coordinates in two-dimensions length of the triangular element side connecting vertexes 2 and 3 distance from vertex 1 to side 23 of the triangular element generic eld variable generic eld variable vector ith component of the eld variable vector {φ} nodal values of the eld variable vector {φ} vector of terms of the variable of the polynomial matrix of nodes coordinates shape functions matrix element conguration at initial time t=0 xiii

C tn+1 θx (x, y), θy (x, y) θz (x, y) βx (x, y), βy (x, y) α

CAP. IV x, y , z î, ˆj, k ˆ ~v vx , vy , vz vi ~g p ρ µ ν ~ ext F Fext,i hi Ei T Sm Sh kef f k kt J~j τ¯ τij ~ ∇ I D Dt

Notation

element conguration at actual time tn+1 rotations of the normal to the undeformed surface in the xz and yz plane ctitious out-of-plane rotation of the normal rotations of the normal to the undeformed surface in the xz and yz plane coecient of DRILL matrix

axes of the global reference Cartesian system unit vectors along x, y , and z axes velocity vector components of velocity vector ~v along x, y , and z axes ith component of velocity vector ~v acceleration of gravity static pressure density dynamic viscosity cinematic viscosity external body vector ~ ext ith component of external body vector F enthalpy internal energy temperature source term source heat of chemical reaction eective thermal conductivity thermal conductivity turbulent thermal conductivity diusion of species j stress tensor ij component of stress tensor τ¯ grad operator unit tensor total derivative

xiv

CAP. V x, y , z t g Fx , Fy , Fz

Notation

axes of the global reference Cartesian system time acceleration of gravity x, y, and z component of the total force acting on the sail

xv

Acronyms and abbreviations used CAD CAE CFD CFL CSD CST DKT DOF DPM e.g. ELDS etc FEA FEM FSI FVM GUI ID IGES NURBS PDE(s) RBF(s) SI STEP STL TL TUI UDF(s) UL VRML

Computer-Aided Design Computer-Aided Engineering Computational Fluid Dynamics Courant-Friedrichs-Lewy Computational Structural Dynamics Constant Strain Triangle Discrete Kirchho Theory (Triangle) Degrees Of Freedom Discrete Phase Model exempli gratia Explicit Large Displacements Solver et cetera Finite Element Analysis Finite Element Method Fluid-Structure Interaction Finite Volume Method Graphical User Interface Identication Data Initial Graphics Exchange Specication Non-Uniform Rational B-Spline Partial Dierential Equation(s) Radial Basis Function(s) International System (of Units) Standard for the Exchange of Product Model Data Stereolithography Total Lagrangian Text User Interface User Dened Function(s) Updated Lagrangian Virtual Reality Markup Language

xvi

Chapter I Introduction Modern engineering and scientic design are more and more driven by a number of evolving operational requirements such as increasingly stringent manufacturing constraints, reduction of costs, and minimization of time-tomarket. As a consequence, all these factors force the designer to operate with a ne balance between the simplicity of approach to be adopted in modelling and the high degree of accuracy to be provided in computational results. In this context, the use of CAE (Computer-Aided Engineering) solutions and related simulation softwares, even more if used since the early stage of the design process, may considerably improve the product reliability more inexpensively and in a shorter time than building and testing physical prototypes. While the simulation of most relevant problems was still far out of reach a couple of decades ago, this possibility is now available thanks to several decisive elements such as the advance in computational power, numerical modelling approaches and methods as well as the considerable decrease of hardware-connected costs. Hence, CAE techniques like CSD (Computational Structural Dynamics) and CFD (Computational Fluid Dynamics) tools for instance, have become extensively utilized and universally accepted as an indispensable means in conceptual design due to their high level of delity in complex systems behaviour and related physical phenomena prediction. This type of numerical technologies is a fairly recent discipline crossing the

1

boundaries of lots of science branches, and it is widely used in all engineering sectors by industrial companies, ordnance, government and corporate laboratories, as well as research organizations and universities. However, it is commonly known that, in order to eectively face non standard physical phenomena and manage complex geometries through the current advanced tools, highly experienced users are strictly requested. As a matter of fact, to prevent the delivering of incorrect evaluations upon which critical decisions will be based, the designer shall be conscious of both the strengths and shortcomings of these computer-based instruments and, in order to attain this professional awareness and the adequate expertise, he shall be suitably trained. Once able to properly manage and exploit at best CAE analyses capabilities and gain an insight into the physics of the phenomena involved, he can simulate the behaviour of complex systems and, thus, be productive by signicantly reducing iterative trials and errors with the consequent saving of time and money. Considering these very attractive advantages, in the last decades plenty of both commercial and open source computational tools have been developed and are insistently proposed in all scientic elds. Part of these softwares aim at solving just a very specic application, whereas others of most recent implementation are designed and built-up to cope with multi-physics problems. In particular, the multi-physics topic has recently caught much attention because the coupled solution of dierent physical problems is required to suitably perform the numerical simulation of various processes in engineering as well as in nature. Given that, the use of a tool designed for satisfying a specic physical task turns out to be inadequate when dealing with analyses involving more than one physical eld. Consequently, to have to recourse to a third party code is often a design necessity. In this case, two psychological attitudes shall be properly considered and not be underestimated. The rst is the user disorientation in choosing the

CHAPTER I.

Introduction

2

I.1

Objective of the document and scope of the work

most eective procedure according to its own design needs and the second, particularly when professional designers are involved, is the predisposition to neither invest eorts in training on new tools nor deepen the knowledge in technical contexts that are far from its own specialization. This is typically the scenario arising for the FSI (Fluid-Structure Interaction) numerical simulation because it implies, at least, the resolution of both structural and uid-dynamic elds and, when temperature eects are relevant, the thermal one as well.

I.1 Objective of the document and scope of the work The aim of the present document, constituting the Doctor of Philosophy thesis, is reporting the numerical implementation and the practical use of a three-dimensional structural solver to be dynamically linked to ANSYS FLUENT (hereafter referred to just as FLUENT) commercial CFD suite in order to enable the FSI analysis of isotropic thin-walled components. This solver, which is referred to as ELDS (Explicit Large Displacements Solver) from now on in the document, is based on the explicit FEM (Finite Element Method) approach and has the main scope to provide the CFD designer with a robust and reliable computational means to face the FSI task in real technical applications also when large displacements of deformable structures occur. The activities required to achieve the objective of the present work have been developed in collaboration with the FSI Project Team at Mechanical Engineering Department of the University of Rome "Tor Vergata".

CHAPTER I.

Introduction

3

I.2

Structure of the document

I.2 Structure of the document The document is arranged in such a manner to ease the comprehension of the modalities according to which ELDS has been developed and can be utilized to full the FSI design requirements of real engineering applications. The theoretical background on which ELDS is founded, the strategical procedure for its embedding into FLUENT frame as well as how its introduction modies the standard CFD solution process are detailed. Because of the strong industrial mark of the ELDS project, some of the principal topics are just treated at high level reporting the fundamental guidelines, with the intention to focus the attention on its practical use to solve real and complex problems. To this end, an entire chapter is devoted to the computational simulation of industrial test cases in view of showing the eectiveness of its utilization and its robustness in managing them. To oer a rst impression of the matters concerning this thesis, a brief summary of the contents of each chapter is presented as follows. After an introduction of the scientic background on which this thesis is in-

chapter I delineates the structure and contents of the document. Chapter II denes the physical phenomenon of FSI evidencing its importended,

tance in modern engineering and technical applied sciences, and the strategies to simulate it. The last part of the chapter contains a review of the international scientic works devoted to the methodologies for its simulation.

Chapter III reports, after an excursus on modern techniques adopted to face structural dynamic problems according to the FEM theoretical approach and an introduction to the two-dimensional FEM calculation, the description of the ELDS implementation. The theoretical background on which ELDS is founded and the methodology employed to derive the structural large displacements solution are identied in detail. The matter of

chapter IV

is the ELDS use. Specically, after a brief in-

troduction on the FLUENT commercial suite, this chapter highlights the

CHAPTER I.

Introduction

4

I.2

Structure of the document

methodological procedure adopted to dynamically link ELDS to FLUENT solver. Particular attention is reserved to the description of the implemented algorithm and the procedure adopted to customize the numerical simulation by means of the feature provided by the CFD solver to enable the FSI analysis.

Chapter V

is devoted to the description of the validation tests for deter-

mining the correctness of the implemented FEM solver. Within

chapter VI

numerical examples of the use of ELDS to solve real

cases of engineering interest are given. Particular emphasis is put on the complete presentation of all modelling and discretisation details in order to evidence the eectiveness in the accomplishment of the FSI task. It is the scope of

chapter VII

to oer an overall summary of the work,

reporting the main conclusions, achieved goals, and future desirable developments and improvements.

CHAPTER I.

Introduction

5

Chapter II FSI statement In this chapter the FSI topic from dierent perspectives is presented in order to provide a general overview of this very challenging and complex issue. In particular, in section II.1 the dierence between the Eulerian and Lagrangian modelling approach is introduced. In sections II.2 and II.3 the denition of the physics of the phenomenon as well as its relevance in industrial engineering and applied sciences, and the numerical methodologies to approach it are respectively described. Finally, in section II.4 a review of the international scientic publications concerning this subject is reported.

II.1 Eulerian versus Lagrangian formulation Referring to CAE methodologies and modelling techniques of direct interest for the thesis, in the past decades many ecient methods for the numerical resolution of various problems in both the CFD and CSD eld have been developed and, thanks to this headway, lots of commercial tools have been created and successfully applied. In particular, as far as concerns CFD, the governing equations have been traditionally written by using Eulerian (spatial) coordinates [5]. From the numerical point of view, the FVM (Finite Volume Method) discretisation has been preferred to the others because of its conservative properties.

6

II.2

FSI physical phenomenon

On the other hand, CSD has also achieved a great advance independently from CFD. Numbers of structural dynamics codes have been developed to solve various structural dynamics tasks. The modelling of a wide range of material laws and structural properties has been made possible by the creation of special nite elements holding desired features. Contrarily to the uid dynamics, Lagrangian (material) coordinates [30] have been selected for the description of the governing equations. As successively stated in paragraphs II.3, this dierence of formulation in treating uids and solids is one of the topics that makes the FSI a very dicult problem to be virtually reproduced by means of the numerical calculation.

II.2 FSI physical phenomenon FSI is the interaction of movable or deformable structures with an internal or a surrounding uid ow. Since this physical phenomenon involves both structural mechanics and uid dynamics eld (in some cases thermodynamics as well), it is a multi-eld system implying that at least two dierent physical subtasks, structural and uid, are needed to be tackled and solved. With regard to the interaction between uid and structural domains, two types of mechanisms may arise. The rst takes place at the domains interface and is due to the uid pressure and shear stress, whereas the second one concerns the eects induced by temperature gradients. In the rst case, the resulting forces of uid ow cause structural domain deformation that, in turn, leads to a change in the ow eld. Therefore, the solution of each subtask has to be considered as a boundary condition on the interface for the other one. In the second case, to properly account for the interaction due to temperature eects, the simultaneous solution of the uid dynamic and energy equation in the ow domain as well as the coupled solution of structural dynamic and temperature-induced structural deformation in structural

CHAPTER II.

FSI statement

7

II.2

FSI physical phenomenon

domain are requested. These solutions shall then be coupled at domains' interface. As such, a generic conguration of an FSI analysis may foresee the presence of either only one of these mechanisms or both simultaneously. Specically, the structural components deformation due to uid loads rather than the uid ow control by modifying wall component shape conguration in the system has received much attention in recent years, and its relevance is still continuously growing because of the large number of practical cases of interest present in modern engineering and other disciplines of applied sciences. As a matter of fact, some of the technical elds and applied disciplines in which FSI can nd application are:

•

Aerospace/Aeronautical for the jet deformations, parachute opening, wing and rudder utter, turbine blade fatigue, rotary wing analysis;

•

Civil Engineering for the wind-induced oscillations of high buildings and long bridges, tent-roofs;

•

Oil&Gas/Power Engineering

for the wind energy, earthquake re-

sponse analysis of liquid storage tanks, pipeline security evaluation;

•

Automotive for the membrane pumps, heat exchangers, pipe-systems, stirring, bueting techniques, turbomachinery, airbags deployment, jet engines, read valves, aeroacoustics, aquaplaning;

•

Medical/Bio

for the analysis of aneurysms, modelling of the arti-

cial heart valves, inner ear motion, blood circulation in human body, pulmonary air ow modelling;

•

Defence for the aircraft fuel tank sloshing, rocket engine components modelling;

•

Marine/Nautical for the vibrations due to blast wave, yacht engineering, water penetration of o-shore structures, submarines deformation

CHAPTER II.

FSI statement

8

II.3

Strategies for simulating FSI problems

during missions, vessels slamming in extreme seas, underwater cable behaviour simulation. Moreover, the impressive number of scientic conferences and multi-eld modules that have been incorporated into the commercial software packages reect the practical importance of the matter, and give further evidence of the high interest on FSI. In the following sections a review of the numerical methodologies to simulate the FSI and international scientic publications produced on this topic are respectively reported.

II.3 Strategies for simulating FSI problems The interaction of a uid with a solid body is a widespread phenomenon in nature occurring at dierent length scales. From the engineering design point of view, to solve some technical applications the structural response to the uid can be neglected and, as such, only the uid dynamic subtask can be solved. In a similar manner, in some other cases the uid forces can be neglected with respect to other external forces and hence it is not necessary to account for ow behaviour response. However, in many physical processes and relevant engineering problems neither uid forces nor structural deformations can be neglected and, consequently, both subtasks must be carried out with an adequate level of detail. According to this latter scenario, the ow path depends on the uid-structure interface deformations and, in turn, the structural displacements are determined by the uid forces at the interface. Since just the analytical solution of uid-dynamic problem is not feasible in nearly all cases, to solve the FSI problem the only possibility is to employ a numerical approach. The detailed coupled numerical reproduction of this interaction turns out to be very challenging mainly for two reasons. At rst, the coupling of uid and structural subtasks leads to a highly non-linear behaviour of the whole system.

CHAPTER II.

FSI statement

9

II.3


Moreover, uid ows are typically given in Eulerian coordinates whereas the structure is treated in a Lagrangian framework (see section II.1). Given that, a dierent high-level modelling shall be carefully managed. The approaches to accomplish the FSI numerical solution can be grouped into two main types of strategy [29] as follows: 1. monolithic methods; 2. partitioned methods. The monolithic (or direct) methods foresee the aggregation of the governing equations of uid and structural subtasks and their simultaneous solution by means of a unique frame. According to this approach, the interface conditions on continuity of velocity and normal stresses become implicit (see section III.1.1) and are automatically satised. Opposite to monolithic procedure, the partitioned methods divide the coupled problem into uid and structural parts which are then solved separately. The term partitioning identies the process of spatial separation of a discrete model into interacting components (partitions) driven by physical or computational considerations. Referring to these methods, the crucial point is the need to transfer information from ow to structural eld and vice versa. Considering this data exchange, they can be divided into weakly (loosely) and strongly (implicit) partitioned algorithms depending on, respectively, whether data exchange at the interface is done only once or more times per solution time step. Figure II.1 illustrates the afore-mentioned classication. Both approaches have advantages and disadvantages. In particular, since the uid and the structural dynamic subproblems are based on dierent types of partial dierential equations, the main advantage of the weak coupling approach is that it allows the use of already developed, ecient, and well validated solvers for each single subtask. Depending on the generality of the two codes, arbitrary complex ows and structures can be considered and

CHAPTER II.

FSI statement

10

II.3


Figure II.1: Methodological approaches for the FSI numerical simulation successfully modelled. The only programming eort lies in creating suitable subroutines for information exchange between solvers. Unfortunately, due to the explicit nature of this coupling (see section III.1.1), convergence problems may arise with the consequent restriction on the choice of the time step (even if implicit time-stepping schemes are used by both solvers) that leads, in turn, to a restriction of the range of approachable applications. On the other hand, the monolithic coupling algorithms are more dicult to build-up and implement. The simultaneous solution of the whole FSI problem normally requires the reformulation of the systems of equations, and sets restrictions on the choice of the numerical methods to be applied. Additionally, special strategies may be needed for modelling the non-linearities in each of the physical domains. Huge programming eorts are usually required to create and validate a new algorithm applicable to various problems. However, because of the simultaneous solution of both the FSI subtasks, there are neither approximation errors nor convergence problems due to data transfer between uid and structural domain. In gure II.2 monolithic and partitioned approaches are compared with regard to their stability, generality as well as programming implementation eorts. Apart from the general subdivision on numerical approaches, there is a

CHAPTER II.

FSI statement

11

II.4

FSI scientic publications review

Figure II.2: Principal characteristics of the FSI methodological approaches widespread variety of themes for the theoretical development and applications related to FSI. In fact, the FSI perspectives may vary depending on types of ow elds covered (compressible, incompressible, laminar, turbulent, etc), types of applications, structural elds (thin-walled, rigid bodies, non-linear material, etc), discretisation schemes (nite volume, spectral methods, multibody dynamics), ow modelling assumptions (continuum, statistical Lattice Boltzmann distribution), and calculation grid treatment (moving grid, xed grid, immersed boundary). In the following section a concise review of the historical development and international scientic works on the numerical FSI is reported.

II.4 FSI scientic publications review In order to provide an overview of the historical development of the FSI numerical simulation, a review of the international scientic contributions was conducted. Due to the vastness of the FSI-related works, only those subjectively retained relevant are reported. As far as partitioned methods are concerned, these were introduced by Park and Felippa [38], and the key idea for these methods is clearly described by

CHAPTER II.

FSI statement

12

II.4


Felippa [24]. Partitioned solutions with staggered coupling algorithms were developed by Farhat [23] to handle aeroelastic wing problems. Strong coupling of partitioned algorithms was applied to large displacements 2D structural problems coupled to viscous incompressible uids by Wall [48] and by Ramm [49]. They also applied the same method for a coupled uid-structure environment with an initially at three-dimensional shell model as given in Wall and Ramm [50]. Other large displacements structural problems dealing with incompressible uids were detailed in Mok and Wall [36], [51], and Tallec and Mouro [42]. FSI with large displacements applied to wind problems was developed by Rossi [39], Wüchner and Bletzinger [55], Badia [3] and Wüchner [56]. More sophisticated developments on strong partitioned method for FSI problems can be found in Steindorf [41], Matthies and Steindorf [34], and Tezduyar [43]. A study on strong coupling partitioned methods for FSI applied to hemodynamic problems can be found in Nobile and Vergara [37], Causin [18], and Fernández and Moubachir [25]. Strong coupling of uid-structure interaction including free surfaces was studied in Dettmer [22] and Wall [52]. Recently, a new approach based on Robin transmissions conditions for uid-structure interaction problems was given in Badia et al. [4]. Referring to monolithic methods, examples can be found in Rugonyi and Bathe [40], and Heil [28], and references therein. Figueroa [26] proposed a simplied monolithic FSI algorithm by embedding the structure into the uid problem. There, the (d-1)-dimensional structure is modelled as a membrane. The same idea of writing the FSI problem only in terms of uid unknowns was presented by Nobile and Vergara [37], where an algebraic law for approximating the structure problem was employed. In any case, the use of non-modular preconditioners for the FSI system has received much less attention. The rst reason is the fact that they are not needed in applications with a negligible added-mass eect because partitioned procedures are very ecient. The second reason is the loss of modularity. Existing uid and

CHAPTER II.

FSI statement

13

II.4


structure codes can still be reused, but the coupling of the codes is more involved than bare interface communication. In fact, uid and structure matrices must be stored in a unique FSI matrix, which has to be accessed to compute the preconditioner. Concerning new research trends, a promising approach that has recently caught the attention of the scientic community is a meshing morphing technique based on RBF (Radial Basis Function) theory [12][13]. This numerical procedure represents a fast and meshless method enabling to tackle FSI by considering a set of deformed congurations of deformable components of the system preliminary calculated by performing a FEM modal analysis under the assumption of linear elasticity.

CHAPTER II.

FSI statement

14

Chapter III ELDS statement This chapter principally pertains to the ELDS implementation. After an excursus on FEM approach to structural dynamic problems, in section III.1 the description of the main advantages and disadvantages of both the implicit and explicit methods for time integration is provided. Then, section III.2 and III.3 respectively report an introduction of the two-dimensional FEM calculation, and fundamental concepts and passages concerning the ELDS algorithm implementation. Finally, in section III.4, an overview of the operations characterizing the ELDS use is presented. As just the basic theory of the explicit nite element method is discussed, an adequate explanation on the specic topics related to this numerical approach can be found in several specialized and known books [5] [19] [33] [57].

III.1 FEM approach to CSD FEM is a numerical technique used for solving the ordinary and partial differential equations that arise in engineering and mathematical physics, for which analytical solution is generally not possible because of the system complexity. This method, based on energy variational principle and discrete interpolation, has been developed since the middle of the 20th century and is applicable

15

III.1

FEM approach to CSD

in several elds of applied science such as frame analysis, stress analysis of solid continua, uid mechanics, heat ow, acoustics, and electromagnetics to name a few. Given this multi-disciplinary analysis capability, the quantities of interest can be diverse according to the specic analyzed system. In the case of structural analysis, context in which FEM constitutes surely the dominant solution procedure, the primary elds of interest of the system are displacements and stresses. The main typologies of structural problems which can be numerically investigated by means of FEM may be classied [33] as follows:

• equilibrium problems (steady state or static); • diusion problems; • transient dynamic (inertial) problems; • Eigenproblems analysis of vibration and stability. Due to the dynamic time-dependent nature of FSI phenomenon, from the structural point of view, the scientic interest in this thesis focuses on the use of FEM for transient dynamic problems solution, namely the determination of the dynamic response of deformable structures under the action of general time-dependent uid loads in incompressible subsonic conditions. This evaluation of structural induced deformations becomes particularly important regardless of engineering design, in the case that vortex shedding matches the resonance frequency of the structure. In a transient dynamic analysis, the distribution of displacements changes according to an acceleration law (Newton's motion law), and the behaviour of the structural system at a specic instant of time t can thus be obtained by solving the following matrix equation

¨ t + [C]{d} ˙ t + [K]{d}t = {Fe }t [M]{d}

CHAPTER III.

ELDS statement

(III.1)

16

III.1

FEM approach to CSD

where [M], [C], and [K], are respectively the mass, damping, and stiness ¨ t , {d} ˙ t , and {d}t are respectively matrices for the complete system and {d} the acceleration, velocity, and displacement vectors at a specic time t comprising values at a number of points (nodes of the mesh) in the system. Finally {Fe }t is a vector of time-dependent loads applied at these nodes. Marching forward in time, the unknown quantities in the equations can be solved by using either explicit or implicit time integration techniques (see section III.1.1). Due to the morphology of the components (thin-walled components) of interest for the present work, another problem that needs to be coped with is the occurrence of non-linear behaviour of the system. In general, the non-linear structural behaviour arises from a number of causes which can be grouped into three principal categories as follows [5]:

• material non-linearities; • geometric non-linearities; • changing status (including contact). When one or more of these non-linearities occurs during structural mechanics calculation, either the nature of the constitutive relationships (Hooke's law) is altered, or the terms in the governing strain-displacement relationships are made non-linear, or dierential equations are made non signicant. In particular, ELDS has been conceived and built to handle a transient dynamic problem even when large displacements of the conguration of study take place (geometric non-linearity) with the fundamental assumption of small strains.

III.1.1 Implicit versus explicit time integration The time integration schemes are numerical iteration techniques employed to obtain a solution of an evolving system at a specic time tn+1 , corresponding

CHAPTER III.

ELDS statement

17

III.1

FEM approach to CSD

to n + 1 steps of calculation subsequent to tn , that is the time at which the solution is known. These computational approaches are commonly classied as implicit and explicit schemes. Referring to a structural dynamic analysis, the implicit schemes are based on nding the position vector {r}tn+1 by using both the known value {r}tn and the value {r}tn+1 that cannot be expressed in closed form, where {r}tn and {r}tn+1 are the position vector {r} at time

tn and tn+1 respectively. For the manner in which this method has been conceived, the time integration of the discrete formulation of the governing relations requires the solution of algebraic equations. In particular, the achievement of the global equilibrium at each time step of calculation turns out to be computationally intensive because it involves the matrix factorization. Providing that equilibrium can be achieved, there is no theoretical limit to the size of the time step that can be used and, therefore, the implicit schemes are termed unconditionally stable. However, the accuracy of the implicit schemes deteriorates as the time step size increases relatively to the period of response of the system. Furthermore, considering the specic numerical needs of implicit solvers, even if ecient parallel algorithms for sparse matrices are currently present, the dimension of the system to be solved is still limited because the building-up and the inversion of matrices for high degrees of freedom (DOF) problems (millions of elements) are very dicult and often not feasible. On the other hand, the explicit schemes mainly aim at nding the position vector {r}tn+1 by using only the known value {r}tn . Since {r} at time tn+1 is dened as {r}tn+1 = {r}tn + {d}tn+1 , the only parameter to evaluate reduces to the displacements vector {d} at time tn+1 . As such, in the explicit methods the time integration of the discrete equations does not require the solution of any equations and, consequently, these local variables are directly evaluated without the need for global equilibrium calculation. Relating to the management and dimension of the system to be solved, these schemes have no limits for DOF handling and turn out to be expecially appropriate

CHAPTER III.

ELDS statement

18

III.1

FEM approach to CSD

for parallel implementation and computing. However, in order to guarantee numerical stability, these techniques need to be performed at small increments and they are consequently termed conditionally stable. In particular, as hereafter described in the document (see paragraph III.3.2), the calculation time step size ∆t must be less or equal to the critical time step size ∆tcrit for all time steps. The critical time step is related to the smallest natural period of the discretised studied system that, in turn, essentially depends upon element smallest size and the wave speed for the material element. The majority of the mainstream and traditional FEM codes use the implicit scheme procedure for time integration. Analyses that are well-suited to implicit solution techniques are for instance the static, low speed dynamic, and steady-state transport studies. By contrast, the simulation performed through the explicit-based codes are typically characterized by a time step even several orders of magnitude smaller than that used in an implicit scheme. The explicit codes are generally better for short transient events where the eects of stress waves are important and the time step has to be mandatory small for some reasons such as the accuracy. Problems theoretically well-suited to explicit solution are those involving, for example, the material non-linearity, large geometric non-linearity, non-linearity in combination with large displacements, simultaneous large displacements and contact problems, rapidly changing or discontinuous loading, and rigid body motion. Considering what afore stated, the time scale then constitutes the reference parameter for a preliminary indication on the time integration scheme advisable to adopt. Figure III.1 depicts some typical transient dynamic applications approachable by means of the FEM simulation and their corresponding characteristic time scale. Since an absolute dividing line about applications that can be approached through one method rather than another is not traceable, which integration scheme is desirable to use often depends upon the specicity of the problem and the available computing power. Focusing on the matter of the present work, typical applications which are

CHAPTER III.

ELDS statement

19

III.1

FEM approach to CSD

Figure III.1: Characteristic time scale of typical transient dynamics analyses

well-suited for the explicit analysis include the following contexts:

•

Defence

for the impact/penetration, armor and anti-armor systems,

kinetic energy and chemical energy devices, underwater shock and explosions, weapon design, blast response;

•

Automotive for the crashworthiness analyses, durability, crash, pedestrian protection, interior safety, restraint systems design, airbag design, sheetmetal stamping, rollover protection;

•

Aeronautics/Aerospace

for the bird strike, material forming, im-

pact, explosion and shock loadings, crack, containment;

•

Petrochemical for the gas and dust explosions, accident simulation, uid sloshing, wall perforation;

•

Nuclear for the de-commissioning, pipe break and whip, jet and missile impingement;

•

Transport for the explosions in vehicle and tunnels, crashworthiness, occupant dynamics, safety;

•

Manufacturing for the stamping, forging;

CHAPTER III.

ELDS statement

20

III.2

Introduction to two-dimensional FEM calculation

•

Structural for the earthquake safety, concrete structures reliability;

•

Electronics

for the drop analysis, package design, thermal manage-

ment. Considering the above-listed applications, it turns out to be evident that the use of explicit schemes can absolutely be considered as a transversal multidisciplinary industrial problem resolution technique.

III.2 Introduction to two-dimensional FEM calculation This paragraph deals with some fundamental concepts of the FEM simulation of plates behaviour. Specically, the plate bending behaviour, twodimensional stresses and strains analysis, stiness matrix calculation, and fundamentals of discretisation are provided with the aim to ease the comprehension of the procedure adopted for the formulation of the ELDS triangular element described in section III.3.3.

III.2.1 Basic concepts of plate bending behaviour A plate is at body with a small thickness compared to its planar dimensions and constitutes one of the most important structural components. Figure III.2 shows, for instance, a representation of a generic rectangular plate having the largest dimensions b and c lying in the xy plane (in-plane directions), the thickness of length t measured in the out-of-plane direction z , the widest surfaces at z = ±0.5t, and the midsurface at z =0. According to the commonly used distinction, a plate is referred to as thin if its thickness is less than about one-tenth its in-plane dimensions b and c, that is t b or c, otherwise it is dened as thick. Relating to structural behaviour, it supports loads transverse or perpendicular to its plane through bending action and may be subjected to small

CHAPTER III.

ELDS statement

21

III.2


Figure III.2: Generic rectangular plate representation

deections (displacements) if the out-of-plane displacement component w is much less than the thickness t, that is, w/t 1. Considering what stated, plates can be classied into three categories depending essentially on the thickness and deection component w as follows [44]: 1. thin plates with small deections; 2. thin plates with large deections; 3. thick plates. In the present work thin plates with small deections have been taken into consideration only.

III.2.2 Two-dimensional stresses and strains Two dimensional elasticity problems typically involve structures that are very thin such as plates. These types of problems are formulated in terms of inplane stresses and strains by assuming that stresses are negligible with respect to the smaller dimension (along z -direction) and are small compared to the in-plane stresses, and by removing from the stress analysis the out-of-plane component of the strain vector arising from the Poisson's eect.

CHAPTER III.

ELDS statement

22

III.2


Considering that, the stresses and strains formulations neglect all terms involving z in the strain-displacement relationship dening the three-dimensional case [33] and, accordingly, the stress and strain vectors are respectively given by

{σ} = {σx , σy , τxy }T {} = {x , y , γxy }T where σx and σy are, respectively, the normal component of stress parallel to x and y axes, τxy is the shearing in-plane stress component in coordinate directions (Cartesian), x and y are the unit elongation along x and y directions, γxy is the shearing in-plane strain component in coordinate directions, whilst the strain-displacement relationships are

x =

∂u ∂x

y =

∂v ∂y

γxy =

∂u ∂v + . ∂y ∂x

(III.2)

For an isotropic material, namely a material with a behaviour with no preferred directions, and in the case that the initial strain vector is null, the stress-strain relationship in matricial form is      σx   x  σy = [E] y     τxy xy where [E] is the material matrix that can be expressed according to the specic two-dimensional condition as dened afterwards (see equations III.3 and III.4). Figure III.3 illustrates the stresses present in an innitesimal thin plate in the case of a two-dimensional elasticity problem. There are two classes of problems in the plane elasticity analysis: 1. plane stress condition; 2. plane strain condition. Assuming that thin plate is in the xy plane, the plane stress condition foresees that stresses σz = τyz = τxz = 0, and z 6= 0 [19]. Then, for an isotropic

CHAPTER III.

ELDS statement

23

III.2


material in isothermal conditions, the material matrix [E] can be expressed as



1 E  ν [E] = 1 − ν2 0

ν 1 0

 0  0 0.5(1 − ν)

(III.3)

where E and ν are the modulus of elasticity (Young's modulus) in tension and compression and the Poisson's ratio of the material respectively. Though the component z exists, this is associated with zero stress and hence needs not to be included in an energy based analysis (see paragraph III.2.3). On the other hand, for the plane strain condition z = γyz = γzx = 0, and

σz 6= 0 [19]. The material matrix [E], for an isotropic material in isothermal conditions, is given by



1−ν E  ν [E] = (1 + ν)(1 − 2ν) 0

ν 1−ν 0

 0 . 0 0.5(1 − 2ν)

(III.4)

Similarly to the plane stress case, though the component σz exists, this is associated with zero stress and, therefore, needs not to be included in an energy based analysis. Considering the gure III.3, when in-plane forces are applied to the structure,

Figure III.3: Stresses within a dierential element of a plate

CHAPTER III.

ELDS statement

24

III.2


the displacements vector at any discrete point of the structure located by the coordinates x and y is

{d} = {u, v, 0}T where u and v are the components of the displacements vector along x and

y axes respectively.

III.2.3 Finite Element stiness matrix calculation Apart from a few very simple cases, the most direct procedure in solid mechanics to explicitly calculate the stiness matrix of an element [k] is the energy approach. According to this methodology, [k] can be determined considering the strain energy of a nite element Ue as Z Z 1 1 1 T T Ue = {} {σ}dV = {d} [B]T [D][B]dV {d} = {d}T [k]{d} 2 2 2 being

Z [k] =

[B]T [D][B]dV

(III.5)

where [B] is the strain interpolation matrix relating element nodal strains vector {} to nodal displacements vector {d}, namely

{} = [B]{d}, [D] is the modulus matrix relating element nodal stresses vector {σ} to nodal strains vector {}, namely

{σ} = [D]{}, and

R

dV indicates the integration over the volume of the nite element. In

the same conditions of study dened in paragraph III.2.2, the modulus matrix [D] has the same form of the material matrix [E]. Considering what above reported, to explicitly gain the mathematical expression of the components of the matrix [k], the integral on the right side of relationship III.5 needs to be evaluated. To this end, apart from a few cases for which analytical integration can be used, in most cases the numerical integration is needed.

CHAPTER III.

ELDS statement

25

III.2


III.2.4 Fundamentals of discretisation Successively the essential concepts on FEM discretisation are supplied. In view of making more compact the notations in some of the mathematical formulations reported hereafter, it will be convenient to utilize a comma to denote the partial dierentiation with respect to the subscript that follows.

Interpolation functions One of the crucial topics in the FEM approach is the choice of the suitable interpolation functions in order to appropriately approximating the variation of a eld variable over an element. As a matter of fact, FEM approximation is due to the lack of knowledge of the displacements eld within the element. To this end, this eld is approximated by means of interpolation functions that operate directly on nodal displacements. In general, to interpolate means to devise a continuous function that satises prescribed conditions at a nite number of points. Among those potentially employable, the polynomial functions are commonly used to interpolate a eld variable because of their dierentiation and integration turn out to be rather straightforward with respect to other types of functions. In the specic case of the FEM analysis, this type of interpolation functions is widely used by assuming that the conditioning points are the nodes of the element, and the prescribed conditions are the nodal values of the eld variables of interest. To guarantee the achievement of a monotonically convergent solution in computing, the approximate functions must full the following requirements [57]: 1. a continuous behaviour within the element; 2. the compatibility along the common nodes, boundaries or surfaces between adjacent elements; 3. the capability to represent a rigid body motion with uniform value of

CHAPTER III.

ELDS statement

26

III.2


eld variables and their spatial derivatives. A complete polynomial of order n in one dimension can be written as

Pn (x) =

n X

(III.6)

ai x i

i=0

whereas in two dimensions as (n+1)(n+2) 2

Pn (x, y) =

X

ak x i y j

i + j 6 n.

k=1

To deep the two-dimensional element analysis, constant, linear, and quadratic complete polynomials in two dimensions can be respectively written as

P0 (x, y) = a0 P1 (x, y) = a0 + a1 x + a2 y P2 (x, y) = a0 + a1 x + a2 y + a3 x2 + a4 xy + a5 y 2 . The Pascal's triangle shown in gure III.4 is useful for including the appropriate number of terms to obtain the complete approximating functions in any order. The coecient ai of the approximating polynomials are named

Figure III.4: Pascal's triangle for the complete polynomial characterization

generalized coordinates and are independent parameters not identied with specic nodes. They have not a particular physical meaning and specify the

CHAPTER III.

ELDS statement

27

III.2


magnitude of the prescribed distribution of the eld variable. If side and internal nodes in a triangular element have a regular pattern, it is feasible generating a family of straight-sided triangles so that each contains a complete polynomial in Cartesian coordinates x and y , as schematically shown in gure III.4.

Coordinate systems Apart from the afore dened generalized coordinates, the global and local coordinates as well as natural coordinates are frequently met in the FEM context. In particular, the global coordinates are suitable for identifying the position of each node, the orientation of each element, the loads, the boundary conditions, and the solution of the eld variables as well. The local coordinates, whose system origin is positioned within the element, are used to simplify the algebraic manipulation for the formulation of the element characteristic matrixes. The main objective of the use of natural coordinates is to determine the position of a point inside an element in terms of coordinates related to the nodes of the elements by means of dimensionless parameters whose absolute magnitude never exceeds unity. Even if the natural coordinates are function of global coordinates, they are commonly dened with respect to element system and enable to express approximating functions in such a manner to permit the utilization of special integration formulas to evaluate the integrals in the element matrixes formulation.

Natural coordinates in two dimensions Considering a two-dimensional triangular area dened by three nodes respectively named 1, 2, and 3, one at apex, a generic point P belonging to the triangle and denoted by (x, y) global coordinates, divides the triangle into three subareas A1 , A2 , and A3 as shown in III.5. Referring to a triangular

CHAPTER III.

ELDS statement

28

III.2


element, then, the point P (x, y) may be located in terms of parametric coordinates, commonly referred to as area coordinates, dened as ratios of areas as follows

A1 A2 A3 χ2 = χ3 = A A A where A is the total area of the triangle. Since A is the sum of the three subχ1 =

areas, these latter are not independent of each other, but they are supposed to verify the constraint relationship

χ1 + χ2 + χ3 = 1. As depicted in gure III.5, one specic area coordinate has a unit value at one node of the element and zero value at other nodes. Then, χi = constant represents a set of straight lines parallel to the side opposite to the ith corner. The equations of sides 12, 23 and 31 are, respectively, χ3 =0, χ1 =0 and

Figure III.5: Denition of area coordinates χ1 , χ2 , and χ3

χ2 =0. Therefore, the three corners have area coordinates (1,0,0), (0,1,0), and (0,0,1) respectively. In the case that midpoints are present (quadratic triangular elements), they have area coordinates (1/2,1/2,0), (0,1/2,1/2), and (1/2,0,1/2), the centroid (1/3,1/3,1/3), and so on.

CHAPTER III.

ELDS statement

29

III.2


Area coordinates can also be established in terms of ratios of lengths. Let

L23 be the length of side 23, η the distance from the internal point P (x, y) to side 23, and h the distance from the vertex 1 to side 23, then A1 = L23 η/2,

A = L23 h/2, and χ1 = A1 /A = η/h. In such a manner, the relationships between the area coordinates χ1 , χ2 , and χ3 and the area coordinates ξ , η can be done by regarding area coordinates as ratios of lengths as follows

χ1 = 1 − ξ − η

χ2 = ξ

χ3 = η

in the case that both side 12 and 13 are each unit long as illustrated in gure III.6 for a quadratic triangular element.

Figure III.6: Triangular element in the local system and area coordinates

In general, the equations characterizing the parameterization of the element require neither that sides of a triangle be straight nor side nodes be uniformly spaced. However, when these conditions are satised, the Jacobian is constant throughout the element, all terms to be integrated are polynomials in area coordinates, and integration formulas become quite simple [19].

Shape functions The shape functions constitute a subset of the element approximated functions and are such not to be arbitrarily chosen.

CHAPTER III.

ELDS statement

30

III.2


Taking into consideration the one-dimensional element with n nodes, from equation III.6 the approximated function for the generic eld variable φ(x) is

n X

φ(x) =

ai xi

i=0

or alternatively expressed in matricial form is (III.7)

φ(x) = {X}T {a} where

{X}T = {1 x x2 {a}T = {a0

a1

a2

···

xn }

···

an−1 }

with the number of generalized coordinates ai equal to the number of nodes within the element. The eld variable φ(x) can also be established within the element by means of its nodal values φi , and can thus be expressed as

φ(x) =

n X

Ni φi

i=1

or alternatively in matricial form as (III.8)

φ(x) = {N}T {φ} where

{N}T = {N1

N2

N3

{φ}T = {φ1

φ2

φ3

··· ···

Nn } φn }

in which Ni are usually referred to as shape functions or basis functions. The main properties of these latter functions are: 1. Ni = 1 at node i and Ni = 0 at all other nodes; 2.

Pn

i=0

CHAPTER III.

Ni = 1.

ELDS statement

31

III.2


The explicit expression of the shape functions Ni can be gained by solving for the generalized coordinates ai in terms of the nodal coordinates xi and nodal values φi using equation III.7, and then rearranging the resulting expression in the form of equation III.8. At each node j , the eld variable φ(x) is calculated as

φ(xj ) =

n X

ai xij

i=0

or in matrix form as

{φe } = [A]{a}. Solving for the generalized coordinates in terms of generalized coordinates and nodal values of the eld variables yields to

{a} = [A]−1 {φe }. Substituting for the generalized coordinates in equation III.7 results

φ(x) = {X}T [A]−1 {φe } namely

φ(x) = [N]{φe } where the matrix

[N] = {X}T [A]−1 is referred to as shape functions matrix. To summarize, the standard interpolation takes the following subsequent stages: 1. to choose interpolation functions; 2. to evaluate interpolation functions at known points; 3. to solve equations to determine unknown constants.

CHAPTER III.

ELDS statement

32

III.3

ELDS implementation

Dierentiating the shape functions [N], the strain-deformation matrix [B] can be obtained [19] and, then, the stiness matrix can be calculated considering the relationship III.5. Focusing the attention on triangular element, the shape functions of a linear triangle are the same as the area coordinates

N1 = χ 1

N2 = χ 2

N3 = χ 3

whereas for a quadratic triangle, by using the area coordinate ξ and η , are

N1 = 2(1 − ξ − η)(0.5 − ξ − η) N2 = ξ(2ξ − 1) N3 = η(2η − 1) N4 = 4ξη N5 = 4η(1 − ξ − η) N6 = 4ξ(1 − ξ − η). Specically, the shape functions above listed will be successively used in the description of the formulation of the CST (see paragraph III.3.3) and DKT element (see paragraph III.3.3) respectively.

III.3 ELDS implementation ELDS is primarily intended to provide the FLUENT user with a powerful numerical means to cope with the large displacements FSI problems of thin-walled structures. Two important considerations have mainly pushed its development. On the one hand the lack of a really ecient manner to approach this specic task through the CFD commercial tool under discussion and, moreover, the fact that, according to a rough estimate, about 70% of the FEM industrial mechanics problems are concerned with shell structures [9], namely curved components with thickness much smaller compared to their curvilinear dimensions.

CHAPTER III.

ELDS statement

33

III.3

ELDS implementation

For its rst release, ELDS has been thought to run fully embedded into the FLUENT frame and supply the feature to completely set up the FEM solution in the same working GUI (Graphical User Interface). This strategical method to structure ELDS oers the undeniable advantage of the ease of use, avoiding thus to stress the user in training with a third party tool with the consequent saving of time and money. It is worth mentioning that the development of ELDS was scheduled in the framework of the FSI Project started at the Mechanical Department of the University of Rome "Tor Vergata". Since its inception, the FSI Project Team has tackled the challenging topic of the uid-structure interaction producing several publications during years [14] [15] [16] [17]. Specically referring to algorithm implementation, in a sense ELDS can be regarded as a follow-up of a previous consultancy activity carried out to satisfy an industrial requirement concerning the numerical reproduction of the behaviour of paper sheets subjected to air ow in the printer industry. In the following subsections the implementation of the ELDS solver according to the FEM co-rotational technique is described. To comprehend the numerical solution environment established to handle geometrically non-linear analyses of shell structures, the following topics are required to be outlined:

• the FEM co-rotational theoretical approach; • the characterization of the core of the explicit numerical solution; • the formulation of the stiness matrix of the general element used by the FEM solver;

• the characterization of the technique for deformable components nodes management. In the following sections, all above-listed topics are respectively detailed.

CHAPTER III.

ELDS statement

34

III.3

ELDS implementation

III.3.1 Co-rotational formulation The geometrically non-linear analysis of shell structure concerns the nonlinear equilibrium path for the material behaviour. It is commonly recognized that this type of numerical analysis is accomplished iteratively for each load increment with a subsequent updating of the coordinates and internal stresses, because equilibrium is achieved in the deformed conguration and not in the starting conguration as in the case of linear analysis [27]. As such, the numerical simulation of the system behaviour is carried out by solving a set of subsequent linear analyses that approximate the non-linear response of the system itself. The numerical techniques used to describe the structural behaviour of a system, in terms of the displacements response can be divided into two methods which are referred to as TL (Total Lagrangian) and UL (Updated Lagrangian) formulation [5]. These methods dier in the conguration they take as reference, usually termed as reference conguration, to evaluate the response of the structure at each iteration, that is the initial and previous conguration respectively [6]. The pioneering works of Wempner [54], and Belytschko and Hsieh [10] introduced around 70s the co-rotational approach in the nite element analysis. In particular, Wempner formulated the co-rotational procedure according to which the response of the structure at each increment (iteration), is dened utilizing a local co-rotational frame which continuously translates and rotates (co-rotate) with the nite element. Moreover, he proposed to decompose the motion of each single element into a rigid body motion and pure deformation. In such a manner, the geometrical non-linearities are totally assigned to the rigid body rotations whilst the pure deformations are attributed to the constitutive equations of the element. Hypothesizing the pure deformation contribution to be small, in the local system the linear theory can be employed. This approach enables to simplify the TL and the UL analysis and,

CHAPTER III.

ELDS statement

35

III.3

ELDS implementation

due to its generality in the implementation, as well as its relative simplicity in handling large displacements and rotations of the components, it enjoyed the immediate success in the setting up of nite elements for beams and shell element types pushed by the attractive characteristic to use linear elements in a non-linear context. An extensive literature review on at triangular shell nite element for geometrically non-linear analysis according to co-rotational technique can be found in Gal and Levy [27]. Employing the co-rotational approach for FEM solution, in order to discriminate for each node at each calculation time step the part of displacement due to the rigid body motion and pure deformation, ELDS uses the convention shown in gure III.7.

Figure III.7: Transformation from the reference to actual time element conguration

Specically, for a generic triangular element belonging to a deformable component, the following assumptions are made:

• the denition of the element nodes at apexes 1, 2, and 3 adopts the counterclockwise direction;

CHAPTER III.

ELDS statement

36

III.3

ELDS implementation

• the element lies in the xy plane dened by the local system; • the local system has the origin that coincides with the in rst point of the element and the x axis (called xloc in gure III.7) and its corresponding unit vector x ˆ lays on the segment connecting node 1 and 2. The z axis (called zloc in gure III.7) is determined through the vector product between x ˆ and the vector connecting node 1 and 3, whilst the

y axis (called yloc in gure III.7) is identied through its unit vector yˆ calculated by means of the vector product between zˆ and x ˆ;

• the reference conguration C ref is the starting conguration, namely that at time of calculation tref =0. Whatever the calculation time, indicating with C tn+1 the actual conguration at time tn+1 , the pure deformation is obtained by carrying out four subsequent operations. These operations respectively foresee the calculation of: 1. the element nodes location in the actual conguration with respect to t

the global system, that is the vectors ri n+1 for i=1, 2, and 3; t

2. the total displacement for each node, that is the vector ri ref − ri n+1 for

i=1, 2, and 3; 3. the rigid rotation that enables the local system at reference conguration time to coincide with the local system at actual time (the gained conguration is that highlighted with the dashed segments in gure III.7); 4. the pure deformation contribution for each node by subtracting the rigid rotation to total displacement. Since the node 1 is taken as reference for rotation, only nodes 2 and 3 give contribution to pure elastic deformation. The displacement vectors due to

CHAPTER III.

ELDS statement

37

III.3

ELDS implementation

pure deformation referring to node 2 and 3 are indicated in gure III.7, el respectively, as drel 2 and dr3 .

Since according to the explicit method the displacements vector at previous time {d}tn is known, once the components of the stiness matrix of the element [k]ELDS have been calculated (see section III.3.3), employing such a methodology the elastic nodal force vector {Fel }tn can be evaluated by the following relationship

{Fel }tn = [k]ELDS {d}tn . As soon as this vector is determined, the nodal displacements vector {d}tn+1 calculation process can be enabled as described in section III.3.2.

III.3.2 ELDS solution characterization The core of the explicit solution foresees the integration at current time tn+1 of the dierential equation in terms of the second derivatives of nodal dis¨ placements vector {d} to nally gain the vector {d}t , that is the nodal tn+1

n+1

displacements vector at current time tn+1 . Once this latter vector is available, the updating of the shape conguration of the deformable components can be accomplished by adding {d}tn+1 to the nodal positions vector known at previous conguration at time tn .

˙ In order to proceed a step in time, the nodal velocities vector {d} tn+1 and the nodal displacements vector {d}tn+1 in the ELDS algorithm are estimated by means of a central dierence scheme. This method uses the general dynamic equilibrium equation of time tn to predict the solution at current time tn+1 . The equilibrium determines the accelerations vector at the beginning of the increment which is assumed to be constant over the time step. The nodal displacements vector at current time tn+1 can be thus obtained by calculating

CHAPTER III.

ELDS statement

38

III.3

ELDS implementation

in sequence the following equations −1 ¨ {d} tn+1 = [M] {Fnodal }tn

(III.9)

˙ ˙ ¨ {d} tn+1 = {d}tn + ∆t{d}tn+1 ˙ {d}tn+1 = {d}tn + ∆t{d} tn+1 being {Fnodal }tn the force nodal vector. As introduced in section III.1.1, to stabilize and produce accurate numerical results, the explicit codes need that the time step ∆t must be chosen smaller than the smallest natural period in the mesh, namely the length of time it takes a signal traveling at the speed of sound in the material, to traverse the distance between the node points. Therefore, to prevent the achievement of computational errors and numerical instability, the time step used shall be less than the critical time step identied by means of the Courant condition as follows

∆tcrit =

lmin am

(III.10)

where lmin is the smallest distance between two dierent nodes of the surface mesh of the deformable components, and am is the speed of sound in the material. This relationship, also called as CFL (Courant-Friedrichs-Lewy) condition, is a necessary condition for ensuring the convergence of the numerical solution of certain partial dierential equations [20] [21]. Specically referring to shell structures, the equation III.10, used for the ELDS calculation, assumes the following form r ρ(1 − ν 2 ) ∆tELDS = S lmin E

(III.11)

where, to be safer, the safety factor parameter S , that is a coecient minus than unity, has been introduced to further reduce the time step size. Adopting the afore reported relation, the value for the calculation time step fulls the constant acceleration assumption over the timestep as well.

CHAPTER III.

ELDS statement

39

III.3

ELDS implementation

Taking into account equations III.9, once ∆tELDS is determined through equation III.11, since the mass matrix [M] can be calculated as afterwards de˙ scribed in paragraph III.3.4 and both {d} and {d}t are known at each time tn

n

step tn , the only unknown term to carry out the explicit solution at time tn+1 is {Fnodal }tn . This vector can be calculated by exploiting the information available at time tn from the following relation

{Fnodal }tn = {Fext }tn − {Fel }tn being {Fext }tn the external force vector at that time. In turn, these quantities can be respectively evaluated from the following equations

{Fext }tn = {Ffluid }tn + {Fgrav }tn {Fel }tn = [k]ELDS {d}tn where {Ffluid }tn and {Fgrav }tn are the nodal force vector due to the uid loads and gravity respectively. In this version of the algorithm, other types of force such as those due to contact and temperature are not considered.

III.3.3 ELDS element stiness matrix formulation Since all shell structures can be approximated by piecewise triangular at elements, a general at three-noded triangular element has been developed. This element is obtained by properly combining the stiness matrix of three dierent elements, each of which contributes with its characteristic DOF. These elements are listed below and are termed as follows: 1. CST; 2. DKT; 3. DRILL.

CHAPTER III.

ELDS statement

40

III.3

ELDS implementation

In particular, the CST (Constant Strain Triangle) membrane element contributes to two in-plane translations, the plate bending DKT (Discrete Kirchho Theory) element contributes to the out-of-plane translation and two in-plane rotations, whilst the DRILL triangle adds one ctitious rotational degree of freedom to each node for the out-of-plane rotation. Figure III.8 shows the degrees of freedom of ELDS complete triangular element, where the convention adopted for nodes numbering is the counterclockwise direction, and the single arrow and the double arrow represent, respectively, the translation and the rotation in the direction of the local reference system axes.

Figure III.8: Building-up of ELDS general at triangular element

Afterwards, the main theoretical and mathematical passages for the implementation of the CST, DKT, and DRILL element stiness matrix are respectively described as well as their assembling for the nal achievement of the ELDS one.

Constant strains triangle (CST) In this subsection, the development of the CST element stiness matrix following the demonstration reported in [19] is presented.

CHAPTER III.

ELDS statement

41

III.3

ELDS implementation

This membrane element is one of the earliest introduced [45] and, despite that it represents the simplest triangular plane stress element and has limits such as constant stresses and in-plane share locking [32], it is widely used because it can approximate the behaviour of small regions of in-plane plate deforming and, consequently, it covers an important class in engineering design problems. A membrane is dened as a thin plate that has neither rotational stiness nor stiness normal to the plane of the element. It can be situated arbitrarily in space but the resultant forces must lie in the plane of the element. Figure III.9 depicts the CST element where the coordinates of nodes 1, 2, and 3 are (x1 , y1 ), (x2 , y2 ), and (x3 , y3 ) respectively, whilst the corresponding displacements are (u1 , v1 ), (u2 , v2 ), and (u3 , v3 ). As previously mentioned, each node has two in-plane transactional DOF for a total of six DOF. In general,

Figure III.9: Constant Strain Triangle (CST)

this type of element is properly used in the areas where strain gradients are small and the mesh transition is present, whereas it shall not be utilized in critical areas such as stress concentration zones, edges of holes, and corners. Moreover, in the case of a general analysis purpose, a very large number of these elements is required for reasonable accuracy in computational results. As far as concerns the stiness matrix mathematical formulation, the CST is

CHAPTER III.

ELDS statement

42

III.3

ELDS implementation

founded on the hypothesis that the displacements eld over the element can be expressed through a complete linear polynomial in both x and y directions respectively by

u(x, y) = a0 + a1 x + a2 y (III.12)

v(x, y) = a3 + a4 x + a5 y

where ai for i=0 to 5 are constant values. Alternatively, relationships III.12 can be written in matrix form as

u(x, y) 1 = v(x, y) 0

x 0

y 0

0 1

0 x

  a0         a   1     0 a2 y  a3        a   4     a5

(III.13)

or again in compact representation as

{U(x, y)} = [X]{a}.

(III.14)

Considering relationship III.2, by properly deriving and manipulating equations III.12, the following relationships for strains can be obtained

x = a1 y = a5 γxy = a2 + a4 . According to the above relationships, the expression of the strains vector involves only constant terms. For this reason, this element has been called constant strain triangle. Substituting the values of the nodal coordinates and the corresponding displacements into equation III.13, the following system can be obtained      u 1 x y 0 0 0 a0     1 1 1              v1  0 0 0 1 x1 y1   a1               u2 1 x y 0 0 0 a 2 2 2   = v2  0 0 1 x2 y2  a3     0              u 1 x y 0 0 0 a     3 3 3 4         v3 0 0 0 1 x3 y3 a5

CHAPTER III.

ELDS statement

43

III.3

ELDS implementation

which can be more compactly represented in matricial form as (III.15)

{u} = [A]{a}.

Substituting into equation III.14 the coecients vector {a} obtained by inverting equation III.15, the relationship linking element displacements vector with nodal displacements one can be gained, and results to be the following

{U(x, y)} = [X][A]−1 {u} where the inverse of the matrix containing nodal coordinates [A]−1 is

[A]−1

 B11 B21  1  B31 = 2A   0  0 0

0 0 0 B42 B52 B62

B13 B23 B33 0 0 0

0 0 0 B44 B54 B64

B15 B25 B35 0 0 0

 0 0   0   B46   B56  B66

(III.16)

where

B11 = x2 y3 − x3 y2

B13 = x3 y1 − x1 y3

B15 = x1 y2 − x2 y1

B21 = y2 − y3

B23 = y3 − y1

B25 = y1 − y2

B31 = x3 − x2

B33 = x1 − x3

B35 = x2 − x1

B42 = x2 y3 − x3 y2

B44 = x3 y1 − x1 y3

B46 = x1 y2 − x2 y1

B52 = y2 − y3

B54 = y3 − y1

B56 = y1 − y2

B62 = x3 − x2

B64 = x1 − x3

B66 = x2 − x1

and A is the area of the triangle that can be expressed as

1 A = [x1 (y2 − y3 ) + x2 (y3 − y1 ) + x3 (y1 − y2 )] 2 whilst, as previously described in paragraph III.2.4, the matrix product

[X][A]−1 represents the shape functions [N], that is [N] = [X][A]−1 .

CHAPTER III.

ELDS statement

44

III.3

ELDS implementation

The shape function matrix [N] can be thus explicitly formulated by multiplying the terms of [X] of relation III.13 for those of [A]−1 matrix of relation III.16. The latter passage yields to the following relationship     N1 (x2 y3 − x3 y2 ) + (y2 − y3 )x + (x3 − x2 )y N2  = 1 (x3 y1 − x1 y3 ) + (y3 − y2 )x + (x1 − x3 )y  2A (x1 y2 − x2 y1 ) + (y1 − y2 )x + (x2 − x1 )y N3

(III.17)

and, thus, the element displacements can be written in terms of the shape functions as follows

  u1         v   1     u(x, y) N1 0 N 2 0 N 3 0 u2 = . v(x, y) 0 N 1 0 N2 0 N 3  v2       u     3  v3

(III.18)

Considering the system III.18 and the strain-displacement relationships III.2, the components of the strains vector are respectively given by

3 ∂ X x = Ni (x, y)ui ∂x i=1 3 ∂ X y = Ni (x, y)vi ∂y i=1 γxy

3 3 ∂ X ∂ X = Ni (x, y)ui + Ni (x, y)vi ∂y i=1 ∂x i=1

that, in turn, can be written in matrix form as

   N1,x  x   y = 0   γxy N1,y

0 N1,y N1,x

N2,x 0 N2,y

0 N2,y N2,x

N3,x 0 N3,y

   u1      v1     0    u 2  N3,y  v2   N3,x     u3        v3

or in more compact manner as

{}3

CHAPTER III.

ELDS statement

x

1

= [B]3 6 {u}6 1 . x

x

45

III.3

ELDS implementation

For what just reported, by taking the derivatives of the shape functions III.17 with respect to x and y , the terms of the strain-displacement matrix [B] can be explicitly gained   y2 − y3 0 y3 − y2 0 y1 − y2 0 1  B = 0 x3 − x2 0 x1 − x3 0 x2 − x1  . 2A x3 − x2 y2 − y3 x1 − x 3 y3 − y2 x2 − x1 y1 − y2 After performing the afore reported passages, the element stiness matrix can be obtained using the strain-displacement matrix [B] and material matrix [E] as follows

Z [k]6

6

x

[B]T6 3 [E]3 3 [B]3 6 dV.

=

x

x

x

V

For constant thickness element, the volume integral can be reduced to the area integral

Z [k]6

6

x

(III.19)

t[B]T6 3 [E]3 3 [B]3 6 dA

=

x

x

x

A

where t is the thickness of the element. As all terms in the strain-displacement matrix [B] and material matrix [E] are constant, the equation III.19 can be rewritten as

Z [k]6

6

x

=t

[B]T6 3 [E]3 3 [B]3 6 dA x

x

x

A

or as

[k]CST = [k]6

6

x

= tA[B]T6 3 [E]3 3 [B]3 6 . x

x

x

(III.20)

Since the membrane element CST contains constants only, and because of the lack of x and y terms associated to strain eld, the CST element weakly reacts in bending.

Discrete Kirchho Theory (DKT) triangle With the purpose of explaining the meaning of some denitions used afterwards to describe the DKT stiness matrix formulation, a brief introduction of the Kirchho theory is rstly provided. As mentioned before, in the present work thin plates with small deections

CHAPTER III.

ELDS statement

46

III.3

ELDS implementation

have been taken into consideration. In particular, in the theory of bending for thin plates there are three basic assumptions [44]: 1. the mid-surface of the plate remains unstretched during deformations; 2. the points straight and normal to the mid-surface of the plate before bending remain straight and normal to the mid-surface after bending; 3. the transverse shear stresses are small compared to normal stresses and, hence, can be neglected. These hypotheses are known as Kirchho's assumptions and are applicable to the bending of thin plates with small deections. Considering an isotropic plate of uniform thickness t with the xy plane as the principal plane and a small section of the plate of length dx in the x direction as shown in gure III.10, when a load is applied in the z direction, the point O on the mid-surface of the plate moves in z direction as the plate deforms due to bending. According to the Kirchho's assumptions, a

Figure III.10: Kirchho's assumptions

line that is straight and normal to the mid-surface before bending remains straight and normal to the mid-surface after bending. In such conditions, the displacement components u, v , and w at any point of coordinates x, y ,

CHAPTER III.

ELDS statement

47

III.3

ELDS implementation

and z can be represented as [19]

∂w = −zθx (x, y) ∂x ∂w v = −z = −zθy (x, y) ∂y w = w(x, y) u = −z

(III.21)

where, w is the transverse displacement, and θx = ∂w/∂x and θy = ∂w/∂y are the rotations of the normal to the undeformed surface in the xz and yz planes respectively. The displacement-curvature relationships for the thin plate can be respectively written as follows

∂2w = w,xx ∂x2 ∂2w κy = 2 = w,yy ∂y ∂2w κxy = = w,xy . ∂x∂y κx =

(III.22)

From equations III.21, and the strain-displacement relationship III.2, the components of the strains vector can be written as

∂2w = −zw,xx ∂x2 ∂2w y = −z 2 = −zw,yy ∂y ∂2w γxy = −2z = −2zw,xy . ∂x∂y x = −z

(III.23)

According to the theory of bending for a thin plate, the plate is in the plane stress condition and hence all stresses vary linearly over the thickness of the plate. The moments per length can thus be represented as − 2t

Z

σx zdz

Mx =

ELDS statement

t 2

− 2t

Z σy zdz

My =

t 2

CHAPTER III.

− 2t

Z

τxy zdz

Mxy = t 2

(III.24)

48

III.3

ELDS implementation

where Mx and My are the moments in the x and y direction respectively, and

Mxy is the twisting moment. By substituting the strains of relations III.23, the plane-stress expressions III.3 can be written as    1  σx  E  ν σy =   1 − ν2 0 τxy

ν 1 0

  0  −zw,xx   −zw,yy 0   0.5(1 − ν) −2zw,xy

and then the equations of stresses in the plane are Ez ∂ 2 w ∂2w σx = − +ν 2 1 − ν 2 ∂x2 ∂y 2 Ez ∂ w ∂2w σy = − +ν 2 1 − ν 2 ∂y 2 ∂x 2 Ez ∂ w . τxy = − 1 + ν ∂x∂y

(III.25)

(III.26)

By substituting the stresses from equation III.26 to III.24 and by integrating over the thickness of the plate, the following relationships for the moments per length are obtained

2 Et3 ∂ w ∂2w Mx = − +ν 2 12(1 − ν 2 ) ∂x2 ∂y 2 3 ∂ w ∂2w Et +ν 2 My = − 12(1 − ν 2 ) ∂y 2 ∂x 2 3 Et ∂ w Mxy = − . 12(1 + ν) ∂x∂y

(III.27)

The equation III.27 can be represented in matrix form as

   1  Mx   ν My = −D   Mxy 0

ν 1 0

  0  w,xx   w,yy 0   0.5(1 − ν) 2w,xy

(III.28)

where

Et3 12(1 − ν 2 ) is dened the exural rigidity of the plate. D=

From equations III.28 and III.22, the moment-curvature relationships is given

CHAPTER III.

ELDS statement

49

III.3

ELDS implementation

by

   −D  Mx   My = −Dν   Mxy 0

−Dν −D 0

  0  κx   κy 0 .   −0.5D(1 − ν) 2κxy

From computational point of view, since the earliest development of FEM, a valuable amount of research has been devoted to the formulation of plate and shell elements. Among the large number of bending elements that have been built up, the DKT triangle developed by Batoz [7] has been chosen. This at triangular plate bending element, previously shown in gure III.8, has been suitably designed to be eciently implemented to solve thin plate analysis by means of numerical calculation on computers. The DKT element has displacements at three corner nodes only, and it has been found during time to be one of the most reliable, cost-eective, and accurate of its class [19]. As such, this element is extensively used as triangular plate bending element in FEA (Finite Element Analysis) codes. Its development is based on the Kirchho theory assumptions, according to which the bending energy present in the element is much higher compared to the shear strain energy and, thus, the transverse shear energy term can be neglected in the energy equation formulation. Considering that, the bending energy can be represented in the following form [7] Z 1 Ub = {κ}T [D]b {κ}dxdy 2 A where



1 Et  ν [D]b = − 12(1 − ν 2 ) 0 3

ν 1 0

(III.29)

 0 , 0 0.5(1 − ν)

t is the thickness of the plate, and the curvature vector is given by   βx,x κ =  βy,y  , βx,y + βy,x

CHAPTER III.

ELDS statement

50

III.3

ELDS implementation

that constitutes the target for the DKT element stiness matrix formulation. In equation III.29 the only variables in form of rst derivative are βx (x, y) and βy (x, y), whilst the displacement component w does not appear. Taking into account for this latter consideration, in order to correlate the rotations of the mid-surface normal to the transverse displacement w, Batoz [7] made the following assumptions:

• the triangular element should have only nine DOF, namely the transverse displacement w and the rotations βx and βy at each vertex node of the element;

• according to Kirchho's theory, the rotations can be dened as βx = w,y

βy = −w,x

(III.30)

• the Kirchho's theory can be imposed at any discrete point in the element;

• the compatibility of rotations βx and βy can not be lost. The starting point of the implementation of the 9 DOF DKT element is a straight-sided element with corner nodes previously depicted in gure III.6. Two fundamental assumptions were made by Batoz [7], namely: 1. the relationship between the rotations at six nodal points including mid-surface nodes and the shape functions at each of six nodes is in the form of a quadratic

βx =

6 X

Ni (ξ, η)βxi

βy =

i=1

6 X

Ni (ξ, η)βyi

(III.31)

i=1

where, for i =1 to 6, βxi and θyi are the rotations at each node as shown in gure III.13, and Ni (ξ, η) are the shape functions on the master parametric element (see gure III.6) in area coordinates ξ and

η already dened in paragraph III.2.4;

CHAPTER III.

ELDS statement

51

III.3

ELDS implementation

2. lateral deection w along each edge is assumed to be a cubic in an edge-tangent coordinate s. The cubic expression is as follows

w,sk = −

1 1 3 3 wi − w,si + wj + w,sj 2lij 4 2lij 4

where, k is the mid-node of the ij side of the triangle, whilst lij represents the length of the ij side of the triangle. To provide the mathematical information necessary to obtain the stiness matrix elements, the following conditions are imposed: 1. the Kirchho's hypothesis is applied to remove transverse shear strain

γ . According to this assumption, the following equations at the corner nodes 1, 2, and 3 are valid, that is βx + w,x γ= = 0. βy + w,y At the mid nodes 4, 5, and 6, βsk +w,sk = 0 where k is the node number; 2. the variation of the rotations along the sides of the triangle is represented by a linear equation

1 βnk = (βni + βnj ) 2 where, k =4, 5, and 6 represent the mid-nodes of the sides 23, 31, and 21 respectively. The conditions just expressed provide respectively 6, 3, 3 constraints for a total of 12 constraints. The latter four assumptions imply the condition that the transverse shear strain along the sides of the triangle βs + w,s = 0 is satised. It can be demonstrated [7] that βx and βy are given by the following relationships

βx = {Hx (ξ, η)}T {U }

CHAPTER III.

ELDS statement

βy = {Hy (ξ, η)}T {U }

(III.32) 52

III.3

ELDS implementation

where {U} is the displacements vector that can be written at each node as

{U}T = {w1

βx1

βy1

w2

βx2

βy2

w3

βy3 }

βx3

and Hx and Hy are the vectors components of the shape functions. Batoz represented the vectors components of the shape functions in the following form



 1.5(a6 N6 − a5 N5 )  (b5 N5 + b6 N6 )    N1 − c5 N5 − c6 N6     1.5(a4 N4 − a6 N6 )     (b N + b N ) Hx (ξ, η) =  6 6 4 4   N2 − c6 N6 − c4 N4     1.5(a5 N5 − a4 N4 )     (b4 N4 + b5 N5 )  N3 − c 4 N4 − c 5 N5



 1.5(d6 N6 − d5 N5 ) −N1 + e5 N5 + e6 N6     −b5 N5 − b6 N6     1.5(d4 N4 − d6 N6 )     Hy (ξ, η) =  −N + e N + e N 2 6 6 4 4    −b6 N6 − b4 N4     1.5(d5 N5 − d4 N4 )    −N3 + e4 N4 + e5 N5  −b4 N4 − b5 N5

where

ak = −

xij 2 lij

bk =

2 ck = (0.25x2ij + 0.5yij2 )/lij

3 xij yij 2 4 lij

2 dk = −yij /lij

2 ek = (0.25yij2 − 0.5x2ij /lij

xij = xi − xj 2 lij = (x2ij − yij2 )

yij = yi − yj

and k =4, 5, and 6 for the sides ij =23, 31, and 12 respectively. The strain-displacement matrix for the DKT element can be represented in the following form

  1   {Bξ,η } = 2A  

T T y31 Hx,ξ + y12 Hx,η T −x31 Hy,ξ

−

T x12 Hy,η

     

(III.33)

T T T T −x31 Hx,ξ − x12 Hx,η + y31 Hy,ξ + y12 Hy,η

CHAPTER III.

ELDS statement

53

III.3

ELDS implementation

where 2A = x31 y12 − x12 y31 . The derivatives of the vectors components of the shape functions with respect to ξ and η can be described by the following equations III.34, III.35, III.36, and III.37. In particular, the derivatives of the vectors components with respect to ξ are



 P6 (1 − 2ξ) + (P5 − P6 )η   q6 (1 − 2ξ) − (q5 + q6 )η   −4 + 6(ξ + η) + r6 (1 − 2ξ) − η(r5 + r6 )     −P6 (1 − 2ξ) + η(P4 + P6 )     Hx,ξ =  q6 (1 − 2ξ) − η(q6 − q4 )   −2 + 6ξ + r6 (1 − 2ξ) + η(r4 − r6 )      −η(P + P ) 5 4     η(q4 − q5 ) −η(r5 − r4 )   t6 (1 − 2ξ) + η(t5 − t6 )  1 + r6 (1 − 2ξ) − η(r5 + r6 )     −q6 (1 − 2ξ) + η(q5 + q6 )     −t6 (1 − 2ξ) + η(t4 + t6 )     1 + r (1 − 2ξ) + η(r − r6) Hy,ξ =  6 4    −q6 (1 − 2ξ) − η(q4 − q6 )      −η(t + t ) 4 5     η(r4 − r5 ) −η(q4 − q5 )

(III.34)

(III.35)

whilst the derivatives of the vectors components with respect to η are   −P5 (1 − 2η) − ξ(P6 − P5 )   q5 (1 − 2η) − ξ(q5 + q6 )   −4 + 6(ξ + η) + r5 (1 − 2η) − ξ(r5 + r6 )     ξ(P4 + P6 )    Hx,η =  ξ(q − q ) (III.36) 4 6     −ξ(r − r ) 6 4     P (1 − 2η) − ξ(P + P ) 5 4 5     q5 (1 − 2η) + ξ(q4 − q5 ) −2 + 6η + r5 (1 − 2η) + ξ(r4 − r5 )

CHAPTER III.

ELDS statement

54

III.3

ELDS implementation



 −t5 (1 − 2η) − ξ(t6 − t5 ) 1 + r5 (1 − 2η) − ξ(r5 + r6 )    −q5 (1 − 2η) + ξ(q5 + q6 )      ξ(t + t ) 4 6     ξ(r4 − r6 ) Hy,η =     −ξ(q − q ) 4 6    t5 (1 − 2η) − ξ(t4 + t5 )    1 − r5 (1 − 2η) + ξ(r4 − r5 ) −q5 (1 − 2η) − ξ(q4 − q5 )

(III.37)

where

Pk = −6 rk = 3

xij = 6ak 2 lij

qk = 3

yij2 2 lij

xij yij = 4bk 2 lij

tk = −6

yij = 6dk 2 lij

and k =4, 5, and 6 for the sides ij =23, 31, and 12 respectively. The strain-displacement matrix can be determined through relation III.33, and considering the equations III.34, III.35, III.36, and III.37, and by substituting the strain-displacement matrix, the element stiness matrix DKT can be obtained by the following relationship Z 1 Z 1−η [k]DKT = 2A {B}T [D]b {B}dξdη. 0

0

For an element with uniform thickness, [k]DKT is exactly integrated by a three-point quadrature rule. To implement the stiness matrix of the DKT element in ELDS, the FORTRAN coding provided by Jeyachandrabose [31] has been used. Once the global coordinates of the nodes of the generic triangle are passed, this algorithm automatically and explicitly generates the DKT stiness matrix. This goal is achieved by means of a procedure based on a global system formulation [31] which uses a number of operations less than the explicit approach formulated in local coordinates proposed by Batoz [8].

CHAPTER III.

ELDS statement

55

III.3

ELDS implementation

DRILL This matrix provides at each node of the triangular element the out-of-plane ctitious rotation θz with no resistance to any other rigid-body motion. The DRILL matrix is a 3x3 matrix, which multiplied to the element-normal nodal rotations produces the corresponding moments Mz , is dened as follows   1.0 −0.5 −0.5 1.0 −0.5  [k]DRILL = [k]3x3 = αEV  −0.5 (III.38) −0.5 −0.5 1.0 where V is the element volume and α is a coecient lower than 0.3 [19]. This matrix is combined to others previously discussed in order to prevent the potential ELDS stiness matrix singularity that may occur during calculation, when elements belonging to a generic deformable component assume the coplanar conguration.

ELDS stiness matrix assembling The stiness matrix of the ELDS element [k]ELDS can be nally obtained by properly assembling the stiness CST, DKT, and DRILL matrixes. This operation is schematically represented as follows

[k]ELDS = [k]CST + [k]DKT + [k]DRILL Once the nodal displacement vector {d} is known, the nodal elastic force vector at time tn can be determined by multiplying [k]ELDS with {d}, namely

{Fel }tn = [k]ELDS {d}tn

  ui         v   i     wi . = [k]ELDS θxi        θyi        θzi tn

(III.39)

Since a curved thin-walled component, and even more at one, can always be approximated as a faceted surface by connectiong at triangular elements together at vertex nodes, the ELDS approach has a general use.

CHAPTER III.

ELDS statement

56

III.3

ELDS implementation

III.3.4 Explicit method nodes management As described in section III.3.2, according to the explicit scheme approach the solution in terms of actual displacements is achieved by integrating the Newton's equation at each single node. In order to accomplish this task, a numerical support infrastructure to operate nodal forces calculation as well as the knowledge on nodes orientation and rotation to build up the element stiness matrix for each element are required. To this end, the triangular elements conguration of the generic single deformable thin-walled component is suitably manipulated so as to obtain an equivalent nodal masses conguration. This operation is illustrated in gure III.11 in case of an oriented structured grid of a rectangular-shaped shell component. In such a manner, each node can be thought treated as a sin-

Figure III.11: Mesh manipulation to nodal masses conguration achievement

gle sphere with dierent volume and orientation in space depending on the location in the native corresponding triangular surface mesh conguration. Moreover, each node interacts with other nodes sharing the native element through the element stiness matrix assembled as earlier described in section III.3.3. Considering what just stated, the core of ELDS deals with the evaluation of: 1. nodal masses; 2. nodal directions and rotations.

CHAPTER III.

ELDS statement

57

III.3

ELDS implementation

Concerning the determination of nodal masses, the mass of each node is determined by adding the contribution of one-third of the mass of the adjacent elements of the native triangular mesh. The nal result is qualitatively shown in gure III.11. This operation is carried out straightforwardly by means of one specic MACRO function provided by the UDF (User Dened Function) feature of FLUENT (see chapter IV) performed once at FEM initialization (see paragraph III.4). As a matter of fact, because of the fundamental small strains assumption, both the shape and area of the elements do not significantly change and, hence, the nodal masses neither. So, the updating of nodal masses during the FSI simulation does not need to be performed. Referring to the calculation of nodal directions, this task is accomplished in two subsequent phases. In the rst phase the element normal vector embedding the area value is acquired through another MACRO function. Successively, the nodal vector of each node is calculated by adding the one-third of the contribution of the adjacent elements of the native triangular mesh. Figure III.12 depicts the averaged operation for the achievement of the node normal of a generic node. As far as nodal rotations are concerned, these can be obtained once nodal directions are available by using nodal directions and normal vector of the element of the native triangular mesh according to the adopted convention shown in gure III.13. Since during the FSI transient calculation, the triangles dicretizing the deformable components may change their normal vector, the direction of each node belonging to the element changes accordingly. As such, the direction for each node needs to be computed at each FEM time step to gain the actual values of the parameters of interest. In particular, taking into account for co-rotational approach on which ELDS is founded (see section III.3.1), to determine the relative rotation, the native triangular conguration at starting time is taken as reference. It is worth mentioning that an important feature of the ELDS is the capability to manage whatever starting conguration shape of the deformable

CHAPTER III.

ELDS statement

58

III.3

ELDS implementation

Figure III.12: A generic node normal calculation

thin-walled component of the system, unless internal sharp edges are present. To eectively carry out this operation, all data necessary to dene the starting conguration are stored so that the dierential calculation with respect

Figure III.13: Nodal rotation convention

CHAPTER III.

ELDS statement

59

III.4

Overall ELDS operations

to the actual time tn+1 during the dynamic transient simulation can be done.

III.4 Overall ELDS operations The ELDS overall process can be seen as composed by two main macrotasks, which are referred to as ELDS initialization and computing. In turn, these macrotasks are subdivided into several operations, that are performed in sequence as depicted through a ow-chart scheme in gure III.14 and III.15 respectively. In particular, the ELDS initialization task rstly foresees the manual imposition of all parameters required for the set-up of the FEM solution. Once this action has been accomplished, it can be possible to automatically perform by means of an EXECUTE_ON_DEMAND MACRO (see paragraph IV.3.2) the remaining operations of the FEM initialization (to ease the comprehension all these actions are grouped together and indicated as FEM initialization in gure IV.11 of chapter IV). This latter phase specically concerns:

• the acquisition of the parameter and physical values for the FEM explicit solution from the FSI panel (see paragraph IV.3.4);

• the allocation and initialization of all arrays needed for FEM calculation;

• the building-up of data superstructure to enable parallel computing; • the calculation of the stiness matrix for each element of the faceted surfaces of deformable thin-walled components;

• the storing of the starting reference conguration (nodes normal vectors and rotations);

• the calculation of the ELDS time step ∆tELDS and number of FEM subiterations nF EM (see equation III.40).

CHAPTER III.

ELDS statement

60

III.4


Figure III.14: ELDS initialization operations

As stated earlier, the ELDS computing task foresees that the operations reported in gure III.15 have to be performed in sequence. During the FSI computational analysis, this sequence is iteratively performed within a generic time step of the CFD calculation for a total number of times (FEM subtimesteps) nF EM given by the following ratio:

nF EM =

∆tCF D . ∆tELDS

(III.40)

To properly perform the computing task, the ELDS needs, as input, the data of loads at the elements nodes of the deformable structures, and provides, as output, the new position of nodes after every single CFD time step of calculation. Another fundamental assumption is that this latter timestep is supposed to be constant.

CHAPTER III.

ELDS statement

61

III.4


Figure III.15: ELDS computing operations

For each sub-timestep of the FEM solution and for each deformable component, the following calculations are performed in sequence:

• the local reference system (Cartesian orthogonal) for each triangle of the surface mesh (as described in section III.3.1);

• the rotation matrix for each element with respect to the global reference system;

• the nodal uid loads in the global reference system; • the mass of each node and nodal force due to gravity in the global reference system (as described in section III.3.4);

CHAPTER III.

ELDS statement

62

III.4


• the normal vector and rotations for each node in the local reference (as described in section III.3.4);

• the total displacement for each node and the component due to elastic deformation;

• the nodal elastic force vector (as described in section III.3.1); • the shift from local to global reference frame; • the integration of Newton's equation according to the procedure described in section III.3.2;

• the nodes position updating. The interaction of ELDS with the CFD solution procedure within FLUENT is explained in the next chapter.

CHAPTER III.

ELDS statement

63

Chapter IV Coupling of ELDS with FLUENT In this chapter the methodological procedure to embed ELDS in the framework of analysis of FLUENT is detailed. The rst three sections are concerned with the CFD commercial suite. In particular, section IV.1 reports a brief introduction on it, whereas sections IV.2 and IV.3 respectively describe the process that is needed to be carried out in order to perform a standard CFD analysis and how this process is modied by the introduction of ELDS so as to enable FSI studies. Finally, sections IV.4 and IV.5 respectively illustrate the dynamic mesh methods supplied by the commercial tool, and an outlook on the whole ELDS practical use.

IV.1 Introduction to FLUENT suite FLUENT [1] is one of the most known and used CFD commercial software throughout the world both in the industrial and research sector. This numerical suite is currently developed, released, and distributed by ANSYS Inc. and allows the simulation of the uid ow, heat and mass transfer in complex systems as well as related physical phenomena involving turbulence and chemical reactions. It provides the user with a wide range of models to perform steady-state or transient analyses dealing with lots of types of uids

64

IV.2

Standard CFD analysis procedure

and conditions including either incompressible or compressible ows, either laminar or turbulent regimes, either Newtonian or non-Newtonian uids, either ideal or real gases, as well as multi-phase ows. Its nite volume solver can run either serial or parallel computing with lots of solution options on structured, unstructured and hybrid meshes made by almost any combination of mesh types comprising hanging nodes and non-matching mesh interfaces. Moreover, it supplies a solution-based mesh adaption and can handle moving walls exploiting algorithms for mesh smoothing and remeshing. Since it is essentially a CFD solver, in view of carrying out an FSI analysis the recurrence either to third party code or inecient transferring data method is mandatory. ELDS has been conceived and developed to overcome this limit, and the features as well as manners that have been adopted to link the structural solver algorithm to the frame of the CFD one are treated hereafter. In the following section the description of the principal phases of a baseline CFD analysis performed by means of FLUENT is reported.

IV.2 Standard CFD analysis procedure According to the commonly recognized subdivision, three steps are usually required to be executed in sequence when performing a standard CFD analysis. These steps respectively are termed as: 1. pre-processing; 2. computing; 3. post-processing. In turn, the pre-processing and computing can be thought as divided into several subsequent substeps. Specically, in the case where FLUENT is used as solver, the gure IV.1 shows the ow chart of the operating sequence of

CHAPTER IV.

Coupling of ELDS with FLUENT

65

IV.2


actions to be performed to develop a standard transient CFD study, grouping through a dierent type of dashed boxes those referring to each of the afore-cited principal phases.

Figure IV.1: Flow chart of the standard CFD unsteady analysis

In particular, the pre-processing step foresees at rst the acquisition of the CAD (Computer Aided Design) geometrical model either through the import or (alternatively) the direct modelling of the geometries involved in the computational study, the extraction of the wetted surfaces from this model, and their cleaning-up. In fact, in order to gain the so-called cleaned-up model, that is a well-suited model for meshing and computing, a series of modications of the wetted surfaces are usually necessary. In general, both the degree and typology of these adjustments are strictly dependent on the requirements imposed by the specic code used for mesh generation that, in turn, is typically related to the solver used for solution computing. The cleaned-up model supports the subsequent generation of the surface and volume mesh through, respectively, the discretisation of the wetted surfaces into surface elements and of the volume enclosed into the calculation domain into cells.

CHAPTER IV.


66

IV.2


Once the typology of the analysis, the physical models, the parameters values, and the solution schemes have been imposed, the computational model is nally ready for the computing step. Lastly, when the numerical solution is achieved and hence computational outputs are stored and available, during the post-processing phase the results concerning the ow eld can be visualized, analyzed, and appropriately manipulated according to the peculiar needs of the analysis. In following subsections for each of the afore-listed steps, the actions that are specically required to be accomplished to suitably perform a CFD study by means of FLUENT are detailed.

IV.2.1 Pre-processing stage The pre-processing ranges from the acquisition of the geometrical model to the setting-up of the computational model. The following subsections deepen the operations previously pointed out.

CAD modelling/import In most cases, the geometries involved in a CFD study are supposed to accurately represent the objects being simulated. As such, the shapes of these objects need to be suitably modelled together with additional surfaces that are often necessary for the completion of the geometrical model. In case of external aerodynamic studies, for instance, the creation of the surfaces constituting the simulation volume is mandatory in order to enclose the virtual domain of space where the numerical solution is computed, that is the ow eld is simulated. The shapes of interest, particularly in case of complex assemblies, are usually generated by means of solid modellers able to eciently create CAD entities such as NURBS (Non-Uniform Rational B-Spline) curves and surfaces, and solid body via three-dimensional construction modelling features. Figure IV.2 and IV.3 refer to, for instance, the geometrical model of two com-

CHAPTER IV.


67

IV.2


ponents of the mechanical system concerning the rst industrial test case, subsequently presented in paragraph VI.2 of chapter VI. In particular, these components are the diuser and the reed valves housing, and the geometrical model of each of them (illustrated on the right) has been generated according to the corresponding mechanical drawing (illustrated on the left) by means of a three-dimensional parametric CAD modeller.

Figure IV.2: Geometrical model of the diuser

Figure IV.3: Geometrical model of the reed valves housing

Once the system to study has been modelled, its surfaces are then exported into a format suitable for cleaning-up phase. In general, two of the most desirable formats are IGES (Initial Graphics Exchange Specication) and STEP (Standard for the Exchange of Product Model Data). In fact, since these formats maintain the mathematical denition of their geometrical entities, they allow the user to utilize plenty of options in surfaces' modications

CHAPTER IV.


68

IV.2


during cleaning-up. Alternatively, another practicable way is to directly receive the surface mesh in a variety of formats as unstructured triangulated surfaces, that is an unstructured collection of triangles dened by an unstructured collection of vertices, e.g. STL (Stereolithography) and VRML (Virtual Reality Markup Language). In this latter case, however, the process usually becomes more elaborate and less ecient due to the more limited number of modications that can be operated to the surface mesh. In the following subsection the cleaning-up of wetted surfaces of the geometrical model is deepened. The attention is only focused on this scenario because it is the most commonly met in the advanced CFD design.

Wetted surfaces cleaning-up As previously stated, the wetted surfaces belong to the geometrical model and are dened as the boundaries of the objects around which the ow is supposed to develop. As these shapes identify the domain for surface mesh generation, they are expected to adequately represent the virtual skin of the actual components according to the requirements of the algorithm utilised for meshing. Yet, unfortunately, in most cases they turn out to be inappropriate in some way to surface grid generation and, for this reason, they require some work before being discretised into elements. It is worth mentioning that the prevalence of such geometry problems has been the motivation of the investments in the development of dedicated geometry repair applications and tool in last decades. Referring to FLUENT, the primary demands to generate an appropriate mesh for calculation are the avoidance of surfaces holes, overlaps or intersections between surfaces. Furthermore, the edges among surface patches shall be pasted and free edges can appear in the model only in the case that free surfaces are present. The latter is the case of the sails of the racing yacht model shown in gure IV.4, which refers to the second industrial test case, hereafter presented in paragraph VI.2 of chapter VI.

CHAPTER IV.


69

IV.2


Figure IV.4: Cleaned-up model of the second industrial test case

In view of achieving the cleaned-up model, that is a suitable model to be submitted to surface mesh generation, the wetted surfaces shall be selected from the geometrical model, the surfaces gaps shall be closed, the geometrical details considered useless for simulation (de-featuring) as well as the duplicate or self-intersecting surfaces shall be deleted, the missing geometries shall be built, and the surface elements normals shall be made consistent. Successively, the cleaned-up model surfaces are arranged in groups in order to satisfy the specic needs of the volume mesh generation, so as to identify for instance the area for prism layers extrusion or to dene the domains for volume mesh generation, as well as more easily manage the assignment of dierent types of boundary conditions during the computational model set-up.

Surface and volume meshing In this substep of the pre-processing, the surfaces constituting the cleaned-up model are discretised into elements and then volume domains into cells in

CHAPTER IV.


70

IV.2


such a manner to cope with the specic requirements of the CFD computing. In particular, as regards surface mesh, its resolution is expected to rstly guarantee the adequate representation of the real components shapes. Secondly, it shall satisfy at least a threshold quality level. In fact, since the volume mesh is generated starting from surface mesh, the quality of cells adjacent to surfaces is strictly connected to surface elements one. Relating to volume mesh, among others, the requests concern the satisfaction of mesh criteria in terms of proper connectivity between elements, element quality, cell aspect ratio, and cell-to-cell expansion factor. Figure IV.5 and gure IV.6, for example, depict a detail of the mesh of the rst industrial test case hereafter described in paragraph VI.1 of chapter VI. Specically, they respectively refer to the surface mesh of the petals and basement of a reed valve and the longitudinal cutting plane of the tetrahedrons volume mesh.

Figure IV.5: Surface mesh detail of the rst industrial test case

Further requirements may concern the achievement of a determined resolution level in sensitive zones to meet the demands of physical models and to properly capture specic key ow features, such as ow separation and re-attachment, for example, as well as the saving of cells to decrease the

CHAPTER IV.


71

IV.2


computing resources demand. As far as mesh quality is concerned, this is typically quantied in terms of skewness in the FLUENT community and it inuences the numerical stability and accuracy in solution computing.

CFD model set-up The last phase of the pre-processing is the computational model set-up that envisages the choice of the solution type, the uid(s) physical parameters values, the boundary conditions, the contextual models, and the resolution schemes as well as the denition of monitors entities and customized utility parameters. Once all these inputs are supplied, the computational model is ready for computing.

IV.2.2 Computing stage The computing stage foresees two subsequent substeps, namely the activation of the starting ow eld and solution calculation. Both operations may be accomplished in serial or parallel manner. In turn, the rst substep can be

Figure IV.6: Volume mesh detail of the rst industrial test case

CHAPTER IV.


72

IV.2


gained either through an initialization of the computational domain or alternatively by acquiring CFD data through the loading of a result le available through a previous calculation.

IV.2.3 Post-processing stage The post-processing tool provided by FLUENT can be used to generate meaningful graphics, animations and reports that make it easy to convey CFD results. Shaded and transparent surfaces, pathlines, vector plots, contour plots, custom eld variable denition and scene construction are some among the available post-processing features. Solution data can also be exported to several third party graphics packages for additional analysis. Figure IV.7 and gure IV.8 illustrate, qualitatively, the computational results referring to the second industrial test case afterwards described in paragraph VI.2 of chapter VI. In particular, these numerical outputs are the pressure distribution on the sails of a high-performance yacht and the path lines of ow coloured by magnitude of velocity around them in steady condition.

Figure IV.7: Sails pressure distribution in steady ow condition

It is worth highlighting that the operations described for the computing and post-processing phases, are performed within FLUENT workbench through

CHAPTER IV.


73

IV.3

Customization of FLUENT computing

Figure IV.8: Stationary ow path lines around sails

either GUI or TUI (Text User Interface), whereas the previous ones, referring to pre-processing, shall be carried out by means of other codes.

IV.3 Customization of FLUENT computing This section reports the main aspects of the achievement of the CFD solution by means of FLUENT and the manner according to which ELDS is introduced to allow the development of an FSI study.

IV.3.1 FLUENT solver The governing equations of continuum mechanics, representing the kinematic and mechanical behaviour of general bodies, are commonly referred to as conservation laws. These relationships are derived by invoking the conservation of mass (continuity equation), energy, and momentum (Newton's law). Whilst they are equally applicable to solids and uids, their diering behaviour is accounted for through the use of a dierent constitutive equations. Expressed in its most general manner, the continuity equation is valid for incompressible as well as compressible ow, and in dierential form can be

CHAPTER IV.


74

IV.3


written as follows

∂ρ ~ + ∇ · (ρ~v) = Sm (IV.1) ∂t where ρ is uid density, ~v is the velocity vector, Sm is the source term, and ~ is the operator referred to as grad, that is a vector dened to be ∇ ˆ ~ = ∂ î + ∂ ˆj + ∂ k ∇ ∂x ∂y ∂z ˆ are the unit vectors respectively along the global reference where î, ˆj, and k directions x, y , and z . The rst term on the left side of the equation IV.1 is the time rate of change of the density (mass per unit volume), whereas the second is called convective term and describes the net mass ow across the control volume's boundaries. The conservation of momentum in vectorial dierential form in an inertial reference frame is given by the following relationship

∂(ρ~v) ~ ~ +∇ ~ · τ¯ + ρ~g + F ~ ext + ∇ · (ρ~v · ~v) = −∇p ∂t

(IV.2)

or alternatively in components

∂(ρvj ) ∂(ρvi vj ) ∂p τij + =− + + ρgi + Fext,i ∂t ∂xj ∂xi ∂xj where p is the static pressure, τ¯ is the stress tensor, ρ~g is the gravitational ~ ext is the external body force. body force acceleration, and F In Newtonian ows, namely uids for which the viscous stresses are proportional to deformation rates of the uid element, the expression of the stress tensor τ¯ is given by

2 T ~ · ~v + ∇ ~ · ~v ) + ∇ ~ · ~vI τ¯ = µ (∇ 3

(IV.3)

where µ is the dynamic viscosity of uid, I is the unit tensor, and the second term on the right hand side is the eect of the volume dilation. By substituting the relation IV.3 into the equation IV.2, the equation referred to as Navier-Stokes equation can be obtained

D~v ~ 1~ 1 ~ ~ = Fext − ∇p + ν ∇( ∇ · ~v) + ν∇2~v Dt ρ 3

CHAPTER IV.


75

IV.3


where ν = µ/ρ is the cinematic viscosity, and

D~ v Dt

is the total derivative of ~v

dened to be

∂vx ∂vy ∂vz D~v ∂~v = + vx + vy + vz Dt ∂t ∂x ∂y ∂z where vx , vy , and vz are the velocity vector components along the axes x, y , and z respectively. The scalar equation IV.1 and the vectorial equation IV.2 form a system of coupled non-linear partial dierential equations usually referred to as PDEs (Partial Dierential Equations). These equations completely model basic ows and can not be solved analytically for most engineering problems. With the increasing of ow complexity, additional equations need to be solved with the PDEs. For instance, in case heat transfer or compressibility eects need to be included in calculation, an additional equation for energy shall be taken into account. Its formulation in vectorial dierential form is the following

X ∂(ρEi ) ~ ~ ~ ~ + ∇· ~v(ρEi +p) = −∇ (kef f ∇T − hj Jj )+ τ¯ ef f +Sh (IV.4) ∂t j where kef f is the eective conductivity, T is the temperature, and J~j is the diusion of species j . The rst three terms on the right side of IV.4 respectively represent the energy transfer due to conduction, the species diusion, and the viscous dissipation. In particular, kef f is the sum of k and kt where

kt is the turbulent thermal conductivity depending on the specic turbulence model being used, Sh includes the heat of chemical reaction and any other volumetric heat dened sources, whereas Ei is internal energy dened as

Ei = hi −

p | ~v |2 + ρ 2

where hi is enthalpy and | ~v | is the absolute value of the velocity vector. As stated before, when the studied ow involves other physical phenomena such as species mixing, chemical reactions, or turbulence, additional conservation equations need to be solved as extensively described in detail in the

CHAPTER IV.


76

IV.3


FLUENT Users' manual [1]. Whatever the specic ow case is, FLUENT adopts the FVM to solve the PDEs constituting the governing equation of the uid-dynamic problem. The FVM is a class of discretisation schemes that has proven highly successful in approximating the solution of a widespread variety of conservation law systems. It is extensively used in uid mechanics, meteorology, electromagnetics, semi-conductor device simulation, models of biological processes and many other engineering areas governed by conservative systems that can be written in integral control volume form. This approach turns out to be particularly suited to CFD problems due to its capability of conserving solution quantities. Specically, this control-volume base technique consists of three fundamental passages:

• the discretisation of the computational domain into a collection of non overlapping control volumes (cells) that completely covers the computational domain, namely the virtual portion of space within which the behaviour of the ow is numerically reproduced;

• the denition of suitable conditions at computational domain boundaries;

• the integration of the governing equations on each individual control volume to build up the algebraic equations for the discrete dependent variables;

• the linearization of the discretised equations and the solution of the resultant linear system. Since the PDEs are non-linear and coupled, the solution process foresees an iterative calculation that starts from the initialized ow eld wherein the error in solution of the entire set of governing equations is lower that a specic

CHAPTER IV.


77

IV.3


threshold condition. Once the solution has been achieved, the distribution of the physical variables, represented by pressure and velocity for basic ow, within the whole calculation domain is available. Relating to computing, FLUENT allows the user to determine the solution choosing one of the following two approaches: 1. pressure-based solver; 2. density-based solver. These numerical algorithms use a similar discretisation process, but the procedure utilized to linearize and solve the discretised equations diers. According to the pressure-based approach, rstly developed to tackle incompressible ows, the velocity eld is calculated from momentum equation whilst the pressure eld is determined by solving either the pressure equation or the pressure correction equation, that is a relationship gained manipulating the continuity and momentum equations. Two pressure-based algorithms are provided by FLUENT, called segregated (decoupled) and coupled. The rst of these two dierent algorithms foresees that the PDEs are solved one after another, whereas the second envisages the solution of the entire set of equation at once. The coupled approach is recommended if a strong interdependence exists between density, momentum, and/or species, and it is generally characterized by a better solution convergence, but the memory requirement increases because of the data storing need. The segregated approach provides more exibility in solution procedure. On the other hand, the density-based solver approach was rstly intended to handle compressible ows. Also in this case the velocity eld is calculated from momentum equation, whereas the pressure and density elds are determined, respectively, from the equation of state and continuity equation. The density-based algorithm solves the governing equations simultaneously, namely coupled together. Governing equations for additional scalars are

CHAPTER IV.


78

IV.3


solved sequentially, that is segregated form the coupled main set. It is worth mentioning that, though the pressure- and density-based approaches were developed to cope with only low-speed incompressible and high-speed compressible ow respectively, they have been reformulated so as to also operate beyond the limit of the theoretical operating range of ow conditions of the original implementation. To get more exhaustive information on how the solvers operate see the FLUENT manual [1]. In the following subsection the customization of FLUENT computing to enable the coupling with ELDS is described.

IV.3.2 Solution customization by means of the UDF feature A UDF [2] is a feature provided by FLUENT conceived to enhance the standard capabilities of the tool through a suite of functions that can be personalized and dynamically linked to its solver. This suite of functions are algorithms written in C-based programming language which are activated through

DEFINE

macros, and other predened

macro functions supplied by the CFD tool, that, depending on the type, can be hooked to the solver by means of GUI or automatically after their formulation. In such a manner, the user is allowed to customize the solution by adding lots of capabilities to t specic modelling needs. For instance, among the variety of applications, it is feasible to dene non standard boundary conditions and material properties, modify solver variables, exchange data with the solver as well as improve solver models and post-processing. All physical values that are exchanged between a UDF with the solver are specied in the SI (International System of Units). UDFs can be divided in general purpose and model-specic DEFINE macros, and can be either called at predened simulation times in the computing process or executed on demand through a macro. Among the general purpose

CHAPTER IV.

DEFINE_ON_DEMAND

DEFINE


predened

macros are, for instance, the 79

IV.3


following:

•

DEFINE_EXECUTE_ON_LOADING

to operate whenever a UDF com-

piled library is loaded;

•

DEFINE_ADJUST

is a general-purpose macro that can be used to ad-

just or modify the FLUENT variables that are not passed as arguments;

•

DEFINE_EXECUTE_AT_END

to enable actions at the end of the CFD

iteration;

•

DEFINE_EXECUTE_AT_EXIT

to employ functions at the end of a

FLUENT session. Because of the specic working range of the

DEFINE

macros, the knowledge

of the context in which the UDFs are called and operate within the FLUENT's solution process is crucial to determine which data are current and available at any given time and accordingly implement the UDF. In gure IV.9 and IV.10 the ow chart of the solution procedure of FLUENT in the case that UDFs are used for a pressure-based (segregated and coupled) and density-based solver are respectively depicted. As far as concerns the UDFs' management, they can be executed as either intepreted or compiled. Both methods have advantages and disadvantages. To briey summarize the main characteristics of both, the interpreted UDFs do not need a compiler, are portable to other platforms, are restricted to the C programming language only, and can not be linked to compiled system or user libraries. On the other hand, the compiled UDFs need a compiler, execute faster than interpreted UDFs, are not restricted in the use of the C programming language, can call functions written in other languages, and can not necessarily be run as interpreted UDFs in the case that they contain specic elements of C language that the interpreter cannot handle. ELDS has been linked to FLUENT solver by employing an extensive use

CHAPTER IV.


80

IV.3


Figure IV.9: Flow chart of the segregated approach

Figure IV.10: Flow chart of the density-based approach

of compiled UDFs. This coupling between the uid-dynamic and structural code alters the standard framework of the typical transient CFD process pre-

CHAPTER IV.


81

IV.3


viously shown in gure IV.1 and enhances the intrinsic capabilities of analysis of FLUENT beyond its boundaries. The resulting process structure of the coupled approach is depicted in the ow chart shown in gure IV.11, where the dashed boxes highlight the operations due to the contribution of ELDS.

Figure IV.11: Flow chart of a transient FSI study through ELDS introduction

In particular, ELDS modies the computing stage by introducing the ELDS initialization and computing respectively. The ELDS initialization is accomplished after the CFD initialization and before the starting of the CFD calculation. It can be thought as divided into two subsequent actions. These actions are the set-up of FEM parameters through the customized FSI panel of FLUENT GUI (see paragraph IV.3.4) and the FEM initialization. The latter, in particular, is dened through a

DEFINE_ON_DEMAND

CHAPTER IV.

macro and enabled by means of a


DE-

82

IV.3

FINE_ON_DEMAND


feature.

The ELDS calculation runs after every CFD calculation timestep and foresees the execution of the explicit FEM iterations. The core of the FEM solver is automatically enabled by an

DEFINE_EXECUTE_AT_END

macro, and

calculates the new position of nodes belonging to deformable components. These data are then passed to the CFD code that updates the position of deformable parts elements according to specications prescribed in a FINE_GRID_MOTION

DE-

macro. Once the position of nodes has been updated,

a combination of smoothing and remeshing is performed on volume mesh to properly re-establish the mesh connectivity and quality. It is useful mentioning that we essentially refer to unsteady analysis because the FSI investigations of our interest are inherently transient.

IV.3.3 UDFs parallel implementation Considering the target of the research activity, ELDS has been conceived since the beginning to run in parallel not to limit the numerical investigation capabilities. As a matter of fact, the parallel version of FLUENT solver enables to obtain a solution to a large computational grid by exploiting the calculation of multiple processes each operating at the same time on a specic partition, that is a portion of the whole domain. In serial computing only two processes are generated, the Cortex, namely the FLUENT's process responsible for user-interface and graphics functions, and a single compute node process, whereas the parallel FLUENT architecture foresees the creation of the Cortex, a host process, and a set of n compute node processes labelled from 0 to n-1 for a total of n+2 running processes. According to its own internal procedure, the host process receives commands from the Cortex, and transfers commands to compute node process node−0. The compute node process node−0 is the reference compute process and is in charge of passing commands to other compute processes involved in cal-

CHAPTER IV.


83

IV.3


culation as illustrated in gure IV.12. All parallel processes have a single integer ID (as well as in serial computing) and a specic task to execute. In particular, the host process does not have a

Figure IV.12: Scheme of data transfer in FLUENT parallel environment

grid partition and then carries out printing, displaying messages, and printing to a le for instance, whereas compute node processes store and perform calculations on their portion of the mesh. Subsequently, once a synchronization with each other is accomplished, the messages to the host process through compute node process node−0 are passed. To do that, each single computing node process is linked to others by means of a communicator. Considering what afore-stated, in view of customizing the parallel calculation of FLUENT by means of the UDF feature, two topics need to be carefully accounted for. The rst topic is the awareness of the topological structure assumed by cells and faces afterwards the mesh partitioning, and the corresponding coding nomenclature. In particular, in partitioned grid there are two types of cells referred to as interior and exterior, and three types of faces referred to as interior, boundary zone, and external. Besides, each thread, that is a data structure used to take information about a boundary or cell zone, stores the data associated with its cells or faces in a set of arrays.

CHAPTER IV.


84

IV.3


The second topic is the parallel processes management. According to this, depending on the typology, the parallel process can perform a limited number of operations. As stated before, the main task of the host process, which apart from DPM (Discrete Phase Model) shared memory model (see [2]) does not contain grid data, is to interpret commands from Cortex. On the contrary, a compute node process has grid visibility and can access to its own partition data. Taking into account the afore-mentioned topic, in order to enable the UDF to work in parallel, either by developing a parallel UDF or by parallelizing a serial UDF, special macros are provided by the solver. These macros are grouped according to their type in the list below: 1. compiler directives; 2. communication directives; 3. predicates; 4. global reduction macros. Compiler directives are utilized to select which portion of the UDF algorithm is to be performed by the serial process, the host process, or the compute nodes processes. Compiler directories are more frequently used in their corresponding negate forms because most operations specied through UDF are executed by the serial solver and either the host or compute node process. Communications directives are used to exchange data between the host process and the compute node processes. In this data transferring process, the number of variables as well as their name type are specied. These types of directives do not need to be protected by compiler directives because they work according to an internal procedure. Predicates are macros supplied in order to enable logical tests, whereas global reduction macros allow to perform operations that collect data from all compute node processes, and reduce the data to a single value, or an array of

CHAPTER IV.


85

IV.3


values. These operations include global summations, global maximums and minimums, and global logicals and synchronization. Other macros concern the looping over cells and faces, the acquisition of cell and face partition ID, and message displaying. To gain further information on the specications of the UDF working in parallel framework see the UDF manual [2].

IV.3.4 Customized GUI implementation In order to provide the user with a support to quickly set up the main parameters needed to FEM solution and accomplish ELDS initialization as described in paragraph IV.3.2, a customized panel (FSI panel) in the same FLUENT user workbench GUI has been implemented through the scheme programming language. This algorithm generating the FSI panel is loaded before starting the case, namely after the

.cas le loading. Once loaded, one item dened FSI ap-

pears at the end of the items listed under the scroll-down menu model of the FLUENT GUI. When the panel is enabled, what is depicted in gure IV.13 appears. Specically, the FSI panel mainly allows to the user to:

• assign deformable components (maximum 2) by picking them through the classical selection list of the FLUENT GUI;

• dene the material properties and shell thickness for each deformable component;

• assign gravity vector components; • impose clamp boundary condition either by specifying the coordinates and tolerance values or selecting the computational model adjacent parts;

CHAPTER IV.


86

IV.4

Dynamic mesh and update methods in FLUENT

Figure IV.13: FSI panel to set up the FEM parameters

• dene monitoring points by specifying the coordinates and tolerance values;

• assign FEM solution parameters; • initialize FEM solver (by running ELDS initialization operations by means of the on-demand init function). The latter operation listed above is dened and performed through a FINE_ON_DEMAND

DE-

macro, and concerns the operations previously de-

scribed in paragraph III.4 for ELDS initialization denition.

IV.4 Dynamic mesh and update methods in FLUENT An FSI analysis imperatively implies the modication of the position of nodes belonging to deformable parts that, in turn, makes necessary the switching

CHAPTER IV.


87

IV.4

Dynamic mesh and update methods in FLUENT

on the dynamic mesh model of FLUENT to allow the mesh to properly change its distribution in space during simulation to adapt to surfaces shape evolution. In general, this feature of the commercial code can be utilized to model ows where the shape of the domain is changing in time due to prescribed or unprescribed motion on the domain boundaries. In our case, in particular, the motion is dened by the ELDS solution at the end of each time step of calculation as described in paragraph III.3.2. FLUENT provides the user with three types of methods to move the mesh in view of updating the volume mesh in the deforming regions subject to the motion dened at the boundaries. These methods are: 1. smoothing methods; 2. dynamic layering; 3. local remeshing methods. According to the smoothing approach, FLUENT provides two dierent techniques, the Laplacian and spring-based technique. The Laplacian technique is the most commonly utilized, is numerically inexpensive and represents the simplest among smoothing algorithms as well. This method adjusts the location of each mesh vertex to the geometric centre of its neighboring vertices without guaranteeing an improvement of mesh quality. As such, FLUENT adopts an improvement version, operating on deforming boundaries only, that works just in case the relocation of the vertex to the geometric centre of its neighboring vertices determines an improvement of the cell skewness. The spring-based smoothing method can be used to update any cell or face zone whose boundary is moving or deforming. Regardless of this model, the edges between any two mesh nodes are idealized as a network of interconnected springs. The starting conguration of cells, nodes, and edges constitutes the equilibrium state of the mesh and a generic displacement at a given bound-

CHAPTER IV.


88

IV.5

Outline of overall methodology process

ary node will generate a force proportional to the displacement along all the springs connected to the node. The dynamic layering is used in the prismatic mesh zones to add or remove the layers of cells adjacent to a moving boundary, based on the height of the layer adjacent to the moving surface in a specic operating condition. The get further information on the theoretical background, set-up, and use of the dynamic mesh model in FLUENT as well as its options and parameters refer to the theory manual [1].

IV.5 Outline of overall methodology process The main characteristic of the implemented approach for the solution of an FSI problem by means of the ELDS embedded into FLUENT CFD solution are summarized in table IV.1. Feature type Numerical FSI approach CFD solver Structural solver Interface data exchange Manner for data exchange

Typology Weak Commercial FVM solver FLUENT In-house FEM explicit solver ELDS Through the UDF feature At each CFD time step

Table IV.1: Characteristic parameters for the proposed methodology

To provide an outline of the basic actions required to perform a complete FSI study, the fundamental steps are listed below: 1. run the FLUENT GUI; 2. read the scheme le to load the FEM GUI; 3. prepare the CFD case; 4. load the compiled libraries (previously obtained by compiling the source code of ELDS);

CHAPTER IV.


89

IV.5

Outline of overall methodology process

5. initialize the ow eld or alternatively read a result le; 6. open the FSI panel and impose the FEM parameters of the ELDS solver; 7. initialize the FEM solution by means of the on-demand init function; 8. run CFD solution.

CHAPTER IV.


90

Chapter V ELDS validation In this chapter the development of the validation tests to verify the correctness of the ELDS computational results is reported. In sections V.1 and V.2, the preliminary test cases and a comparison between the numerical outputs gained through ELDS and the commercial FEM solver LS-DYNA on the same conguration of study are respectively detailed. Specically, the preliminary test cases are the free falling test and dependence of starting shape conguration test, whereas the conguration of study for the FEM solvers comparison foresees the reproduction of the large deection of a rectangular sheet clamped along one of its minor edges. Since this chapter is essentially devoted to the ELDS validation, all trials described are carried out by disabling ow calculation so that the whole solution frame works as a pure FEM solver.

V.1 Preliminary validation tests Two dierent preliminary validation trials have been performed. Specically, the purpose of the rst test is the simulation of the free falling of a rectangular sheet, whereas the second one mainly aims at showing the eectiveness of an important feature of ELDS, namely the capability to treat whatever starting shape of the deformable components accounting for that during FEM

91

V.1

Preliminary validation tests

calculation. In both preliminary trials the sheet is supposed to be made of paper and the physical parameters values utilized for the numerical solution are reported in the SI in table V.1. Physical parameter Young's modulus Poisson's ratio Density

Value 7.7 GP a 0.3 1350 kg/m3

Table V.1: Physical parameters values for both preliminary validation trials

In the following subsections both preliminary trials set-up and computational results are detailed and discussed.

V.1.1 Free falling test As previously stated, the rst preliminary validation test is intended to verify the proper reproduction of the free drop of an object only subjected to gravity force without the contribution of the uid resistance. Figure V.1 shows the model used for this test from a top view. The model consists of a rectangular sheet of 0.4 mm of thickness with edges length of 2 and 20 cm respectively. The model is discretised with an oriented structured triangular mesh for a total amount of 80 elements and 63 nodes. Geometrical dimensions characterizing this model are summarized in table V.2. The simulation volume is constituted by a parallelepiped box with dimen-

Figure V.1: Mesh of the sheet model of the rst preliminary validation trial

CHAPTER V.

ELDS validation

92

V.1


Parameter Length Width Thickness

Value 0.2 m 0.02 m 0.0002 m

Table V.2: Geometrical dimensions of the sheet model sions 0.3, 0.08, and 1.2 m respectively along x, y , and z axes with the major edge then parallel to vertical z -axis. Its constituting surfaces are discretised by an unstructured triangular mesh with quite constant resolution as shown in gure V.2 on the left. Since the gravity acceleration direction is toward the negative values of z -axis, the sheet component is envisaged to cover almost the entire vertical dimension of the simulation volume. As such, both the smoothing and remeshing algorithms are enabled. By disabling the visualization of lateral walls of the simulation volume, the gure V.2 depicts both the top and the bottom extremity surfaces, and the sheet representing the component of interest for the numerical calculation. Since the ow calculation is disabled in order to make ELDS work only, the typology of the boundary condition does not inuence the nal result. However, a standard wall condition was applied to all zones. By neglecting the resistance of air, the displacement at specic instant t of time z(t) for a falling component that starts at rest is given by the following formula

1 (V.1) z(t) = gt2 2 where g is the absolute value of the gravity acceleration imposed equal to 9.81 m/s2 .

To monitor the distance covered by the component during falling, the position of the mid-side node of the external edge is recorded. Figure V.3 illustrates the position of the rectangular sheet during simulation. Specically, several dierent images stored from 0.5 to 4.5 seconds of simu-

CHAPTER V.

ELDS validation

93

V.1


Figure V.2: Mesh of the parts constituting the free falling test case

lation with step 0.5 second are placed side by side and visualised from the same perspective.

Figure V.3: The shape assumed by the sheet during simulation

Figure V.4 shows the comparison between the ELDS numerical output of the vertical displacement (z -component) of the mid-side node of the external edge of the rectangular sheet and the theoretical parabolic prole dened by equation V.1. As the ELDS calculation output coincides with the theoreti-

CHAPTER V.

ELDS validation

94

V.1


Figure V.4: Comparison between ELDS and theoretical falling test prole

cal one, this trial test evidences the correctness in reproducing the physical phenomenon of the free drop of an object.

V.1.2 Initial conguration inuence verication The objective of second validation test is to evidence the capability of ELDS to manage the no-at initial conguration of deformable components. Although the surfaces having internal sharp edges can not be treated by the current release of the FEM solver, however this feature provides the proposed numerical methodology with a high degree of generality. In view of showing the operative correctness of this feature, the deformation of a rectangular convex sheet free from external forces and having the horizontal at shape as reference conguration has been monitored during time. As a matter of fact, considering this specic set-up, the sheet is expected to tend to naturally assume the starting (reference) conguration due to elastic recovery of the material. Hiding lateral surfaces of the simulation volume to ease the visualization, the model used to carry out this validation test is depicted in gure V.5 from a front and an oblique view, respectively reported on the left and on the right. To numerically apply the condition above described, the initialization

CHAPTER V.

ELDS validation

95

V.1


Figure V.5: Model of the second validation test case function of the source code has been properly modied so as to assign, at each node, the value of zero to the components of both the nodal normal and nodal rotation vector except for the vertical component of nodal normal vector. Moreover, during simulation, both the gravity force and CFD solution have been disabled. As far as concerns other aspects of the computational model, all components including deformable sheet have been discretised by means of an unstructured triangular mesh, two opposite extreme nodes lying in the valley of the convex sheet are clamped, and a wall condition was applied to lateral surfaces, and to both the top and bottom surfaces. Regarding material properties values, those reported in table V.2 have been utilized. Figure V.6 and gure V.7 show the most relevant computational outputs of this validation test. Specically, in the rst image several snapshots of the shape of the sheet recorded during simulation time are depicted, whereas in the second one the prole of the vertical coordinate of the lateral midside node of the sheet mesh versus simulation time is plotted. The sequence illustrated from a front view in the rst image, in particular, macroscopically highlights how the sheet oscillates around its reference at conguration decreasing during time the amplitude of oscillation because of the material

CHAPTER V.

ELDS validation

96

V.1


Figure V.6: Congurations of the sheet during simulation damping ξ imposed equal to 0.5. This latter consideration can also be made

Figure V.7: Vertical coordinate prole of the monitoring node

considering the prole of the second image, where the oscillation of the mon-

CHAPTER V.

ELDS validation

97

V.2

Comparison with LS-DYNA code

itoring node is reported in a range of time from 0 to 0.6 s. The maximum amplitude of the oscillation occurs at the beginning of the simulation and is about 0.09 m and the asymptotic vertical position value is about 0.05 m. Taking into account the considerations made, the validation test can be judged positively accomplished and then ELDS properly includes in calculation the starting curved conguration of the deformable components.

V.2 Comparison with LS-DYNA code This validation test concerns the analysis of the deformation of the rectangular sheet used for the free falling test. To this end, in the model used to develop the rst validation test, the nodes belonging to one of its short edges have been clamped. Then the position of the external mid-side node of the sheet has been monitored during simulation, and this computational result has been compared with the corresponding one gained through the commercial FEM explicit solver LS-DYNA. In table V.3 the value of the material physical properties and main parameters for the simulation are summarized. Variable Young's modulus Poisson's ratio Component thickness Density Gravity acceleration Time step size

Typology 7.7 GP a 0.33 0.2 mm 1350 kg/m3 9.81 m/s2 0.01 s

Table V.3: Parameters of the computational model

Figure V.8 depicts the proles of the vertical position of the monitoring node during time for dierent values of the damping coecient, namely ξ =0.5,

ξ =1.0, ξ =2.0, and ξ =3.0. The value ξ =1.0 represents the critical value that

CHAPTER V.

ELDS validation

98

V.2


makes the sheet reach its maximum elongation without oscillating. The reported proles show that the higher the damping coecient value, the higher the simulation time value at which the nal deformed conguration is reached.

Figure V.8: Monitoring point displacement z-component prole

Figure V.9 shows the comparison between serial (1 computing process) and parallel (2 and 8 computing processes) computing results in terms of the monitoring node vertical displacement. Since these proles are coincident, the distinctiveness of FEM explicit solution and its independence on the number of chosen computing processes are evidenced. The gure V.10 reports the comparison between the results gained by means of ELDS and LS-DYNA for a value of the damping coecient equal to 1.0. The results obtained through both codes compare favorably well although the ELDS prediction slightly overestimates the deection of the sheet with respect to LS-DYNA. Figure V.11 illustrates, from the side view, the sheet in the deformed nal conguration, namely at about 1 s of simulation time. The deformed conguration, in particular, shows large deections, as well as coherency in deformation.

CHAPTER V.

ELDS validation

99

V.2


Figure V.9: Comparison between serial and parallel computing results

Figure V.10: Numerical proles comparison

CHAPTER V.

ELDS validation

100

V.2


Figure V.11: Final deformed conguration of the sheet

CHAPTER V.

ELDS validation

101

Chapter VI Industrial case studies With the objective to highlight the capability in coping with industrial needs, in this chapter the employment of ELDS coupled with FLUENT to perform the FSI analysis of real meaningful applications is presented. In particular, these analyses concern the simulation of the motion of the petals of a two-stroke engine reed valve, and the blowing of sails of an America's Cup Yacht. In the following paragraphs, the simulation of the just cited applications are respectively described.

VI.1 Case study 1 The rst case study aims at simulating the dynamic of a six petals reed valve due to the interaction with the uid ow passing through the duct of a twostroke engine. Reed valves (some commercial examples are reported in gure VI.1) are often used as pressure-driven ow stopper in complex mechanical systems such as two-stroke and jet engines, and compressors. From the functional operating point of view, this mechanical system is made up by a hollow metal basement which has some openings (windows) that are partially covered by exible petals that are typically mounted on the basement through screws.

102

VI.1

Case study 1

Figure VI.1: Examples of commercial reed valve pack

During valve operation, the petals are expected to be deected by pressure from the uid ow passing through the windows. The basement acts as a barrier for the valve closure and petal movement and, in some applications, there is also a curved stop component xed similarly to petals which has the double function to prevent an excessive opening of the petals and uiddynamically optimize their shape in the open channel conguration. Since the dynamic response of the valve, determined by the interaction between the pulsating ow and petals movement, inuences the eciency of the whole system it is placed into, to get the knowledge of the reed valve timing, as well as the inuence of the main geometrical and physical parameters on its behaviour are crucial and need to be deepened for engineering design.

VI.1.1 Case study 1: model set-up The studied model, courtesy of the Industrial Engineering Department of the University of Rome "Tor Vergata", is an assembly reproducing the inlet system of Maxter 125 cc two-stroke engine used for high performance go

CHAPTER VI.

Industrial case studies

103

VI.1

Case study 1

karts. According to the path covered by uid ow, the actual components composing the system respectively are the carburettor, the diuser, the reed valves pack, and the primary inlet duct of the crankcase. The mathematical representation of the wetted surfaces of the system have been extracted from the complete model created by means of a three dimensional CAD modeller adopting a reverse engineering technique. Figure VI.2 shows, on the left, the carburettor through a mid-section cutting sketch and, on the right, its corresponding simplied representation, namely the wetted surfaces of the convergent duct where valve body surfaces have been neglected.

Figure VI.2: Generation of the model of the convergent cone of the carburettor Other components dening the system are the diuser and the reed valve housing that have been already illustrated together with the images of the corresponding two-dimensional mechanical drawing in chapter IV in gure IV.2 and gure IV.3 respectively. The shape of the diuser is properly designed to allow to gently deform the ux from the circular cross section of the carburettor to the rectangular section presented by the reed valve pack, whereas the reed valve housing is conceived to split the ow into the six windows that may be partially opened by petals movement during the transient loading of the typical working cycle of the engine. The last component of the system, namely the primary inlet duct of the crankcase, is represented in gure VI.3 and has a twofold function. At rst,

CHAPTER VI.


104

VI.1

Case study 1

it conveys the ow downstream to the six windows and, at second, it acts as valve stopper by accommodating petals at full opening position. In this latter gure, the model of the complete assembly and a detail of reed petals at rest (straight conguration) through a transparent visualization are respectively reported on the left and on the right.

Figure VI.3: Solid model of the primary inlet duct of the crankcase Once all surfaces needed to detail the system objective of the analysis have been assembled, the whole model has been cleaned-up and then discretised into triangles. The algorithm adopted to generate the most part of the surface mesh was a free generation algorithm. Figure VI.4 illustrates from an oblique view a detail of the surface discretisation of the basement in which the starting conguration of petals as well as the components connecting petals to basement are also visible. The total area of the six petals is 0.003468

m2 , they are made of the same material, and each one is 34 mm length, 17 mm width, 0.4 mm thick. Because of the proximity of petals to basement, particular attention has been taken to prevent the generation of cells with a high value of skewness nearby these connections. To this end, an advancing front type algorithm has been utilized for the petals surface mesh generation in order to achieve a ner resolution in proximity of the extreme short edges. The obtained result, that is a discretisation gradient along petals length, is well visible from a top view in gure VI.5. The nal computational grid has about 640000 triangles and 310000 tetrahedral cells. Table VI.1 summarizes the value of the main parameters of the material

CHAPTER VI.


105

VI.1

Case study 1

Figure VI.4: Detail of the surface mesh of the basement

which all petals are made of, whilst the uid is assumed to be approximated by air with a constant density of value 1.22 kg/m3 . As far as boundary conditions are concerned, a pressure condition has been imposed at inlet and outlet surfaces with a gauge pressure of 10 kP a, whereas all remaining surfaces have a standard wall condition. Regarding petals, these are clamped at nodes lying along one of the minor sides by selecting in the FSI panel their corresponding supports. The solution computing is divided into two subsequent phases: 1. phase1: steady CFD calculation by disabling the FEM solution in order to achieve the stationary ow eld;

Figure VI.5: Discretisation of the upper three petals

CHAPTER VI.


106

VI.1

Case study 1

Physical parameter Young's modulus Poisson's ratio Density Damping coecient

Value 21.2 GP a 0.31 1350 kg/m3 0.5

Table VI.1: Physical parameters values for the industrial test case 1

2. phase2: transient FSI calculation initialized by means of the steady ow eld gained at the end of the previous computing phase. To accomplish the phase2 of computing, a timestep of 0.01 ms with 20 subiterations per timestep was used by including the gravity force contribution, and by enabling both smoothing and remeshing algorithm because of the expected displacements degree of petals. In the following paragraph some numerical outputs of interest for the FSI calculation phase are reported.

VI.1.2 Case study 1: computational outputs The results described below refer to the FSI analysis and, in particular, are focused on the petals opening under the loading of air ow due to dierence of pressure established between the entrance and the exit surface of the model. Under this condition, the ow pushes petals and determines their deformation from the straight shape characterizing the initial conguration to the deformed one. As already stated, the entity of the deformation aects the eciency as well as the behaviour of the whole mechanical system. To demonstrate this inuence, the rst set of snapshots, collected in gure VI.6, shows the ow distribution for dierent instants of simulation time ranging from 0.1 to 5 ms, by visualizing the uid path lines coloured according to the velocity magnitude and tracked through spheres. Moreover, the path lines are visualized with an autoscale range of values starting from 0 m/s. This sequence of images

CHAPTER VI.


107

VI.1

Case study 1

Figure VI.6: Flow path lines for 10 kP a gauge pressure case highlights the initial deformation of the petals and the larger displacement interesting the central ones going on with the progress of the simulation. Figure VI.7 depicts, from a lateral view and in the same range of time as the previous sequence, a series of snapshots of the velocity vectors projected onto the symmetry plane and coloured by the longitudinal component of velocity. The visualization of velocity values, xed between 0 and 150 m/s, shows the alignment of vectors successively to the initial deformation of the petals. The central portion of the duct is interested by the major values of velocity that stabilizes during time.

CHAPTER VI.


108

VI.1

Case study 1

Figure VI.7: Opening of petals for 10 kP a gauge pressure case

CHAPTER VI.


109

VI.2

Case study 2

To detail the behaviour of petals, the position of the minor free edge midnodes of the central petals has been monitored and recorded during simulation. Referring to them as top and bottom, the gure VI.8 depicts the proles of the vertical component of displacement vector for both nodes respectively.

Figure VI.8: Monitoring nodes vertical displacement prole

In this latter gure, the transient proles of displacements vector vertical component of monitoring nodes are plotted. It can be seen that, because of constant gauge pressure condition, proles tend to assume an asymptotic value of about 4 mm.

VI.2 Case study 2 The second industrial case study is concerned with the nautical eld and, in particular, with an old conguration America's Cup Class Yacht. In this sector in the last decade almost all competing teams have invested lots of economical resources and technical eorts in the CFD simulation-based design because of its ever-growing reliability in determining sail aerodynamic indices, hull and appendages hydrodynamic, as well as the eects of small

CHAPTER VI.


110

VI.2

Case study 2

changes on sailing performances. Besides, with respect to wind tunnel trials on scaled model of the boat, the virtual computational analyses are capable of providing an extensive range of information through the numerical solution and, by controlling all characteristic parameters at the same time, to predict the behaviour of the boat even in unusual sailing condition scenarios. The analyzed racing yacht model, already shown in gure IV.4, is by courtesy of Ignazio Viola and represents a coarser version of that employed in 2009 to break the limit of 1 billion computational cells using a software [46]. Specifically, this model reproduces the Luna Rossa scaled wind tunnel model [47] sailed with the mainsail and asymmetrical spinnaker referring to a downwind sailed conguration of study. The FSI analysis, in this case, principally aims at investigating the sails swelling up due to light wind at 45◦ apparent wind angle. As most data are condential, parameters' values are not reported and the results are illustrated qualitatively only.

VI.2.1 Case study 2: model set-up Most of the geometrical parts of the simulation volume and the Luna Rossa model are reported in gure VI.9. As far as boundary conditions are concerned, the same adopted in the previous study [47] were used for the FSI computing. As for the computational grid, a tetrahedral volume mesh was generated starting from one single oriented domain exclusively discretised with triangles. The k- turbulence model has been employed whereas both the dimensions of the sails and the values of the physical properties of the material they are made of are condential. As regards structural constraints, the main sail has both the vertex and the nodes shared with the mainmast clamped, whereas the spinnaker has the three extreme vertexes clamped. These constraints are respectively high-

CHAPTER VI.


111

VI.2

Case study 2

Figure VI.9: Parts constituting the simulation volume

lighted in gure VI.10 where the sails are represented according to their conguration in the native model. As far as boundary conditions are concerned, a velocity inlet and a pressure outlet condition has been imposed respectively at inlet and outlet surface. All surfaces of the model of the yacht as well as the sails have a standard wall condition, whereas the remaining surfaces have a wall condition without resistance (free slip). Referring to uid, the air with constant density of value 1.22 kg/m3 has been considered. The computing study has been divided into the three subsequent phases as follows: 1. phase1: a pure explicit FEM simulation by disabling the CFD solution was performed to make the sails go soft under the only eect of gravity force to obtain the shape of sails without wind; 2. phase2: steady CFD calculation by disabling the FEM solution was carried out to reach the starting ow eld corresponding to the starting ow conguration of the FSI analysis;

CHAPTER VI.


112

VI.2

Case study 2

Figure VI.10: Clamped nodes of main sail (left) and spinnaker (right)

3. phase3: transient FSI calculation initialized by means of the steady ow eld gained at the end of the previous phase. In the computing phases where the FEM calculation is enabled, that is phase1 and phase3, the gravity force is included in calculation. In the phase1 the timestep used was of 0.05 s, whilst in the phase3 a timestep of 0.5 ms with 24 subiterations per timestep was used. The gure VI.11 shows a sequence of several snapshots of the shape of the sails during the rst phase of computing. In particular, the sequence regards the congurations ranging from 0.05 and 1.2 s of simulation time, and evidences the fall of sails coherently with the structural constraints. Since during the sails blowing a large displacement of sails is expected, both the smoothing and remeshing algorithm have been enabled to allow the accomplishment of the solution.

VI.2.2 Case study 2: computational outputs The gure VI.12 depicts, for dierent instants of the FSI computing, the shape of sails coloured by pressure xing the minimum and maximum for

CHAPTER VI.


113

VI.2

Case study 2

Figure VI.11: The sails blowing down sequence visualization. Due to its extensive area and since its shape is expected to sensitively change with respect to the starting conguration, the spinnaker shows more explicitly the geometrical modication consequent to the air blowing. In particular, in the rst part of the simulation, the blowing up mainly interests the bottom area of the sail, and then gradually involves the remaining part of the sail till it is completely blown up. In this conguration the pressure assumes a steady distribution and the pressure contours become quite concentric.

CHAPTER VI.


114

VI.2

Case study 2

Figure VI.12: The sails blowing sequence On the contrary, at least for this type of computational results, the displacements of the main sail are more dicult to be identied. Relating to other achievable results, some of the most useful yacht design parameters are the proles of total force components acting on the sails. Such

CHAPTER VI.


115

VI.2

Case study 2

data are reported in gure VI.13 and gure VI.14 for the main and spinnaker respectively. Specically, for both, the components along the three axes of

Figure VI.13: Proles of total force components acting on main sail

Figure VI.14: Proles of total force components acting on spinnaker

the global reference system of the total net force are detailed from 0 to 0.4 s of simulation time. The corresponding components of these proles have qualitatively the same behaviour. The Fy component, for instance, presents a low value until about

CHAPTER VI.


116

VI.2

Case study 2

0.035 s, and afterwards has a sudden increment in term of absolute value. Similarly acts the Fx component, whereas the Fz component oscillates around lower values than others components.

CHAPTER VI.


117

Chapter VII Conclusions The development and use of an explicit FEM solver to enable three-dimensional FSI studies of deformable isotropic shell structures by means of the commercial CFD software FLUENT has been presented. This algorithm, called ELDS, has been written in a C-based language and can be dynamically embedded into the CFD commercial solver by means of the UDF feature. The developed FEM solver is based on the co-rotational formulation and can be used to numerically reproduce large displacements of deformable structures under the assumptions of small strains eld. To this end, a general three-noded triangular element has been implemented to deal with whatever shape of deformable components. Besides, in order to eciently support the user in the setting up of the study, an intuitive GUI within the same workbench framework of FLUENT has been built as well. The theoretical background and the methodologies adopted to implement the structural solver have been described, as well as the validation cases to evaluate its correctness. Several test cases have been performed and detailed to show the eectiveness in handling and solving real applications in dierent scientic areas. The challenging trait of these test cases evidences the maturity reached by ELDS and its eectiveness in providing solutions to complex problems in the real

118

world of technique and scientic design. Since it has been implemented in parallel, ELDS can be used to carry out FSI studies of small models, as well as of large ones with millions of cells. Potential improvements to be developed in next future deal with the accomplishment of further validation tests, the formulation of other types of element, and the implementation of a contact algorithm.

CHAPTER VII.

Conclusions

119

Bibliography [1] ANSYS, (2009),

ANSYS FLUENT User's Guide - Release 12, ANSYS

Inc. [2] ANSYS, (2009),

UDF User's Guide - Release 12, ANSYS Inc.

Stabilized Pressure Segregation Methods and Their Application to Fluid-Structure Interaction Problems, PhD thesis, Tech-

[3] Badia, S. I., (2006),

nical University of Cataluña.

Fluid-structure partitioned procedures based on robin transmissions conditions, Journal of

[4] Badia, S., Nobile, F., and Vergara, C., (2007), Computational Physics, submitted. [5] Bathe, K. J., (1996),

Finite Elements Procedures, Prentice-Hall, New

Jersey. [6] Bathe, K.J., Ramm E., and Wilson E.L., (1975),

Finite element for-

mulations for large deformation dynamic analysis, Int. J. Num. Method Engng., 9, 353-386. [7] Batoz, J. L., Bathe K. J., and Ho L. W., (1980), A Study of Three-Noded

Triangular Plate Bending Elements, International Journal for Numerical Methods in Engineering, Vol.15, 1771-1812. [8] Batoz, J. L., (1982),

An explicit formulation for an ecient triangular

plate-bending element, International Journal for Numerical Methods in Engineering, 18:1077-1089. 120

BIBLIOGRAPHY

Formulation and evaluation of new triangular, quadrilateral, pentagonal and hexagonal discrete Kirchho plate/shell elements, International Journal for Numerical Methods in Engineering,

[9] Batoz, J. L., (2001),

52:615-630. [10] Belytschko, T., and Hsieh B. J., (1997),

Non-linear transient nite ele-

ment analysis with convected coordinates, International Journal for Numerical Methods in Engineering, 7:255-271. [11] Belingardi, G., (1998),

Il Metodo degli Elementi Finiti nella Proget-

tazione Meccanica, Levrotto&Bella. [12] Biancolini M, (2011),

Analisi Fluido Struttura di una Ala Anteriore di

Formula 1, XL Convegno Nazionale AIAS, Palermo. Fluid Structure Interaction (FSI) with RBF Morph: a Generic Formula 1 Front End, EnginSoft International Con-

[13] Biancolini M, (2011), ference 2011.

[14] Baudille, R. and Biancolini, M. E. (2002),

Dynamic analysis of a two

stroke engine reed valve, Proceedings on CD-ROM of XXXI AIAS Conference, Parma, Italy.

Modelling FSI problems in FLUENT: a dedicated approach by means of UDF programming, Pro-

[15] Baudille, R. and Biancolini, M. E., (2005),

ceedings on CD-ROM of the European Automotive CFD Conference, Frankfurt, Germany. [16] Baudille, R. and Biancolini, M. E., (2008), A general approach for study-

ing the motion of a cantilever beam interacting with a 2D uid ow, Interaction and Multiscale Mechanics, Vol. 1, No. 4.

BIBLIOGRAPHY

121

BIBLIOGRAPHY [17] Baudille, R., Biancolini M. E., and Mottola E.,(2009), Non-linear models

of reed valve dynamics, Int. J. Vehicle Systems Modelling and Testing, Vol. 4, No. 3.

Added-mass eect in the design of partitioned algorithms for uid-structure problems, Com-

[18] Causin, P., Gerbeau, J. F., and Nobile, F., (2005),

puter Methods in Applied Mechanics and Engineering, 194:4506-4527. [19] Cook, R. D., Malkus D. S., and Plesha M., (1989),

Concepts and Appli-

cations of Finite Element Analysis, Third edition, John Wiley & Sons. [20] Courant, R., Friedrichs K., and Lewy H., (1928),

Über die partiellen

Dierenzengleichungen der mathematischen Physik, Mathematische Annalen, vol. 100, no. 1, pages 32-74. [21] Courant, R., Friedrichs K., and Lewy H.,(1967), On the partial dierence

equations of mathematical physics, IBM Journal, pp. 215-234, English translation of the 1928 German original.

Finite Element Modelling of Fluid Flow with Moving Free Surfaces and Interfaces Including Fluid-Solid Interaction,

[22] Dettmer, W. G., (2004),

PhD thesis, School of Engineering, University of Wales, Swansea. [23] Farhat, C., (1997), Parallel and distribuited solution of coupled nonlinear

dynamic aeroelastic response problems, Papadrakakis M. editor, Parallel Solution Methods in Computational Mechanics, pages 243-301. J., Wiley and Sons. [24] Felippa, C. A., Park, K. C., and Farhat, C., (1998),

Partitioned anal-

ysis of coupled systems, Oñate, E. and Idelsohn, S., editors, IV World Congress on Computational Mechanics, Barcelona, CIMNE.

BIBLIOGRAPHY

122

BIBLIOGRAPHY

A newton method using exact jacobians for solving uid-structure coupling, Computers and

[25] Fernández, M. A., and Moubachir, M., (2005), Structures, 83:127-142.

[26] Figueroa, C., Vignon-Clementel, I., Jansen, K., Hughes, T., and Taylor,

A coupled momentum method for modeling blood ow in three-dimensional deformable arteries, Computer Methods in Applied C., (2006),

Mechanics and Engineering, 195:5685-5706.

Geometrically Nonlinear Analysis of Shell Structure Using a Flat Triangular Shell Finite Element, Archives of

[27] Gal, E., and Levy R., (2006),

Computational Methods in Engineering, Vol.13, 3, 331-388.

An ecient solver for the fully coupled solution of large-displacement uidstructure interaction problems, Comput. Meth-

[28] Heil, M., (2004),

ods Appl. Mech. Engrg., 193:1-23. [29] Hou, G., Wang J., and Layton A., (2012),

NumericalMethods for Fluid-

Structure Interaction - A Review, Commun. Comput. Phys., Vol. 12, No. 2, pp. 337-377, August. [30] Kundu, P. K., and Cohen, I. M., (2008),

Fluid Mechanics (4th revised

ed.), Academic Press. [31] Jeyachandrabose, C., Kirkhope, J., and Babu, C. R., (1985),

An alter-

native explicit formulation for the DKT plate-bending element, International Journal for Numerical Methods in Engineering, Volume 21, Issue 7, pages 1289-1293. [32] Logan, Daryl L., (2012),

A First Course in the Finite Element Method,

Global Engineering. [33] Mohr, G. A., (1992), Finite elements for solids,

uids, and optimization,

Oxford Science Publications. BIBLIOGRAPHY

123

BIBLIOGRAPHY [34] Matthies, H. G., and Steindorf J., (2004),

Partitioned Strong Coupling

Algorithms for Fluid-Structure-Interaction, Institute of Scientic Computing, Technical University Braunschweig, Brunswick, Germany.

Numerical Approximation of Fluid-Structure Interaction Problems with Application to Haemo-dynamics, PhD thesis,

[35] Nobile, F. (2001),

D'epartement de Mathématiques, École Polytechnique Fédérale de Laussane.

Partitioned analysis schemes for the transient interaction of incompressible ows and nonlinear exible structures, Wall, W. A., Bletzinger K. U. and Schweizerhof K. ed-

[36] Mok, D. P. and Wall, W. A., (2001),

itors, Trends in Computational Structural Mechanics, pages 689-698, Barcelona, CIMNE.

An eective uid-structure interaction formulation for vascular dynamics by generalized Robin conditions,

[37] Nobile, F., and Vergara C., (2008),

SIAM J. Sci. Comp., 30(2):731-763. [38] Park, K. C., and Felippa, C. A., (1983),

Partitioned analysis of coupled

systems, Belytschko, T. and Hughes, T. J. R. editors, Computational Methods for Transient Analysis, pages 157-219, Elsevier Science Publishers. [39] Rossi, R., (2005), Light-weight

Structures. Numerical Analysis and Cou-

pling Issues, PhD thesis, University of Bologna, Italy. [40] Rugonyi, S., and Bathe, K. J., (2001), On nite elements analysis of uid

ows coupled with structural interaction, CMES - Computer Modeling in Engineering and Sciences, 2:195-212.

Partitionerte Verfahren für Probleme der FluidStruktur Wechselwirkung, PhD thesis, Technischen Universität Braun-

[41] Steindorf, J., (2002), schweig. BIBLIOGRAPHY

124

BIBLIOGRAPHY [42] Tallec, P. L., and Mouro J., (2001), Fluid structure interaction with large

structural displacements, Computer Methods in Applied Mechanics and Engineering, 190:3039-3067.

Solution techniques for the fully discretized equations in computation of uid-structure interactions with spacetime formulations, Computer Methods in Applied

[43] Tezduyar, T. E., Sathe, S., and Stein, K. (2006),

Mechanics and Engineering, 195:5743-5753. [44] Timoshenko, S., and Woinowsky-Krieger S., (1959),

Theory of Plates

and Shells, McGraw-Hill International Editions. [45] Turner, M. J.,Clough R. W., Martin H. C., and Topp L. J., (1956),

Stiness and Deection Analysis of Complex Structures, J. Aero. Sci., Vol. 23, No. 9, pp. 805-823. [46] Viola, I. M., Ponzini, R., and Passoni, G., (2009),Sailing

Past a Billion,

ANSYS Advantage, Volume III, Issue 2. [47] Viola, I. M., (2009), Downwind Sail Aerodynamics:

A CFD Investigation

with High Grid Resolution, Journal of Ocean Engineering. [48] Wall, W. A., (1999),

Fluid-Struktur-Interaktion mit stabilisierten

Finiten Elementen, PhD thesis, Institut für Baustatik, Universität Stuttgart.

A strong coupling partitioned approach for uid-structure interaction with free surfaces,

[49] Wall, W. A., Genkinger, S., and Ramm, E., (2007), Computers and Fluids, 36:169-183.

A coupled uid structure environment with a three-dimensional shell model, ECCOMAS 2000, Barcelona,

[50] Wall, W. A. and Ramm, E., (2000), CIMNE.

BIBLIOGRAPHY

125

BIBLIOGRAPHY

Partitioned analysis schemes for the transient interaction of incompressible ows and nonlinear exible structures, Trends in computational structural mechanics, Barcelona,

[51] Wall, W. A. and Mok, D. P., (2001),

CIMNE.

A strong coupling partitioned approach for uid-structure interaction with free surfaces,

[52] Wall, W. A., Genkinger, S., and Ramm, E., (2007), Computers and Fluids, 36:169-183. [53] Weaver, W. Jr., and Johnston P. R., (1984),

Finite Elements for Struc-

tural Analysis, Prentice-Hall. [54] Wempner, G., (1969), Finite

elements, nite rotations and small strains of exible shells, International Journal of Solids and Structures, 5:117153.

[55] Wuchner, R., and Bletzinger K. U., (2005),

Stress-adapted numerical form nding of pre-stressed surfaces by the updated reference strategy, International Journal for Numerical Methods in Engineering, 64:143166.

[56] Wuchner, R., (2006),

Mechanik und Numerik der Formndung und Fluid-Struktur-Interaktion von Membrantragwerken, PhD thesis, Tech-

nischen Universitt München. [57] Zienkiewicz, O. C., (1971),

The nite element method in engineering

science, McGraw Hill.

BIBLIOGRAPHY

126

Index computational, i, ii, 14, 10, 18, 39,

acceleration, 16, 17, 38, 39, 75, 93,

42, 50, 61, 66, 67, 70, 72,

98 algorithm, 5, 10, 11, 13, 15, 18,

73, 77, 83, 86, 91, 92, 96,

34, 38, 40, 55, 65, 69, 78,

98, 105, 107, 111, 113, 115

79, 85, 86, 88, 93, 105, 107,

computing, 19, 26, 6062, 6567, 7174, 78, 79, 8284, 99,

113, 118, 119

100, 106, 111, 113

approach, ii, 1, 3, 4, 6, 915, 18, 25, 26, 3336, 5558, 77

constraint, 1, 29, 52, 111, 113

79, 81, 82, 88, 89

coupled, 2, 710, 13, 64, 7680, 82, 102

bending, 21, 41, 4648, 50

Courant, 39

boundary (condition), 7, 28, 70, 72,

CSD, 1, 6, 7, 15

79, 86, 93, 106, 111, 112

CST, 33, 4043, 46, 56 curvature, 4850

C (programming language), i, 79,

customization, 74, 79

80, 118 CAD, 6668, 104

deformable, i, 3, 7, 14, 16, 34, 36,

CAE, 1, 2, 6

38, 39, 5658, 6062, 83,

Cartesian, 23, 28, 62

86, 87, 91, 95, 96, 98, 118

CFD, i, ii, 1, 37, 33, 61, 6367,

deformation, 79, 16, 33, 3538, 47,

69, 7174, 77, 7983, 89,

63, 75, 95, 98, 99, 107, 108

96, 106, 110, 112, 118, 132

density, 75, 7881, 92, 98, 106, 107,

CFL, 39

112

clean-up, 66, 6870, 105

displacement, ii, 3, 4, 9, 13, 16

co-rotational, i, ii, 3436, 58, 118

19, 22, 23, 25, 26, 33, 35

127

INDEX 38, 4248, 50, 51, 53, 55 58, 63, 88, 89, 93, 94, 99,

feature, i, 5, 7, 34, 58, 65, 79, 83, 84, 88, 89, 91, 95, 118

107, 108, 110, 113, 115, 118

FEM, i, 35, 1416, 19, 21, 26, 28,

DKT, 33, 40, 41, 46, 50, 51, 53, 55,

33, 34, 36, 50, 58, 6062,

56

82, 83, 86, 87, 8991, 95,

DOF, 18, 40, 42, 51

98, 99, 106, 112, 113, 118,

DPM, 85

132

DRILL, 40, 41, 56 dynamic, ii, 1, 35, 7, 9, 10, 12, 15

ow, 710, 12, 16, 34, 64, 65, 67, 69, 7179, 88, 90, 91, 93,

20, 38, 60, 64, 75, 77, 79,

102108, 112, 113

81, 8789, 102, 103, 110,

FLUENT, i, 35, 33, 34, 58, 6365,

118 elasticity, 14, 2224 ELDS, 35, 15, 17, 21, 33, 34, 36, 3841, 5558, 6065, 74, 79, 80, 82, 83, 8691, 9395, 98, 99, 102, 118, 119 element, i, 1, 3, 7, 15, 18, 19, 21, 2431, 3338, 4046, 5053, 5558, 6062, 66, 6971, 75, 83, 92, 118, 119 energy, 7, 8, 15, 20, 24, 25, 50, 74, 76

67, 69, 7274, 7780, 82 84, 8689, 102, 118, 132 uid, i, ii, 1, 3, 711, 13, 14, 16, 20, 34, 40, 62, 64, 65, 72, 74, 75, 77, 81, 92, 102104, 106, 107, 112 FORTRAN, 55 FSI, i, 314, 16, 33, 34, 58, 60, 61, 64, 65, 74, 82, 83, 86, 87, 89, 90, 102, 106, 107, 111 113, 118, 119, 132 FVM, ii, 6, 77, 89

Eulerian, 6, 10

global, 18, 28, 37, 55, 62, 75, 85, 86

explicit, i, ii, 3, 11, 15, 1721, 25,

GUI, i, 34, 74, 79, 82, 86, 89, 118

32, 34, 3840, 45, 46, 55, 57, 60, 83, 89, 98, 99, 112, 114, 118 FEA, 50 INDEX

Hooke (law), 17 ID, 84, 86 IGES, 68

128

INDEX implicit, 10, 11, 15, 1719 initialization, 58, 60, 61, 73, 82, 86, 87, 90, 95 inlet, 106, 112 interpolation, 16, 25, 26, 32 isotropic, i, 3, 23, 24, 47, 118 iterative, 2, 17, 35, 60, 61, 77, 80, 83 Kirchho, 41, 46, 47, 5052 Lagrangian, 6, 7, 10, 35 linear, 14, 27, 33, 35, 36, 43, 48, 52, 77, 78 local, 18, 28, 30, 35, 37, 41, 55, 62, 63, 88 LS-DYNA, 98, 99 Luna Rossa, 111 MACRO (function), 58, 60, 79, 80, 82, 83, 8587

motion, 8, 16, 19, 26, 35, 36, 56, 88, 102 Navier-Stokes (equations), 75 Newton (equation), 16, 57, 63 Newton (uid), 65, 75 node, 17, 2632, 34, 3642, 5053, 5563, 65, 8385, 8789, 92 94, 9699, 106, 110, 111, 113, 118 non-linear, i, 9, 11, 12, 17, 19, 34 36, 76, 77 NURBS, 67 outlet, 106, 112 panel, 60, 82, 86, 87, 90, 106 parallel, i, 18, 19, 23, 29, 60, 65, 72, 8386, 99, 100, 119 Pascal (triangle), 27 PDE(s), 7678

mass, 13, 40, 57, 58, 62, 74, 75

plate, 2124, 41, 42, 4650

membrane, 8, 13, 41, 42, 46

Poisson (eect), 22

mesh, 14, 17, 39, 42, 57, 58, 62,

Poisson (ratio), 24, 92, 98, 107

6466, 6972, 83, 84, 87 89, 9294, 96, 105, 106, 111 model, 1, 58, 1013, 64, 6672,

position, 18, 28, 38, 61, 63, 83, 87, 93, 98, 105, 110 post-processing, 65, 67, 73, 79

76, 77, 79, 85, 86, 88, 89,

pre-processing, 6567, 70, 72, 74

92, 93, 95, 96, 98, 103105,

pressure, 7, 73, 75, 7880, 102, 103,

107, 111, 112, 119

106110, 112114

momentum, 74, 75, 78

INDEX

129

INDEX 7780, 83, 8693, 96, 99,

procedure, 35, 10, 13, 14, 16, 19,

105, 106, 111113, 118

21, 25, 35, 55, 6365, 78, 80, 83, 85 RBF, 14 remeshing, 65, 83, 88, 93, 107, 113 rigid, 12, 19, 26, 3537, 49, 56 rotation, 3537, 41, 42, 48, 51, 52, 5660, 62, 63

steady, 16, 19, 64, 73, 106, 107, 112114 STEP, 68 stiness, 17, 21, 25, 33, 34, 38, 40 42, 46, 51, 52, 5557, 60 STL, 69 strain, ii, 17, 2125, 33, 4143, 45, 46, 4850, 52, 53, 55, 118

scheme, 11, 12, 1719, 21, 38, 57, 60, 67, 72, 77, 84, 86, 89

stress, 7, 10, 16, 19, 2125, 34, 35, 42, 4749, 75

segregated, 7881 shape, 8, 38, 58, 67, 69, 71, 88, 91, 95, 96, 103, 104, 107, 112 114 shape function, 3033, 4446, 51, 53, 54 shear, 7, 23, 47, 50, 52 shell, 13, 3336, 39, 40, 50, 57, 86,

strong (coupling), 10, 13 temperature, 3, 7, 40, 76 thickness, 21, 22, 33, 4650, 55, 86, 92, 93, 98 TL, 35 transient, 16, 17, 19, 20, 58, 60, 64, 66, 8183, 104, 107, 110,

118 SI, 79, 92 simulation, i, 15, 9, 11, 12, 19 21, 35, 58, 60, 64, 67, 70, 77, 79, 88, 9199, 102, 107, 108, 110114, 116 skewness, 72, 88, 105

113 translation, 41 triangle, 2730, 33, 4144, 46, 50, 52, 55, 58, 62, 69, 105, 111 TUI, 74 UDF(s), i, 58, 7981, 8386, 89,

smoothing, 65, 83, 88, 93, 107, 113 solution, i, 14, 611, 13, 1519, 21, 26, 28, 34, 36, 3840,

118, 132 UL, 35 unsteady, 66, 83

57, 60, 62, 6567, 7174, INDEX

130

INDEX validation, 5, 92, 95, 96, 98, 118, 119 valve, 8, 68, 71, 102105 velocity, 10, 17, 73, 75, 76, 78, 107, 108, 112 visualization, 67, 9395, 105, 107, 108, 114 VRML, 69 weak (coupling), 10 wetted (surfaces), 66, 69, 70, 104 workbench, i Young (modulus), 24, 92, 98, 107

INDEX

131

Appendix 1 Abstract (Italian) Nella presente tesi viene descritto lo sviluppo di un'ecace ed eciente metodologia per svolgere analisi computazionali tri-dimensionali di interazione uido-struttura di componenti sottili deformabili in materiale isotropo tramite l'utilizzo del codice CFD commerciale ANSYS FLUENT. Con l'obiettivo di estendere le capacità standard di analisi del codice uidodinamico e trattare le non-linearità geometriche di strutture deformabili, grazie alla funzionalità UDF è stato implementato in un linguaggio di programmazione basato sul C un solutore FEM parallelo fondato sulla formulazione co-rotazionale. Inoltre, per supportare l'utente durante tutta la preparazione dello studio numerico FSI, è stato implementato un pannello apposito per impostare i parametri della soluzione FEM all'interno dell'interfaccia graca del solutore uido-dinamico. Diversi casi test sono stati eettuati per validare la bontà e l'adabilità dei risultati ottenuti con il solutore FEM esplicito, comparando questi con quelli di soluzioni note o di altri solutori FEM commerciali consolidati. Inne è descritto l'utilizzo della metodologia proposta nella tesi per la soluzione di casi industriali reali con la nalità di mostrare le reali potenzialità del solutore FEM e candidarlo come un maturo e robusto strumento di indagine numerica per soddisfare le sempre crescenti esigenze della progettazione ingegneristica.

132

Appendix 2 Abstract (French) Dans la présente thèse on décrit le développement d'une méthodologie de calcul ecace et eciente pour l'analyse en modélisation numérique tridimensionnelle des interactions uide structure concernant structures minces déformables en matériau isotrope utilisant le code CFD commercial ANSYS FLUENT. Dans le but d'étendre les capacités traditionnelles de l'analyse de code de dynamique des uides et de traiter la non-linéarité géométrique des composants déformables, grâce à la fonction UDF on a mis en ÷uvre, dans le langage de programmation C basé sur un parallèle FEM (Méthode des Eléments Finis), un solveur basé sur la formulation de co-rotation. En outre, pour supporter l'utilisateur lors de la préparation de l'étude numérique FSI, on a mis en place un panneau de conguration spécique permettant le réglage des paramètres de la solution aux éléments nis à l'intérieur de l'interface graphique du solveur de dynamique des uides. Un nombre considérable de cas ont été évalués an de valider la qualité et la abilité des résultats obtenus avec le solveur aux éléments nis explicite, en les comparant avec les résultats obtenus à travers des solutions connues ou d'autres solveurs FEM commerciaux consolidés. Enn, on illustre l'application de l'approche numérique proposée dans cette thèse à la solution de cas industriels réels dans le but de montrer le vrai potentiel du solveur aux éléments nis et de poser sa candidature à instrument de sondage

133

numérique complet et solide pour répondre aux besoins croissants de la conception mécanique.

Appendix 3

134

Appendix 3 Abstract (Spanish) La presente tesis describe el desarrollo de una metodología eciente y ecaz sobre el análisis computacional en tres dimensiones de la Interacción Fluido-Estructura (conocida por sus siglas en inglés como FSI) en componentes sólidos y deformables en materiales isótropos mediante la utilización del software comercial ANSYS FLUENT CFD. Con el objetivo de ampliar las capacidades de análisis estándar del código uido dinámico y tratar la no linealidad geométrica en componentes deformables, gracias a la funcionalidad UDF, se desarrolla una solución FEM paralela basada en la formulación co-rotacional, implementado en un lenguaje de programación basado en C. Adicionalmente, con la intención de facilitar la conguración del estudio FSI al usuario, se ha implantado un Interfaz Gráco de Usuario (GUI) intuitivo y fácil de usar dentro del mismo entorno de pruebas de la solución CFD. Varios casos de prueba han sido desarrollados para validar la veracidad y abilidad de los resultados obtenidos de la solución FEM. Estos resultados han sido comparados tanto con soluciones contrastadas como a través de otras soluciones FEM comercialmente consolidadas. Finalmente, se describe la aplicación de la metodología propuesta en esta tesis para la resolución de casos reales en la industria, con la nalidad de mostrar el auténtico potencial de la solución FEM, siendo ésta un claro candidato maduro y robusto para satisfacer las crecientes necesidades de diseño

135

impuestas por la ingeniería de requisitos.

Appendix 2

136