Copy No. _____ Defence Research and Development Canada

Recherche et développement pour la défense Canada

DEFENCE

&

DÉFENSE

Large Deformation Fluid-Structure Interaction and Interface Tracking for Chinook Derrick Alexander Martec Limited Prepared by: Martec Limited 1888 Brunswick Street, Suite 400 Halifax, Nova Scotia B3J 3J8

PWGSC Contract Number: W7707-115146, Task 1 Contract Scientific Authority: Malcolm J. Smith, Group Leader, NPSS, 902-426-3100 ext 383

The scientific or technical validity of this Contract Report is entirely the responsibility of the contractor and the contents do not necessarily have the approval or endorsement of Defence R&D Canada.

Defence R&D Canada – Atlantic Contract Report DRDC Atlantic CR 2011-343 February 2012

This page intentionally left blank.

Large Deformation Fluid-Structure Interaction and Interface Tracking for Chinook Derrick Alexander Martec Limited Prepared By: Martec Limited 1888 Brunswick St, Suite 400 Halifax Nova Scotia, Canada B3J 3J8 PWGSC Contract Number: W7707-115146, Task 1 CSA: Malcolm J Smith, Group Leader, NPSS, 902-426-3100 x383

The scientific or technical validity of this Contract Report is entirely the responsibility of the contractor and the contents do not necessarily have the approval or endorsement of Defence R&D Canada.

Defence R&D Canada – Atlantic Contract Report DRDC Atlantic CR 2011-343 February 2012

Principal Author Original signed by Derrick Alexander

Derrick Alexander Senior Research Scientist Reviewed by Original signed by Neil Pegg

Neil Pegg Head/ Warship Performance Section Approved by Original signed by Calvin Hyatt

Calvin Hyatt Chair/ Document Review Panel

© Her Majesty the Queen in Right of Canada, as represented by the Minister of National Defence, 2012 © Sa Majesté la Reine (en droit du Canada), telle que représentée par le ministre de la Défense nationale, 2012

Abstract …….. The ability to assess the loading, response and survivability of naval structures and equipment to withstand the effects of underwater weapons containing non-ideal explosives is of fundamental importance; particularly for determining the threat damage potential of new weapon technologies, and developing effective counter-measures. Further advancement of the Martec Limited’s Chinook computational fluid dynamics software is detailed, with more advanced methods for material interface tracking and structural large deformation capability for fluid structure interaction (FSI). The local mesh adaptation technique modifies the CFD mesh in the cells adjacent to the structure, such that the cell faces correspond to the location of the structure, allowing for precise determination of pressure loads and interface velocities. Details of the numerical implementation and theory are provided. A series of validation cases were performed in order to verify the FSI and material interface tracking functionality and validate the results. The large deformation capability, initially detailed here, can be applied for any FSI interaction, such as applications in underwater explosions, particle break-up, aerodynamic ablation, and structural failure. The large deformation capability allows for assessment of failure properties and the behaviour of structures under high-rate impulse loading.

Résumé …..... La capacité à évaluer le chargement, la réponse et la survivabilité de structures et d’équipement navals aux fins de résistance aux effets des armes sous-marines contenant des explosifs hétérogènes est d’une importance fondamentale, en particulier pour la détermination des dommages potentiels de la menace représentée par les nouvelles technologies d’armes et pour la mise au point de contre-mesures efficaces. L’avancement du logiciel Chinook de Martec Limited, qui utilise la dynamique numérique des fluides (DNF), est expliqué en détail, avec des méthodes plus avancées de suivi de l’interface des matériaux et de déformation structurale importante pour l’interaction fluides-structures (IFS). La technique d’adaptation locale de maillage modifie le maillage de DNF dans les cellules adjacentes à la structure, de manière à ce que les faces de cellule correspondent à l’emplacement de la structure, ce qui permet une détermination précise des charges de pression et des vitesses d’interface. Des détails sur l’implémentation numérique et la théorie sont fournis. Une série de validations ont été effectuées afin de vérifier les fonctions d’IFS et de suivi de l’interface des matériaux et afin de valider les résultats. La capacité de déformation importante, qui était expliquée en détail dans le présent document à l’origine, peut être appliquée pour toute IFS, comme les applications lors d’explosions sous-marines, de rupture de particules, d’ablation aérodynamique et de défaillance structurale. Elle permet l’évaluation des propriétés de défaillance et du comportement des structures sous l’effet de chargement d’impulsion à grand débit.

DRDC Atlantic CR 2011-343

i

This page intentionally left blank.

ii

DRDC Atlantic CR 2011-343

Executive summary Large Deformation Fluid-Structure Interaction and Interface Tracking for Chinook Derrick Alexander; DRDC Atlantic CR 2011-343; Defence R&D Canada – Atlantic; February 2012 Introduction or background: The ability to assess the loading, response and survivability of naval structures and equipment to withstand the effects of underwater weapons containing nonideal explosives is of fundamental importance; particularly for determining the threat damage potential of new weapon technologies, and developing effective counter-measures. Currently Martec Limited utilizes the computational fluid dynamics (CFD) code Chinook to simulate underwater explosion (UNDEX) events. Chinook has been specifically designed for blast investigations and previously validated for a range of blasts, underwater explosions, and fluid structure interactions. Results: Further advancement of the Chinook software is detailed in this report with more advanced methods for material interface tracking and structural large deformation capability for fluid structure interaction (FSI). The local mesh adaptation technique modifies the CFD mesh in the cells adjacent to the structure, such that the cell faces correspond to the location of the structure, allowing for precise determination of pressure loads and interface velocities. The local mesh adaptation also allows for free surface tracking, such as bubble surfaces or wave movement, without the diffusion generated by mixed material models. The scheme is fully conservative. A series of validation cases were performed in order to verify the FSI and material interface tracking functionality and validate the results. These consisted of a moving plate piston problem, a shock tube, expansion of a casing containing high pressure gas, and shock propagation in a casing. UNDEX specific cases included are a bubble expansion, bubble merging, and cases with steel plate deformation due to underwater explosions. Results show a dramatic improvement in material interface diffusion compared to mixed cell formulations. Computational times are generally increased, but for some cases actually decreased by several orders of magnitude. Significance: The large deformation capability, initially detailed here, can be applied for any FSI interaction, such as applications in underwater explosions, particle break-up, aerodynamic ablation, and structural failure. The large deformation capability allows for assessment of failure properties and the behaviour of structures under high-rate impulse loading. Future plans: In order to make the Chinook FSI and material tracking simulations a truly robust tool a few issues still remain to be resolved, including dealing with stiff equations of state, interface smoothing, and the addition of surface tension terms. Further validation test cases for larger scale UNDEX test problems need to be performed. Results should be compared to further experimental data to confirm the simulated results are appropriate.

DRDC Atlantic CR 2011-343

iii

Sommaire ..... Large Deformation Fluid-Structure Interaction and Interface Tracking for Chinook Derrick Alexander; DRDC Atlantic CR 2011-343; R & D pour la défense Canada – Atlantique; Février 2012 Introduction ou contexte : La capacité à évaluer le chargement, la réponse et la survivabilité de structures et d’équipement navals aux fins de résistance aux effets des armes sous-marines contenant des explosifs hétérogènes est d’une importance fondamentale, en particulier pour la détermination des dommages potentiels de la menace représentée par les nouvelles technologies d’armes et pour la mise au point de contre-mesures efficaces. Martec Limited utilise actuellement le logiciel Chinook à dynamique numérique des fluides (DNF) pour simuler des explosions sousmarines. Ce logiciel a été spécialement conçu pour des enquêtes sur une explosion et a servi précédemment à valider une diversité de souffles, d’explosions sous-marines et d’interactions fluides-structures. Résultats : L’avancement du logiciel Chinook est expliqué en détail dans le présent rapport, avec des méthodes plus avancées de suivi de l’interface des matériaux et de déformation structurale importante pour l’interaction fluides-structures (IFS). La technique d’adaptation locale de maillage modifie le maillage de DNF dans les cellules adjacentes à la structure, de manière à ce que les faces de cellule correspondent à l’emplacement de la structure, ce qui permet une détermination précise des charges de pression et des vitesses d’interface. L’adaptation locale de maillage permet également un suivi de surface libre, comme les surfaces de bulles ou le mouvement des vagues, sans la diffusion générée par des modèles de matériaux mixtes. Cette mesure est entièrement conservatrice. Une série de validations ont été effectuées afin de vérifier les fonctions d’IFS et de suivi de l’interface des matériaux et afin de valider les résultats. Ces validations ont touché un problème de piston de plaque mobile, un tube à chocs, l’expansion d’un boîtier contenant du gaz haute pression et la propagation de chocs dans un boîtier. Les cas d’explosions sous-marines comprises sont une expansion de bulles, une convergence de bulles et des cas de déformation de plaque en acier en raison d’explosions sous-marines. Les résultats démontrent une amélioration spectaculaire de la diffusion de l’interface des matériaux en comparaison aux formulations de cellules mixtes. Les temps de calcul sont généralement augmentés, mais ils ont diminué d’un certain nombre d’ordres de grandeur pour certains cas. Importance : La capacité de déformation importante, qui était expliquée en détail dans le présent document à l’origine, peut être appliquée pour toute IFS, comme les applications lors d’explosions sous-marines, de rupture de particules, d’ablation aérodynamique et de défaillance structurale. Elle permet l’évaluation des propriétés de défaillance et du comportement des structures sous l’effet de chargement d’impulsion à grand débit. Perspectives : Il faut effectuer d’autres essais de validation pour des problèmes d’essais d’explosions sous-marines à plus grande échelle. Les résultats doivent être comparés à d’autres données expérimentales afin de confirmer que les résultats simulés sont appropriés. iv

DRDC Atlantic CR 2011-343

Table of contents Abstract …….. ................................................................................................................................. i Résumé …..... ................................................................................................................................... i Executive summary ........................................................................................................................ iii Sommaire ....................................................................................................................................... iv Table of contents ............................................................................................................................. v List of figures ................................................................................................................................. vi List of tables ................................................................................................................................. viii 1

Introduction............................................................................................................................... 1 1.1 Tracking Methodologies ................................................................................................ 2 1.2 Overview of Chosen Methodology................................................................................ 2

2

Numerical Method .................................................................................................................... 4 2.1 FSI Initialization ............................................................................................................ 5 2.2 FSI Update ..................................................................................................................... 7 2.3 Altered Mesh ............................................................................................................... 15 2.4 Usage ........................................................................................................................... 25

3

Validation Cases ..................................................................................................................... 28 3.1 Moving Plate ............................................................................................................... 28 3.2 Shock Tube .................................................................................................................. 30 3.3 Casing Expansion ........................................................................................................ 31 3.4 Meso Scale Shock........................................................................................................ 34 3.5 Bubble Expansion........................................................................................................ 36 3.6 Bubble Merging ........................................................................................................... 38 3.7 Plate Deformation........................................................................................................ 40

4

Conclusions and Recommendations ....................................................................................... 44

References ..... ............................................................................................................................... 45 List of symbols/abbreviations/acronyms/initialisms ..................................................................... 49 Distribution list .............................................................................................................................. 51

DRDC Atlantic CR 2011-343

v

List of figures Figure 2-1: FSI coupling types ........................................................................................................ 4 Figure 2-2: Behaviour of altered CFD mesh (black) due to structure(in red) ................................. 5 Figure 2-3: FSI initialization flow chart .......................................................................................... 6 Figure 2-4: FSI update flow chart.................................................................................................... 8 Figure 2-5: Determination of intersected fluid faces ....................................................................... 9 Figure 2-6: HLLC Riemann fan with two intermediary states ...................................................... 11 Figure 2-7: Fluid mesh alteration to align with structure .............................................................. 15 Figure 2-8: Mesh alteration flow chart .......................................................................................... 16 Figure 2-9: Projection of point onto triangle ................................................................................. 17 Figure 2-10: Restricted boundary vertex movement ..................................................................... 19 Figure 2-11: Projection of point onto line segment ....................................................................... 19 Figure 2-12: Cell detected on opposing sides of structure ............................................................ 20 Figure 2-13: Cell detected on opposite side of curled structure .................................................... 21 Figure 2-14: Cell side determination for unstructured grids ......................................................... 21 Figure 2-15: Collapsed cell generation.......................................................................................... 22 Figure 2-16: Mesh movement and conservative passing for switching cells ................................ 23 Figure 2-17: Issues with passing conservative information through faces .................................... 23 Figure 3-1: Plate in duct schematic ............................................................................................... 28 Figure 3-2: Pressure contours for 2D triangular grid .................................................................... 29 Figure 3-3: Pressure distribution at 0.8 ms .................................................................................... 30 Figure 3-4: Pressure distribution at 0.55 ms .................................................................................. 31 Figure 3-5: Casing expansion initial regions ................................................................................. 32 Figure 3-6: Pressure contours at 0.3 ms ........................................................................................ 33 Figure 3-7: Casing mass fraction contours at 0.3 ms..................................................................... 33 Figure 3-8: Total mass and casing mass throughout simulation.................................................... 34 Figure 3-9: Density and pressure contours illustrating propagation of shock and rarefaction waves in time with material interface tracked ............................................................. 35 Figure 3-10: Density and pressure contours illustrating propagation of shock and rarefaction waves in time with the mixed cell formulation ........................................................... 35 Figure 3-11: UNDEX bubble shape and pressure history ............................................................. 36 Figure 3-12: UNDEX maximum bubble radius ............................................................................ 37 vi

DRDC Atlantic CR 2011-343

Figure 3-13: Initial high pressure region mass fraction contours and mesh distribution .............. 38 Figure 3-14: Merging of high pressure regions (red) .................................................................... 39 Figure 3-15: Single high pressure region with material interface tracking ................................... 39 Figure 3-16: Single high pressure region without material interface tracking .............................. 39 Figure 3-17: Charge under steel plate schematic........................................................................... 40 Figure 3-18: Explosive mass fractions (colour), structural plate (black line), and explosive material interface (white line) at 0.1398 for 0.001 m stand-off ................................. 41 Figure 3-19: Structural plate vertical velocity contours at 0.1398 for 0.001 m stand-off ............. 42 Figure 3-20: Pressure (colour) and material interface contours (black lines) for 0.02 m standoff ................................................................................................................................ 42

DRDC Atlantic CR 2011-343

vii

List of tables Table 3-1: Casing expansion initial conditions ............................................................................. 32 Table 3-2: UNDEX bubble simulation comparisons ..................................................................... 38 Table 3-3: Maximum plate displacements..................................................................................... 43

viii

DRDC Atlantic CR 2011-343

1

Introduction

The ability to assess the loading, response and survivability of naval structures and equipment to withstand the effects of underwater weapons containing non-ideal explosives is of fundamental importance; particularly for determining the threat damage potential of new weapon technologies, and developing effective counter-measures. Since 1993 Martec Limited has been developing blast loading and response programs for analysis of naval structures and targets subjected to underwater explosions (UNDEX). Continued development has examined the effects of: hull shell structures with both air and water backed conditions, structural buoyancy, hull whipping response, and for non-UNDEX problems modeling the effects of metalized explosives and warhead casings. Currently Martec utilizes the computational fluid dynamics (CFD) code Chinook to simulate UNDEX events. The Chinook platform is ideal for analyzing sub-surface, near-surface, and above-surface blast threats to Canadian Forces vessels. Chinook has been specifically designed for blast investigations and previously validated for a range of blasts, underwater explosions, and fluid-structure interactions [1]-[5]. Martec has demonstrated the ability to model multidimensional, multi-phase explosive events. Chinook includes the capability to model non-ideal equations of state that can be critical to model key experiments. Further advancement of the Chinook software is detailed in this report with more advanced methods for material interface tracking and structural large deformation capability for fluid structure interaction (FSI). Chinook can be coupled with the finite element (FE) codes LS-Dyna [6] or Abaqus [7] in order to simulate fluid structure interaction. FSI is required for investigations where the structural and fluid responses are strongly coupled. When underwater explosions occur in close-proximity to a structure, large deformation coupling, where the structural movement is accounted for in the fluid domain, is required since the structural response significantly influences subsequent fluid flow. An immersed FSI approach is used in which the structure and fluid are meshed independently, and the interaction location is dynamically determined. Appropriate boundary conditions are applied to the fluid boundary as the structural mesh moves across the fluid cells. Previous the location of information passed between the solvers was the closest fluid cell face to the structural surface. When the structural surface moved a “cell erosion” technique was utilized in which complete fluid cells were uncovered or covered up by the movement. Comparisons of prior UNDEX simulations to experiments indicated that the numerical simplifications made to model FSI were not necessarily appropriate. When a structure moved, the material on either side of the structure did not respond appropriately. When the structure moved a large distance compared to the mesh size, discrepancies arose because the material interface, separating the fluid on either side of the structure, did not follow the deforming surface of the structure. In addition to the FSI implementation, there was previously no way to track material interfaces. For UNDEX problems this is required to avoid numerical diffusion on the explosive bubble interface and the free surface interface. Without the bubble interface tracking, the diffusion of material properties will degrade the results for the bubble movement beyond the initial expansion.

DRDC Atlantic CR 2011-343

1

1.1

Tracking Methodologies

There exist a number of interface tracking methodologies designed to improve numerical simulations. Reviews of the competing methodologies are given in [8]-[10]. For implementation in Chinook it was desired to find a methodology that worked for both FSI with thin and thick shell structures as well as for material interface tracking. Tracking methodologies fall broadly into two categories; front tracking where the movement of the interface is advected directly or interface reconstruction where materials are advected to mixed cells and then the location of the interface is reconstructed. Numerous parties have worked on these methods and there is blurring of the methods such that distinct categorization is not always possible. Front Tracking Methods Lagrangian methods [11]-[13] distort the mesh such that the fluid in a cell never leaves the cell. Large distortions can induce errors and often tangled meshes will result that require remeshing and mapping of the solution onto the new grid. Arbitrary Lagrangian-Eulerian (ALE) methods [14]-[16] are Eulerian away from the interface and have a moving Lagrangian mesh near the interface. This removes a number of the deficiencies of Lagrangian schemes; however their implementation can be complex in multiple dimensions or with multiple interface intersections. Level set methods [17]-[19] capture the location of the interface by solving a governing equation for the advection of a characteristic material variable. Front tracking methods are similar to level set methods except advection of an interface mesh is done by a Lagrangian process. Velocities at the interface are based on the solution of a Riemann problem which exactly captures interface discontinuity. For both these methods where the interface is tracked separately from the fluid cells, the issue remains how to appropriately treat the cells which the interface intersects. Interface Reconstruction Methods In a Volume of Fluid (VoF) method [20]-[22] the volume of each material is tracked. Cells that contain more than one material contain a portion of the interface. At each iteration, the material fractions in these cells are used to reconstruct an approximation to the interface and this approximate interface is then used to update values of the fluid volumes in each cell at the next iteration.

1.2

Overview of Chosen Methodology

Since the previously implemented FSI routines provide structural interface advection information directly, front tracking was chosen for the material interface tracking as well. While level set methods are popular for determination of interface locations in compressible flow fields, they do not necessarily handle non-enclosed regions like an immersed plate. A commonly used approach for interface tracking is to simply have the structure “cut” the fluid cells into multiple smaller cells [23][24]. This method has several drawbacks including being non-trivial for unstructured three-dimensional meshes, generating complex cell and face shapes, and potentially generating vanishingly small cells. The small cells limit the global time step 2

DRDC Atlantic CR 2011-343

through Courant–Friedrichs–Lewy (CFL) conditions. To avoid these restrictions, merging criteria must be used to arbitrarily join the small cells to neighbouring cells. A more accurate locally adapting mesh technique has been developed to model large FSI deformations in structures and capture material interfaces. Similar to the ideas behind ALE methods, the grid in the immediate vicinity of the interface can instead be altered such that the cell faces lie upon the interface itself. In the vicinity of the fluid-structure interface, the fluid mesh vertices will be altered to lie exactly on the structural elements. The structural velocity is therefore correct for the mesh region on the interface. The work of Sachdev [25] applied the local mesh adaptation routine to two-dimensional, structured, Cartesian grids. This idea has been extended to work with three-dimensional, unstructured grids with arbitrary cells types. The FSI technique developed is applicable to a wide variety of problems including underwater explosions with motion-induced cavitation and structural failure, surface ablation, and even combustion of solids. The local mesh adaptation also allows for free surface tracking, such as bubble surfaces or wave movement, without the diffusion generated by mixed material models. .

DRDC Atlantic CR 2011-343

3

2

Numerical Method

Chinook is a three-dimensional, unstructured, finite volume code which reduces cell count and allows for boundary-fitted meshes. When simulating a fluid structure interaction (FSI) problem, Chinook can be dynamically coupled with a commercial finite element solver such as LS-Dyna [6] or Abaqus [7]. The CFD fluid solver is used to determine the fluid load on a structure and the finite element solver is used to determine the response of the structure. An immersed approach is utilized in which the structural and the fluid domains are meshed independently in each software suite and are not required to align. Coupling between the fluid and structure is accomplished by dynamically determining the intersection of the structure and fluid meshes (the wetted surface) as the simulation progresses, and using the intersected surface to calculate structural loads and apply fluid boundary conditions. Depending upon the structural properties and the expected range of motion different types of FSI coupling can be utilized, as noted in Figure 2-1. The first option is one-way coupling of the pressure load from Chinook to the finite element (FE) solver. This option should be used when loading from the fluid domain is desired but the structure has little influence on the fluid domain, such as encountered with a rigid body. The second option couples the Chinook pressure load to the FE solver and couples the structural velocities to Chinook. Structural velocities are applied as boundary conditions at the non-moving wetted surface in the fluid domain. This option should be used for small deformation scenarios where the structural displacement has little influence on the fluid domain but the velocity of the structure influences the fluid, such as cavitation with small structural movements. The third coupling option transfers the Chinook pressure load to the FE simulation and transfers the structural node positions, node velocities and element failure states from the FE simulation to the fluid simulation. These updated positions, velocities, and failure states are used to calculate an updated intersection surface and to apply velocity boundary conditions on the fluid. This option should be used for large deformation problems, where the structure has a large influence on the fluid. Fluid Explosive Structure Original Position

Two-way Coupling Pressure and Velocity

Two-way Coupling Pressure, Velocity and Displacement

Figure 2-1: FSI coupling types In previous small-deformation analyses [1][2], velocity coupling was sufficient to model the cavitation in the water caused by the structural response, because of the relatively low magnitudes of deformation. In contrast, large-deformation simulations require the fluid to move with the structure as it deforms, particularly for multi-material calculations. This report details an extension of the large deformation FSI modeling capability in Chinook. In prior work [3][4] the 4

DRDC Atlantic CR 2011-343

grid location of information passed between the solvers was the closest fluid cell face to the structural surface. When the structural surface moved a “cell erosion” technique was utilized in which complete fluid cells were uncovered or covered up by the movement. The current technique improves on this calculation by modifying the CFD mesh in the cells adjacent to the structure, such that the cell faces correspond to the location of the structure. The local mesh adaptation also allows for free surface tracking, such as bubble surfaces or wave movement, without the diffusion generated by mixed material models. The CFD and FE solvers each have a numerical mesh. Throughout this document the CFD mesh is referred to as having cells, faces, and vertices, while the structural mesh imported into Chinook is referred to as having elements and nodes. In addition for material interface tracking a mesh is generated to represent the interface. This is referred to as the structure (similar to the FE mesh) and is treated in the same manner except when explicitly stated otherwise. Figure 2-2a shows a CFD mesh in black and a structural cylinder in red. The CFD mesh is altered to capture the cylindrical shape, Figure 2-2b. As the cylinder propagates from right to left in time, the altered CFD cells return to their original shape once the structure is no longer in their vicinity, Figure 2-2c.

a) Original CFD mesh

b) Initial altered CFD mesh

c) Subsequent CFD mesh

Figure 2-2: Behaviour of altered CFD mesh (black) due to structure(in red) The FSI routines in Chinook are largely independent of the main flow solver. Details of the FSI initialization, FSI update functions, and altered mesh routines used in Chinook are given in Sections 2.1 to 2.3. Details of the Chinook keyword usage for running a FSI and/or material interface tracking case are given in Section 2.4.

2.1

FSI Initialization

A flow chart of the FSI initialization routines is given in Figure 2-3.

DRDC Atlantic CR 2011-343

5

FE solver started: grid and initial conditions determined Chinook started: grid and initial conditions determined

Wait for FSI initialization

Start FSI initialization Generate cell space tree based on center and extent Flag boundary cell movement restraints

Material interface locations from different state conditions

Exchange initial FSI data with FE solver (elements nodes, nodal positions and velocities)

Interface velocities from Riemann problem at corresponding faces

Generate interface elements and nodes from faces Element areas, normal vectors Generate element neighbours Perform FSI update Figure 2-3: FSI initialization flow chart An FSI calculation is performed by first starting a version of the FE solver (LS-Dyna or Abaqus) that has been modified to accept FSI communications. Chinook and the FE solver communicate information at user defined time intervals through socket connections. After start up the FE solver will wait for communication from Chinook. Currently the FE solver is run in serial and only communicates to one Chinook process. For a parallel CFD calculation the information about the structural mesh is passed to the other Chinook processes within Chinook itself. 6

DRDC Atlantic CR 2011-343

As detailed in Section 2.3 the vertices that lie on the boundary of the fluid domain are not allowed to move off of the boundary. To facilitate this restriction all boundary cells are initially checked and their vertices are flagged to have no displacement normal to the boundary surface. In order to maintain fluid domain features vertices on the junction between two different boundary conditions are flagged to remain restricted to the line separating the two boundaries. This retains the initial geometric features of the Chinook domain, provided sufficient initial boundary conditions are defined. For material interface tracking the initial location of the interface is defined by searching through all faces and checking for those that bound a user specified state condition (volumetric initial condition region). The velocity at which these faces propagate is determined from the solution of a Riemann problem and a “structural” mesh is generated from these faces, which is detailed further in Section 2.2. The structural elements, node positions, and node velocities are initially passed to Chinook, to allow for a representation of the structural surface to be assembled in Chinook. While Chinook can perform 2D and 3D calculations the structural model must be a 3D shell. Solid structural components can be used by applying a shell region around the solid [4]. A data structure is built containing all of the structural mesh elements and nodes. If the structure is not in the same coordinate frame as the fluid, its location may be translated and rotated to align correctly. For multiple FSI components or a combination of material interface and FSI components each element and node is marked by the “region” it belongs to. This allows for independent tracking of multiple structures, such as required for a bubble interface and FE plate. All structural elements are assumed to be triangles. If the FE solver passes quadrilaterals or the material interface has quadrilateral faces, they are divided into two triangles along the node 1-3 line segment. For 2D material interfaces a third node is added at the center of the two element nodes, offset in the third dimension by the half the distance between the first two nodes. This greatly simplifies the routines that deal with the structural element information as they do not have to contain contingencies for multiple shapes. Structural element neighbours are determined by storing a table of elements for each node. Neighbour elements are those which share a common node. The FSI neighbours are only determined at the start of the calculation and are assumed not to the change throughout the calculation. Material interface neighbours must be determined whenever the interface is regenerated.

2.2

FSI Update

For each time iteration in Chinook, the cell fluxes and source terms are determined, the conservative values for the next time step are determined, and the primitive variable values are updated. A FSI update occurs subsequent to these other steps during each iteration, as shown in the flow chart of Figure 2-4. Further information on each step in the flow chart is given below.

DRDC Atlantic CR 2011-343

7

Chinook iteration performed

FSI update

Update nodal positions

FE iteration performed

Remove prior intersections

Coupling time?

Determine new intersections Small deformation?

Yes: Wait for FSI exchange Yes

No

FSI coupling time?

Alter Mesh

Yes

Material update freq? No Determine velocity of faces on interface Average velocity on material element

No

Yes Material interface locations from faces on interface Interface velocities from Riemann problem at corresponding faces Generate material interface element neighbours

Material node velocities

Determine pressure force across faces on interface Maximum pressure on FSI element Force diffusion Exchange FSI data (nodal forces, new velocities) FSI face velocity from elements

Continue to next iteration Figure 2-4: FSI update flow chart

8

DRDC Atlantic CR 2011-343

Updated Elements The positions of the structural nodes are updated at each iteration from the prior iteration nodal velocities.

xi = xi , prior + vi ∆t

(1)

The element areas and normal directions are then recalculated. Remove Prior Intersections References to any faces or cells that were previously flagged as intersected or having their geometry altered are removed. Face velocities are set to zero (non-moving). Intersection Detection Subsequent to the update of the structural element positions, the location where the structure intersects the fluid cells (the wetted surface) is determined. The unstructured meshes utilized in Chinook complicate the calculation of the intersections. An intersection occurs where there is a fluid face having adjacent fluid cells on opposite sides of a two-dimensional structure, i.e. that has a dual segment that intersects the structure, Figure 2-5b. With structured meshes, determining the intersections is trivial since the nearest fluid cell/face to a given structural element may be determined directly based on the fluid mesh extents and spacing. Unstructured grids require a search of the fluid mesh potentially involving millions of fluid faces and structures with tens of thousands of elements. Such a search would be intractable if performed directly, since the cell/element counts on meshes of such sizes could exceed ten billion intersection tests per time step if a brute-force approach was used. Fluid grid Structure Fluid cell duals a)

b) Figure 2-5: Determination of intersected fluid faces

Rather than perform a brute force search, a multi-step process is used to determine intersections. First the fluid mesh is added to a space tree based on the cell center and extent. The tree branches are divided by the average of the cell coordinates (in the largest coordinate extent direction). At each time step, every structural element searches through the space tree to find a set of candidate fluid cells that are within the element’s bounding box. For each candidate cell the dual line segments to the face neighbour cell centers are then checked for intersection with the element geometry using the segment-triangle test [26]. If an intersection exists, the face between the two DRDC Atlantic CR 2011-343

9

cells sharing the intersected dual segment forms part of the wetted surface and is said to be intersected by the structure. Each intersected face is flagged and the structural element stores a link to the face. With the location of fluid solver interaction determined as the bold faces shown in Figure 2-5b the geometric mismatch between the fluid and structural boundaries can be seen. In the previous FSI implementation in Chinook this created imprecise application of the fluid loads and structure velocities. More importantly the changes in force and velocity were not correct when the structure moved relative to the fluid since volumetric changes would only occur when the structure moved far enough to be detected at a different fluid face. In order to produce the correct fluid motion the actual wetted surface location, as specified by structure, should be tracked by the fluid. FSI Face Velocity If the mesh is not being altered, as is the case when small deformation coupling is assumed and no material interfaces are present, the velocity at the fluid face, vn,face, (in the face normal direction) can be determined as the average of the intersected element values. The FSI element velocities are the average of the nodal values.

vi ,element = v n , face =

1 ∑ vi,node 3 3

cv   velement ⋅ n face ∑ # intersections # intersections

(2)

Where cv is a user defined scaling coefficient (1 for full velocity applied) and nface is the face normal vector. The nodal velocities to be applied at the next Chinook iteration are determined by one of four methods. Regenerated Material Interface Nodal Velocity For material interfaces, the velocity of the fluid face can be determined via a solution of the Riemann problem at the intersected face. Using the notation of [27] for the HLLC approximate Riemann solver (which is utilized in the solution of the advection fluxes in Chinook), the velocity of the intermediary states vL* = vR* = SM, Figure 2-6, can be determined via the left and right states. The calculation of the interface velocity assumes that the face is non-moving and hence, vface = 0.

10

DRDC Atlantic CR 2011-343

t

SL

vface

L* L

SM

SR

R* R x

Figure 2-6: HLLC Riemann fan with two intermediary states

SM =

ρ R v R (SR − v R ) − ρ L v L (SL − v L ) + p L − p R ρ R (SR − v R ) − ρ L (SL − v L )

(3)

Where ρ is density, v is velocity (normal to face), and p is pressure. The left and right wave speeds are approximated via the maximum and minimum wave speeds.

SL = min (v R - a R , v L - a L ) , SR = max(v R + a R , v L + a L )

(4)

Use of a Roe averaged state [28] for the wave speed determination was not found to provide any appreciable benefit. For different equations of state on either side of the interface, the HLLC interface velocity calculation, given by equation (3), may not be accurate [29][19]. The velocity of the center state can be approximated by a linearization of the wave curves for each state (Young-Laplace Law at planar phase boundaries). This uses the fact that the velocity of the L* and R* states are the same and the pressures are the same except for a surface tension term.

v L* = v L −

aL

ρL

ε L = vR +

aR

ρR

ε R = v R*

p L = p L + a ε L = p R + a ε R − (d − 1)σκ *

2 L

(5)

2 R

Where d is problem dimension, σ is surface tension, and κ is the mean curvature of the interface. These equations can be solved for εL and εR,

DRDC Atlantic CR 2011-343

11

a 2 ε L  1  R ε  = a  a 2  R  a L a R  a R + L ρ R   L  ρL

− aR aL



vL − vR    p − p − (d − 1)σκ  R   L 

ρ R 

ρL

(6)

The interface velocity, vL* = vR*, is then determined from equation (5). If the Chinook iteration corresponds to the frequency at which the material interface is regenerated then an element is created for each face that intersects the interface. The interface outward normal vector and velocity are added to data arrays. Each vertex on the face is added to a hash table which has a unique key based on the global vertex number. Local vertex numbers are sequential on each parallel process, while global numbers are unique for each vertex across all processes. This is done to avoid duplicating vertices from neighbouring faces. The arrays of face and vertex information are sent to all other parallel processes so that each process has a complete description of the entire interface. For each intersected face an element is generated with the appropriate nodes and velocity. Structural element neighbours are determined by having each node record which elements use it. The neighbour elements are then all those which share a node. Nodal velocities are determined from the area weighted average of the velocities from the elements that share the node (in this case velement = vface).

v node =

∑v ∑A

elements

element

elements

Aelement (7)

element

Non-Regenerated Material Interface Nodal Velocity If the material interface is not being regenerated at the current iteration then the velocity of the interface elements is determined from an area weighted average of the fluid face velocities that intersect the element.

velement =

∑v A ∑A

face int er sec tions

face

(8)

face int er sec tions

If the structure is not being regenerated there may be elements that do not intersect a fluid face (this also occurs for FSI structures). These may sit fully within one cell, with no obvious front and back states that can be used to determine the element velocity. To solve this problem, we use a diffusion operation that conceptually is equivalent to a linear interpolation between known neighbour values. The velocity component in each dimension for the non-intersecting elements is simply an average of the neighbour values.

12

DRDC Atlantic CR 2011-343

vi ,element =

1 ∑ vi,element # neighbours neighbours

(9)

The diffusion process is repeated until the change in the velocity between iterations has fallen below a tolerance for all elements. The tolerance is determined as the maximum change in velocity over all the dimensions:

 vi − vi , prior  tol = max tol , vi  

    i =1, 2 , 3 

(10)

Nodal velocities are determined from the area weighted average of the velocities from the elements that share the node as given by equation (7). Updated FSI Nodal Velocity For FSI elements, if the solution time corresponds to the coupling frequency with the FE solver then information is exchanged with the FE solver. It should be noted that typically the FE solution proceeds much faster than the CFD solution and the FE solver will usually be waiting for Chinook, otherwise Chinook will wait for the FE solver. The force at each node is determined from the fluid pressure as:

Fnode =

1 c force ∆pA elements 3

∑

(11)

Where cforce is a user defined scaling coefficient (1 for full force applied), ∆p is the pressure difference from the fluid cells on either side of the element, and A is the area of the element. The 1/3 factor is to distribute the force to each node of the triangular element. If a structural element intersects multiple fluid cell dual lines, the maximum pressure difference is used. If a structural element does not intersect any fluid cell dual lines, the pressure difference is determined by a diffusion operation from those that do. For 3D solutions, or 2D solutions with structural elements in the 2D plane, the pressure difference for the non-intersecting elements is simply an average of the neighbour values. If an element has failed it is not a valid neighbour and is not used in this summation.

∆p =

1 ∑ ∆p # neighbours neighbours

(12)

For 2D planar solutions the pressure difference for an element not in the 2D plane is determined from the neighbour element that is closest to the 2D plane, i.e. has the smallest out-of-plane z DRDC Atlantic CR 2011-343

13

coordinate magnitude. The diffusion process is repeated until the change in the maximum pressure difference between iterations has fallen below a tolerance for all elements. The tolerance is determined as:

 ∆p − ∆p prior tol = max tol ,  ∆p 

   

(13)

While the number of loops required for the pressure diffusion is problem dependent, four or five loops were generally required in the test cases examined in this report. For 2D axisymmetric solutions the pressure differences in the 2D plane are sorted by radial distance. For unknown elements the value is determined via linear interpolation of the sorted data. This is only valid for simple geometries that have a constant variation with radius, such as a plate. If a more complex geometry, such as a cylinder, is required the simulation must be performed in three dimensions. Updated nodal velocity values and element failure flags are received from the FE solver. Non-Updated FSI Nodal Velocity For FSI elements, if the solution time does not correspond to the coupling frequency the previous nodal velocity values are retained. Flux at Interface Faces During the next Chinook iteration, any faces marked as on the interface do not use the normal approximate Riemann solver to directly calculate flux from both neighbour cells. Instead two ghost cells are created for each face and a flux is determined separately for each interior cell. A slip wall condition is applied that uses the same fluid properties as the interior cell (mass fractions, pressure, temperature), but the velocity in the ghost cell is adjusted to generate the specified interface velocity.

   v ghost = 2vint erface − vint erior

(14)

The momentum and energy of the boundary cells are thus altered.

QρE = h − p + 1 ρv ghost 2

(15)

The flux at the interface is determined via a moving mesh formualtion where the flux through a moving face is determined via equation (16), that subtracts the amount of the conservative quantity, Q, that would move to the neighbouring cell due to the movement of the face.

14

DRDC Atlantic CR 2011-343

∂ [Q]δV + ∫ ([F ] − [Q]vint erface )⋅ δA = ∫ [S ]δV ∂t V∫ A V

(16)

Where t is time, V is cell volume, F is conservative flux, A is face area, and S is a source term. For a specified interface velocity, vinterface, this formulation generates zero mass flux through the interface. There can be a change in momentum and energy as the structure alters the velocity of the surrounding fluid. There is therefore no mechanism for material on one side of the interface to travel to the other side, as the cell on each side has its flux determined independently.

2.3

Altered Mesh

For material interface tracking or large deformation FSI the fluid mesh is altered to align with the structural mesh. The concept is shown schematically in Figure 2-7, where the vertices of intersected fluid faces are projected to the closest point on the structure. The intersected faces are then aligned with the structure allowing for a more appropriate determination of force at the structure and velocity. If a cell has been previously altered but is no longer within a distance of one cell from the structure, it returns to its original mesh coordinates. A flow chart of the mesh alteration routines is shown in Figure 2-8. Fluid grid Structure Fluid cell duals a)

b)

c)

Figure 2-7: Fluid mesh alteration to align with structure

DRDC Atlantic CR 2011-343

15

For all vertices on intersected faces flag and project vertex to structure Determine interface side for all cells with altered/prior altered vertex Check for collapsed cells Recalculate altered cell and face geometries Alter cell conservative quantities Determine face velocities Figure 2-8: Mesh alteration flow chart Vertex Projection Each vertex on an intersected face is flagged as lying on the interface. The vertex is then projected to the closest position on the structure (normal projection). When using FSI and material interfaces in the same simulation, the FSI face is always given precedence. This mean that if both a FSI structure and a material interface structure intersect a face, its vertex will be moved to the closest FSI structure. The movement of the vertices is restricted if they lie on the CFD domain boundary. If the vertices are interior to the domain their movement is unrestricted and are projected normal to the surface (to lie on the closest structure). If the vertex lies on a single boundary condition, its movement is restricted to lie on the same boundary surface. Vertices that share two boundary conditions are restricted to lie on the line joining those boundary conditions in three-dimensional problems or are restricted to not move in two-dimensional problems. Vertices that share three boundary conditions are restricted to not move in three-dimensional problems. Projection of a vertex onto the structure is accomplished via the following four steps. First, the structural elements are recursively searched, starting from the one the face intersects and then through its neighbours, until the element with the closest center location to the vertex is found. Second, the vertex location is projected onto the element. The projected location must be within the element. Projection of point p (the vertex) onto a triangle (the element) with nodes a, b, and c, Figure 2-9, can be accomplished via the following [30].

16

DRDC Atlantic CR 2011-343

p a

c

pstructure

b Figure 2-9: Projection of a point onto a triangle Taking the following dot products:

d1 = ab ⋅ ap , d 2 = ac ⋅ ap , d 3 = ab ⋅ bp d 4 = ac ⋅ bp , d 5 = ab ⋅ cp , d 6 = ac ⋅ cp

(17)

If d1 and d2 are zero or less, the projected position of p is the position of vertex a. If d3 is zero or greater and d4 is zero or less, the projected position of p is the position of vertex b. If d6 is zero or greater and d5 equal to less than d6, the projected position of p is the position of vertex c. The following quantities are used to determine if the projected position of p is closest to a triangle edge:

va = d 3 d 6 − d 5 d 4

, vb = d 5 d 2 − d 1 d 6

, vc = d1 d 4 − d 3 d 2

(18)

If vc is zero or less, d1 is zero or greater, and d3 is zero or less, the projected position is given by:

 d1  d1  + b p structure = a1 − d1 − d 3  d1 − d 3 

(19)

If vb is zero or less, d2 is zero or greater, and d6 is zero or less, the projected position is given by:

 d2 p structure = a1 −  d2 − d6

 d2  + c d2 − d6 

(20)

If va is zero or less, (d4-d3) is zero or greater, and (d5-d6) is zero or greater, the projected position is given by:

DRDC Atlantic CR 2011-343

17

 d4 − d3 p structure = b1 −  d4 − d3 + d5 − d6

 d4 − d3  + c d4 − d3 + d5 − d6 

(21)

If none of the prior conditions occur, the projected point is inside the triangle and is position is given by:

p structure =

av a + bvb + cvc v a + vb + v c

(22)

Third, the distance between the original vertex location and the projected location is recorded. To ensure the closest position to the structure has been determined the vertex is also projected onto the element neighbours. If a neighbour element has a shorter distance between the vertex and projected location then its neighbours are checked in a recursive process. Fourth, if the projected distance found from the fluid face is the smallest of all the values for the vertex then the vertex is moved to the projected location. For parallel calculations, vertices that are used for multiple processes exchange the flag denoting that they lie on the interface and the minimum distance to structural regions from all processes. The closest location amongst all processes is used. For a restricted boundary vertex the closest location on the boundary to the structure is determined by first finding the position projected normally onto the structure and then projecting this position back onto the boundary. In this case the minimum distance used for determining where the closest projected point is located is the distance between the point on the structure and the projected position on the boundary. This method of boundary restriction will not always generate the vertex in the correct spot on boundary. See for example the case illustrated in Figure 2-10a, where the normal projection of the vertex onto the structure and back to the boundary does not correspond with the actual intersection of the structure and the boundary. This mismatch does not arise in most cases, but to reduce this error a loop that successively projects to the structure and back is added to the boundary, with a maximum of 10 projections, as shown in Figure 2-10b.

18

DRDC Atlantic CR 2011-343

poriginal

pboundary

poriginal

pstructure pboundary

pdesired

a) Mismatch of initial guess and desired location

b) Successive projections

Figure 2-10: Restricted boundary vertex movement For a boundary vertex on two boundary conditions the closest location on the boundary to the structure is determined by first finding the position projected normally onto the structure and then projecting this position back onto the line segment separating the boundary conditions, Figure 2-11. In this case the minimum distance used for determining where the closest projected point is located is the distance between the point on the structure and the projected position on the line segment. The projected position on the line segment is found as follows. pline

b

a

pstructure

Figure 2-11: Projection of point onto line segment Taking the following dot products:

d1 = ap ⋅ ab , d 2 = ab ⋅ ab

(23)

If d1 is zero or less, the point is the position of vertex a, if d1 is equal or greater than d2 the point is the position of vertex b, otherwise the point on the line is given by: DRDC Atlantic CR 2011-343

19

pline = a +

d1 ab d2

(24)

Interface Side During the determination of fluid-structure intersections, as described in Section 2.2, the cells on either side of the intersected faces are marked as being on the inside or outside of the interface, depending on a dot product between the cell center dual line segment vector and the element normal direction. This side is recorded for all the different structural regions for which the cell has an intersecting face. In some instances a single cell can be found to lie on opposing sides of the same structural region. In Figure 2-12 the dual line between cell A and B will intersect structural element 1, 2, and 3. Cell A will be determined to be on the inside of element 1 and 3, but on the outside of element 2. If a cell is determined to lie on an opposite side of the structure from the side given by a previous element the fluid face is removed as an intersection. This is shown in Figure 2-12b where the face separating cell A and B would not be taken as an intersection. This allows for a smooth interface to be retained and provides a mechanism for merging/separating interfaces. If the fluid face is aligned with more structural elements in one direction than the other, the structure is assumed to lie in this direction. In Figure 2-12a for example, the dual line between cells A and B intersects structural elements 1, 2, and 3, but the normal direction at the face is assumed to be in the element 3 direction. B

3

B

2 2

1

1

A A

a) 3 elements intersect cell dual line

b) 2 opposite elements intersect cell dual line

Figure 2-12: Cell detected on opposing sides of structure In some cases the side of the fluid cells can still be determined to lie on opposing side of the interface. In Figure 2-13 for example, cell A is on the inside of element 1, but the outside of element 3. The mismatch will not be removed since the side determination occurs through two separate faces. In these instances the side of the fluid cell is determined from the sides of its neighbour cells that do not lie through an intersected face. 20

DRDC Atlantic CR 2011-343

1

A

2 3

B

Figure 2-13: Cell detected on opposite side of curled structure For cells that have an altered vertex, but do not have an intersected face, a loop is performed through their neighbour cells that do not lie through an intersected face, until a side can be found. This requires a loop, through all the cells that share a vertex in order to ensure the side is determined. In Figure 2-14 cells A and F are known to be on the outside of the structure since they have intersected faces. The side of cell D may not be initially determined as the sides of cells C and E are possible unknown until all the cells have been examined at least once. A second loop through the cells would then pass the outside value from cell E to D.

B A

C D

E F

Figure 2-14: Cell side determination for unstructured grids In parallel the neighbour may lie in a neighbouring process, so the sides of halo cells (neighbour cells on other processes) must be exchanged. A similar process is performed for all cells that share a non-interface vertex that previously lay on an interface. If a cell is on the opposite side of the interface compared to the previous iteration it is marked as having switched sides. Collapsed Cells Due to the fact that all vertices can be moved independent of one another, the situation can arise where a cell has approximately a zero volume. In Figure 2-15 all the vertices of the center cell are modified, resulting in a very small cell, with close to zero volume. To avoid the severe time step restrictions associated with such small cells, they are considered “collapsed” when all vertices have been modified. Collapsed cells are marked as inactive and not used in the CFD calculation. In Figure 2-15 this means the fluid solution sees the structure as being slightly thicker in the location of the center cell.

DRDC Atlantic CR 2011-343

21

Collapsed cells that have an intersected face still require correct conservative information in order to determine the pressure difference on FSI faces and the Riemann velocity for material interfaces. These are determined from neighbour cells that do not lie through an intersected face.

a) Structural location (red) through original grid

b) Intersected faces in bold

c) Altered fluid mesh

Figure 2-15: Collapsed cell generation If two interfaces are merging together a cell may have all vertices altered, but have a volume that is near its original volume. Currently these cells will be marked as inactive, which will induce an error in the surrounding cells which receive the conservative quantities. Cell and Face Geometry Recalculation The fluid cell volume, centroid, and extent as well as the fluid face area, centroid, and normal direction are recalculated for all cells and faces that have altered vertices. Fluid face normals are generated from vectors created from the face vertices. Chinook uses the convention that the fluid face normal points from the left cell to the right cell. For collapsed cells the face normal direction may not point from the left to the right cell if the vertices are in an arbitrary order. When first reading the computational mesh in Chinook the dot product is taken between the face normal direction determined from the vertex vectors and the vector connecting the cells centers. If the dot product is negative the ordering of the face vertices is reversed. It should be noted that a similar procedure cannot be used during the calculation with collapsed cells due to the possible misalignment between the normal vector and the vector connecting the cells centers. For axisymmetric cases the cell centroid location can be calculated with a radial weighting such that the centroid is the center of volume, not the geometric center of the vertices. In an attempt to improve the solution of an axisymmetric bubble expansion this calculation was added, but was found to provide no benefit. Altered cells are generally at least half their original size, since the structure will be detected on another face (different set of altered cells) once the structure passes over the original cell centroid location. Smaller cells can arise for complex interface shapes (near a single cell). Reduced cell sizes will result in smaller time steps (due to CFL restrictions) and overall a greater number of iterations per simulation compared to the non-altered mesh simulations.

22

DRDC Atlantic CR 2011-343

Determination of Conserved Quantities The changing cell shapes and volumes must be correctly accounted for, in addition to appropriate handling of cells that switch sides of structure, such that the fluid materials on either side remain separate and the conservative quantities are actually conserved. When a fluid cell switches sides of the structure its conservative quantities are passed to non-switching cells, Figure 2-16. In the simplest case this results in the conservative information from the switching cell and its nonswitching neighbour cell being averaged into the neighbour (which now has a volume approximately the same as the two cells at the previous time step). Structure location on original grid Iteration 1

Iteration 2

Iteration 3

Altered grid

Movement of conservative quantities for side switching cell Figure 2-16: Mesh movement and conservative passing for switching cells In an initial implementation of the logic, when a cell switched sides of the interface, an attempt was made to pass conservative data through the intersected faces. The drawback to this method is that a switching cell may not share a face with those cells on the other side of the interface. This is illustrated in Figure 2-17, where the structure (in red) traverses diagonally across cell C. The conservative information from cells D and E are passed through the structural face to cell F, but cell C does not share a face and therefore the conservative information would be lost. A B

C

A D

B C

D

E F

E

F

Figure 2-17: Issues with passing conservative information through faces

DRDC Atlantic CR 2011-343

23

Recursive search algorithms could be used if a cell is detected that has no valid neighbour to pass to, i.e. cell C passes to D and then to F. The issue is that for unstructured grids there may be 20-30 intermediary cells, not just one. Since the code is domain decomposed in parallel, some or all of the 20-30 neighbours may be on different processors and this greatly increases the search methodology and time required to pass information. Instead conservative information is first compiled at the vertices which lie on the interface, or were previously on the interface. Any vertices that are used in more than one parallel domain simply have their conservative information from each process added together. The conservative information at the vertices is then used to generate the conservative information at the adjoining cells. The number of vertices on the interface, or previously on the interface, are determined for each cell. The prior conservative information at the vertices is determined by the volume weighted average with the previous cell volumes.



Q

⋅ Vol



prior  (Q ⋅ Vol )vertex = ∑  prior # interface vertices  cells



 cell

(25)

In terms of density, for example, this provides the mass value at the vertex. In addition the total new volume (with the newly adjusted cell sizes) at the vertex is determined.

Vol new   Vol vertex = ∑   cells  # interface vertices  cell

(26)

The new conservative values at each cell can then be determined by the average of the vertex values.

Qnew =

(Q ⋅ Vol )vertex # interface vertices Vol vertex ⋅ (# interface vertices)

∑

(27)

Face Velocities For material interfaces the velocity at the face is determined via the Riemann problem as noted in Section 2.2. For FSI interfaces the velocity is determined as described for small deformation cases in Section 2.2. Monitor Point Adjustment If monitor points are being utilized and occur in an altered cell, all of the cells which share a vertex are checked to see if any have a cell center closer to the desired monitor point location. This check currently does not adjust the monitor point if the closer cell lies in a different partition.

24

DRDC Atlantic CR 2011-343

2.4

Usage

Chinook problems are defined via a keyword file. For an FSI problem a special version of LSDyna or Abaqus must be used that includes functionality to couple with Chinook. The FE solver must be started prior to starting Chinook. The interface tracking routines currently only function with Eulerian phases, not Lagrangian particles. The keywords for simulating an FSI problem coupled with a FE solver or material interface tracking are as follows. FSI Keyword Format *FSI ! deform freq [Int] [Real] ! handle [Real] [Real] ! translation [Real] [Real] ! rotation [Real] [Real] ! length scale [Real] ! smoothing loops [Int]

host [String]

port [Int]

[Real] [Real] [Real] velocity scale [Real]

load scale [Real]

The data values are defined as follows: deform

A flag to enable/disable large deformation. Activation of this flag causes the deformed structural positions to be updated throughout the simulation. [Int] 0 = large deformation disabled (small deformation) 1 = large deformation enabled

freq

Frequency at which data is coupled between Chinook and the FEA code. [Real]

host

The IP address of the computer running the FEA code. For loopback to the local machine use “127.0.0.1”. [String]

port

Port number for coupling data between computers. Chinook and the FE code communicate through ports. Identical port numbers are required in both the Chinook and FE input files. Typical values used are above 10 500. The port number that the user specifies in the input files may be in use by another program on the computer. If this occurs, change the port number in the input files. If desired, contact your system administrator for a more detailed list of available ports, as this may be unique for each computer. [Int]

DRDC Atlantic CR 2011-343

25

handle

The x, y, and z coordinates of the point about which the structural model is rotated and translated in the fluid domain. [Real]

translation

The distances in x, y, and z by which the structure is translated in the fluid domain. [Real]

rotation

The angles at which the structure is rotated about the x, y, and z axes. Note that the handle point specifies the centre of rotation. The order of rotation is rotation about x, then about y, and then about z. [Real]

length_scale

Factor used to scale the structure model to the units of the fluid domain. [Real]

velocity_scale Activates a percentage of the structural velocity to be coupled back the fluid. Expected values range between 0.0 – 1.0. 0.0 = structural velocity is not passed to the fluid domain 1.0 = 100% of the structural velocity is passed to the fluid domain This can also be used to scale the structure model velocity units to those of the fluid domain. [Real] load_scale

Activates a percentage of the fluid load to be coupled to the structure. Expected values range between 0.0 – 1.0. 0.0 = fluid load is not passed to the structural domain 1.0 = 100% of the fluid load is passed to the structural domain This can also be used to scale the structure model force units to those of the fluid domain. [Real]

smoothing_loops Number of loops for which the conservative values next to the interface are averaged with their neighbour values on the same side of the interface. 0 = no smoothing. [Int] FSI Example The example below activates coupling of Chinook with LS-Dyna. This example uses large deformation. Data is passed through port 10500 to the computer with the IP address 192.168.100.47. Coupling occurs at a frequency of 0.25 μs. The structure and fluid models have the same units (i.e. no scaling). The structure is positioned in the same coordinates as the fluid domain (i.e. no rotation or translation). *FSI ! deform 1 26

freq 2.5e-7

host "192.168.100.47"

port 10300 DRDC Atlantic CR 2011-343

! handle 0.0 0.0 ! translation 0.0 0.0 ! rotation 0.0 0.0 ! length scale 1.0 ! smoothing loops 0

0.0 0.0 0.0 velocity scale 1.0

load scale 1.0

Material Interface Keyword Format *MAT_INTERFACE ! n_states [Int] ! state_num [Int] !(repeat as necessary) ! freq smoothing loops [Int] [Int] The data values are defined as follows: n_states

Material interface location is initially determined by the edge of states. This defines the number of states on the interior of the interface. [Int]

state_num

State numbers of the regions on the interior of the interface [Int]

freq

The frequency at which the material interface is regenerated from fluid faces aligned to the interface. [Int]

smoothing_loops Number of loops for which the conservative values next to the interface are averaged with their neighbour values on the same side of the interface. 0 = no smoothing. [Int]

DRDC Atlantic CR 2011-343

27

3

Validation Cases

A series of validation cases were performed in order to verify the FSI and material interface tracking functionality and validate the results. These consisted of a moving plate piston problem, a shock tube, expansion of a casing containing high pressure gas, and shock propagation in a casing. UNDEX specific cases were a bubble expansion, bubble merging, and cases with steel plate deformation due to underwater explosions. All FSI problems utilized LS-Dyna for the FE solver.

3.1

Moving Plate

A simple fluid structure interaction test case is a structural plate moving with a fixed 200 m/s velocity through a duct filled with air, as shown in Figure 3-1. This creates a compression region ahead of the plate and expansion region behind the plate. The structural plate was modelled in LS-Dyna as a 0.1 m by 0.1 m plate meshed with 10 by 10 quadrilateral elements and a fixed nodal velocity of 200 m/s in the y dimension. The duct is 1 m long in the y dimension with a square 0.1 m by 0.1 m cross-section. The plate is initially located at y = 0.5 m. The air was initially at 101325 Pa and 300 K. Coupling between Chinook and LS-Dyna occurred at a frequency of 1×10-4 s.

Figure 3-1: Plate in duct schematic This simulation was performed several times with a CFD mesh utilizing 3D hexahedral cells, 3D tetrahedral cells, 2D quadrilateral, and 2D triangular cells. Cells had side lengths of 0.01 m in all simulations. For the 2D simulations the z dimension was not utilized in the CFD solution. Figure 3-2 shows the grid and pressure contours at 0 ms and 0.8 ms. The high pressure region can be seen ahead of the plate and low pressure region behind the plate. The altered fluid mesh perfectly captures the structural plate location (y = 0.66 m) and sharp interface between the two sides.

28

DRDC Atlantic CR 2011-343

Figure 3-2: Pressure contours for 2D triangular grid The pressure vs. distance along the duct is shown in Figure 3-3 for the 2D quadrilateral grid (structured), the 2D triangular grid (unstructured), and the analytical solution at a time of 0.8 ms. In the quadrilateral (structured) case the movement is aligned with the fluid cell faces and a perfect jump across the plate is achieved. In the triangular (unstructured) case the unaligned flow generates some scatter in pressure about the correct value on either side of the plate. The 3D cases provide similar results to the 2D cases for the structured and unstructured grids.

DRDC Atlantic CR 2011-343

29

Figure 3-3: Pressure distribution at 0.8 ms

3.2

Shock Tube

The previous moving plate case was repeated with a pressure of 1.012325×107 Pa on the y < 0.5 m side and the plate velocity not specified. Coupling between Chinook and LS-Dyna occurred at a frequency of 1×10-5 s (every 5-10 Chinook iterations). The pressure difference across the plate imparts a force which drives it down the duct. The plate is a 0.01 m thick plastic with a density of 1270.0 kg/m3. The pressure as a function of distance along the tube is shown in Figure 3-4. The high pressure region must move the plate and compress the air on the low pressure side of the plate.

30

DRDC Atlantic CR 2011-343

Figure 3-4: Pressure distribution at 0.55 ms This case was repeated with the material interface location between the high and low pressure states tracked. There is no structural plate present in this simulation and hence a uniform pressure is observed across the interface, Figure 3-4. Also shown is the solution when the interface is not tracked. While the pressure field is better without tracking, there is diffusion of the materials from one side of the interface into one another. While this is a very simple test case which may not be indicative of realistic simulation performance, the material interface tracking added 50% to the computational run time.

3.3

Casing Expansion

The material interface tracking for the expansion of a high pressure region confined by a high density “casing” was simulated. The two dimensional domain is a 1 m by 1 m box with rigid walls, with a 0.05 m thick rectangular casing region around a 0.2 m by 0.2 m high pressure region centered in the domain, see Figure 3-5. This simulation is meant to demonstrate the material tracking implementation only, and does not represent a realistic test scenario. An ideal gas equation of state was used for all regions. The casing region is a purely fluid region and has no strength. The initial conditions for the regions are given in Table 3-1.

DRDC Atlantic CR 2011-343

31

Figure 3-5: Casing expansion initial regions

Table 3-1: Casing expansion initial conditions Region High Pressure Casing Exterior

Pressure (Pa) 1.01×109 2.98×108 1.01×105

Temperature (K) 3000 385 300

Density (kg/m3) 1177 2700 1.177

The computational grid utilized cells of 0.01 m side length. A quadrilateral grid with 10 000 cells and a triangular grid with 22 842 cells were utilized, but both show similar results. The initially square high pressure region and casing deform and become rounded. The pressure contours on the triangular grid at 0.3 ms have similar levels with the material interface tracking and without. Figure 3-6 shows the shock from the high pressure expansion after it has just reflected from the walls. The material interface case shows a more rounded shock front. The benefit of material interface tracking is to retain sharp resolution of the distinct materials. Figure 3-7 shows the mass fraction of the casing with and without the interface tracking. The non-tracked solution has diffused the interface over approximately 8 cells by 0.3 ms.

32

DRDC Atlantic CR 2011-343

a) Interface tracking (interface shown as bold black lines)

b) Mixed cell formulation

Figure 3-6: Pressure contours at 0.3 ms

a) Interface tracking (interface shown as bold black lines)

b) Mixed cell formulation

Figure 3-7: Casing mass fraction contours at 0.3 ms A conservative scheme will maintain the energy and the mass of each material throughout the simulation. The mass of air and casing in the entire domain can be determined by a summation of the cell density times its volume. Figure 3-8 shows the total mass of air and casing in the

DRDC Atlantic CR 2011-343

33

simulation. There is no variation throughout the simulation, indicating the altered mesh scheme is conservative. The total energy is also constant throughout the simulation.

Figure 3-8: Total mass and casing mass throughout simulation

3.4

Meso Scale Shock

When examining shock propagation from fluids and through the interior of solids, numerical diffusion of the interface between the materials is not desired. This can lead to excessively high grid resolutions in order to accurately capture the interface when material interface tracking is not used. To simulate detonation a linear shock wave propagates through a liquid explosive with an initial density of 1128 kg/m3, and then proceeds through two regions of aluminum, 2785 kg/m3, separated by more liquid explosive. A Mie-Gruneisen (shock Hugoniot) equation of state is used for both the liquid explosive and aluminum regions. 2000 cells were used over the 1×10-4 m domain in order for the simulation without interface-tracking to provide acceptable results. Interaction of the shock with the boundaries between the aluminum and explosive will generate a series of reflected shocks and expansion fans, as illustrated in Figure 3-9.

34

DRDC Atlantic CR 2011-343

Figure 3-9: Density and pressure contours illustrating propagation of shock and rarefaction waves in time with material interface tracked

Figure 3-10: Density and pressure contours illustrating propagation of shock and rarefaction waves in time with the mixed cell formulation DRDC Atlantic CR 2011-343

35

Figure 3-10 shows the same simulation when the material interface is not tracked. While hard to discern due to the high grid density, the material interface in this case is smeared over approximately 20 cells by the 2×10-8 s. When the interface is not tracked, cells can contain more than one material. An iterative process is used in these cells to determine pressure based on the mass fractions and equations of state of the materials. For this test case the computational cost of the iterative process outweighed the cost of the material interface tracking. The non-tracked case took 4.2 times longer to simulate.

3.5

Bubble Expansion

Overpressure

For an underwater explosion the expansion and contraction of the explosive product bubble must be correctly determined. A typical UNDEX bubble cycle is shown in Figure 3-11. Upon detonation the explosive products form a high pressure region. This sends a shock through the water as it expands. Due to the momentum imparted to the water the products and the water near them will be come over-expanded. After reaching its maximum radius the bubble will start to collapse down and form another high pressure region. This process repeats itself with the maximum bubble radius and high pressure falling in amplitude during each cycle.

Time

Shock and Bubble Formation

Overexpansion

1st Collapse

Overexpansion

2nd Collapse

Figure 3-11: UNDEX bubble shape and pressure history A bubble expansion generated by 1.03 kg sphere of TNT (0.05 m radius) is simulated and compared to the analytical correlation given in Cole [31]. The greater the depth, the shorter the cycle time, so depth of 1000 m is used to reduce total simulation time. The grid contained 6500 cells in the radial direction to a radius of 20 m. The TNT is modeled with a Jones-Wilkins-Lee (JWL) equation of state and the water with a stiffened gas equation of state. The stiffened gas parameters are gamma = 4.4 and p∞ = 5.5×108 Pa. The initial conditions are representative of a 1000 m water depth and a TNT density of 1600 kg/m3. The maximum bubble radius vs. time is shown in Figure 3-12.

36

DRDC Atlantic CR 2011-343

Figure 3-12: UNDEX maximum bubble radius Cole gives expressions for the maximum bubble radius, r, and time of pressure peak from the first bubble collapse, t, as:

r = K6

w1/ 3 (d + 33)1/ 3

w1/ 3 t = K5 (d + 33)5 / 6

(28)

Where w is the charge weight (in lb), d is the charge depth (in ft), and the K coefficients for a TNT explosive are K6,TNT = 12.67 ft4/3/lb1/3 and K5, TNT = 4.268 s⋅ft5/6/lb1/3. The maximum radius and time of first pressure pulse are compared in Table 3-2. The maximum bubble radius is well captured by the interface tracking, while the non-tracked simulation has the interface smeared over a region that is more than 5% of the radius. None of the smeared region matches the analytical radius. Both the tracked and non-tracked interface simulations predict a time for the maximum pressure from the first bubble collapse that is longer than the analytical value. However the point at which the pressure in the bubble is greater than the water pressure agrees well with the analytical value for both simulations.

DRDC Atlantic CR 2011-343

37

Table 3-2: UNDEX bubble simulation comparisons

3.6

Analytical

Interface Tracking

Mixed Cell Formulation

Maximum Radius (m)

0.340

0.342

Time of First Pressure Pulse (ms)

6.537

6.51 for pTNT > pwater 6.84 for pTNT,max

0.317 at YTNT = 90% 0.325 at YTNT = 50% 0.334 at YTNT = 10% 6.53 for pTNT > pwater 7.08 for pTNT,max

Bubble Merging

To verify whether interfaces can merge together two cylindrical high pressure air bubble regions, 0.05 m in radius initially separated by 0.35 m, are allowed to expand and interact. The material interface between the high pressure and low pressure regions is tracked. The high pressure regions are initially at 100 atm and 300 K, compared to the ambient 1 atm and 300 K air. A single bubble will expand to a radius of 2.3 m under these conditions. An unstructured grid is used with cell clustering around the initial high pressure regions. The domain is a circular region that extends to a radius of 5 m, with 14 025 cells. This is relatively coarse for a two-dimensional simulation, but indicative of the cell sizes that would be used for a three-dimensional simulation. The initial mass fraction of the high pressure regions and mesh in their vicinity is shown in Figure 3-13.

Figure 3-13: Initial high pressure region mass fraction contours and mesh distribution As the high pressure regions expand they will start to merge once the resolution of the separating low pressure region is less than grid 1 cell in size, as shown in Figure 3-14. For all the FSI and material interface tracking cases there are no features resolved below the grid resolution. Conservation may not be enforced when there is interface merging/separation since some cells may be “orphaned” and not have any valid neighbors to pass their information to. This occurs at the blue region enclosed by the red bubbles in Figure 3-14. Once the interfaces move close enough that all the low pressure cells will switch sides, they will have no low pressure neighbours to pass their conservative values to. This information will therefore be lost. With proper grid resolution the non-conservation effects should be minor since it will be lost from a small number 38

DRDC Atlantic CR 2011-343

of reduced size cells. The high pressure regions in this case will continue to merge and eventually form a single high pressure bubble. Figure 3-15 and Figure 3-16 illustrate the mass fractions of the single bubble with and without interface tracking. It can be seen that interface tracking is critical to retaining sharp interfaces for the level of grid resolution utilized in this case.

Figure 3-14: Merging of high pressure regions (red)

Figure 3-15: Single high pressure region with material interface tracking

Figure 3-16: Single high pressure region without material interface tracking DRDC Atlantic CR 2011-343

39

3.7

Plate Deformation

Previously the deformation of steel plates subject to the load from a RP-83 detonator (1.1 g of RDX) was simulated with Chinook [32]. These simulations utilized the previous large deformation implementation, which did not deform the fluid mesh and was non-conservative. The domain is illustrated in Figure 3-17. The galvanized steel plate is 0.0013 m thick, with a 0.127 m radius. In the trial series [33], the charge standoff and depth of a water-backing layer were varied in order to examine the effect of both on structural failure. The plates exhibited a range of damage and failure patterns, including large-deformation, and rupture.

Figure 3-17: Charge under steel plate schematic For the two cases simulated in this report no water backing thickness was used. The water, air, and explosive were simulated in the CFD domain. The Tait equation of state was used to model the water with a cavitation pressure cut-off at 5000 Pa. Simulations with a stiffened gas equation of state gave similar results. An ideal gas equation of state was used to model the air and explosive products. The initial ambient conditions were used with a pressure of 101 325 Pa for the air and water and an initial pressure of 4.2×109 Pa in the explosive. The structural model is fully clamped at the outer edges. The circular plate was made from 18 gauge galvanized steel, and was modeled using the LS-DYNA material type *MAT_PLASTIC_KINEMATIC. This material type models isotropic and kinematic hardening plasticity with the option of including rate effects. Material properties for the plate included a yield strength of 250 MPa and a failure strain of 0.25, with strain rate effects not considered. It is expected that more accurate results could be obtained with more complete material models, particularly for the ruptured plate results. The plate was initially modeled with quadrilateral elements with edge lengths of 0.001 m at the center of the plate growing to 0.0075 m at the outer edge.

40

DRDC Atlantic CR 2011-343

Two cases were simulated: a 0.02 m standoff where the target plate did not rupture, and a nearcontact 0.001 m standoff where the plate ruptured and formed three petals. It should be noted that this work only investigated the early-time response of the plates, while damage is compared to results obtained following the completion of the experiment. Experimentally it was found that the structural response and damage was driven by the shock, and not by bubble collapse. Dynamic measurements of the wall displacement for an UNDEX-loaded cylinder [34] show that deformation due to bubble collapse decreases rapidly relative to the shock at standoff distance below 0.6R (where R is the maximum free-field bubble radius), hence the deformation is shock dominated. At a standoff distance of 0.2R, the bubble collapse deformation is reduced to 20% of the total deformation. It is therefore likely that for contact charges where the standoff distance is below 0.01R, the influence of bubble collapse is very small. For the 0.001 m stand-off case both the structure and the bubble interface were tracked. The grid was relatively coarse, utilizing 0.001 m cells near the bubble and plate. The mass fractions of the explosive at 0.1398 s are shown in Figure 3-18. The plate has ruptured at this point and the front of the explosive products propagates through the air faster than the water. The vertical velocity contours on the structural plate at 0.1398 s is shown in Figure 3-19. While five petals are formed compared to the experiment which had three, the form of the petal failure is purely a function of structural meshing for this level of modelling detail.

Figure 3-18: Explosive mass fractions (colour), structural plate (black line), and explosive material interface (white line) at 0.1398 s for a 0.001 m stand-off

DRDC Atlantic CR 2011-343

41

Figure 3-19: Structural plate vertical velocity contours at 0.1398 for 0.001 m stand-off The pressure and material interface contours along the symmetry plane for the 0.02 m stand-off case are shown in Figure 3-20 for the previous [32] and current results. For this case only the structure was tracked, not the bubble interface. Contour levels are comparable to the prior result, with more refined gradients due to more grid refinement. One difference to note is that plate in the current result is more deformed than in the prior result. This reduces the pressure experienced by the plate and hence the loading.

a) Prior result [32]

b) Current simulation

Figure 3-20: Pressure (colour) and material interface contours (black lines) for 0.02 m stand-off The maximum plate displacement values are provided in Table 3-3. The current result is further from the experimental value than the prior less rigorous result. This suggests that there may be part of the simulation that does not correspond to the experiment, such as the yield strength which was assumed, or that there are errors in the numerical implementation. Further experimental comparisons are warranted that include comparisons of pressure histories and a time history of the plate displacement.

42

DRDC Atlantic CR 2011-343

Table 3-3: Maximum plate displacements Case Experiment Prior Simulation Current Simulation

DRDC Atlantic CR 2011-343

Maximum Plate Displacement (m) 0.0230 0.0211 0.0185

43

4

Conclusions and Recommendations

An immersed fluid structure interaction (FSI) methodology was demonstrated for the coupling of CFD fluid domains and FE structural domains. A local mesh adaptation method has been implemented for large deformation FSI calculations. In this implementation the mesh adaptation routines are independent of the CFD solver and they are designed to function with commercial FE tools (LS-Dyna and Abaqus). The FSI routines were extended so Chinook can track material interfaces (no FE solver required). A number of issues still remain to be resolved to provide a truly robust large deformation capability. The small fluctuations in cell volume can give rise to large pressure variations for stiff equations of state. This includes water, and for some cases results in unphysical low pressures in the cells surrounding the interface. For some three-dimensional material interfaces bumpy surfaces were generated. A smoothing routine may be required in addition to the fluid feedback, similar to what is used in a number of Lagrangian interface tracking codes. While the material interface velocity calculation has been coded for interface surface tension this feature is not yet implemented. Calculation of the interface curvature, required for the tension, is simplified due to the use of a mesh to define the interface. Validation test cases detailed in this report show the capabilities of the improved FSI and material interface routines; however larger scale test problems need to be performed. Results should be compared to further experimental data to confirm the simulated results are appropriate. The large deformation capability, initially detailed here, can be applied for any FSI interaction, such as applications in underwater explosions, particle break-up, aerodynamic ablation, and structural failure.

44

DRDC Atlantic CR 2011-343

References ..... [1] Link, R., F. Lin, D. Whitehouse, and J. Slater, “CFD Modeling of Close Proximity Underwater Explosions”, 73rd SAVIAC Shock and Vibration Symposium, Newport, Rhode Island, USA, 2002. [2] Gregson, J., R. Link, J. J. Lee, and M. Smith, “Coupled Simulation of the Response of Targets to Close Proximity Underwater Explosions in Two and Three Dimensions”, 77th SAVIAC Shock and Vibration Symposium, Monterey, California, USA, 2006. [3] Gregson, J., T. Dunbar, J. Lee, and R. Link, “Simulation of Structural Failure from Contact Underwater Explosions”, 78th SAVIAC Shock and Vibration Symposium, Philadelphia, Pennsylvania, USA, 2007. [4] Dunbar, T., L. Jiang, J. Szabo, and L. Cheng, “Fully-Coupled CFD-FEA Simulation of Polymer Treated Steel Plates Subject to Close-Proximity UNDEX Loading”, 79th SAVIAC Shock and Vibration Symposium, Orlando, Florida, USA, 2007. [5] Ritzel, D. V., D. Whitehouse, J. Crocker, and T. E. Dunbar, “Response of a Flexible Structure to Shock, Blast, and Seismic Loads”, presented at the 80th SAVIAC Shock and Vibration Symposium, San Diego, California, USA, 2009. [6] Livermore Software Technology Company (LSTC), 7374 Las Positas Road, Livermore, California, USA, 94551, www.lstc.com [7] Simulia, Rising Sun Mills, 166 Valley Street, Providence, Rhode Island, USA, 02909-2499, www.simulia.com [8] Hyman, James M., “Numerical Methods for Tracking Interfaces”, Physica, Vol. 12D, pp. 396-407, 1984. [9] Shyy, W., M. Francois, and H.S. Udaykumar, “Cartesian and Curvilinear Grid Methods for Multi-domain, Moving Boundary Problems”, Thirteenth International Conference on Domain Decomposition Methods, Lyon, France, 2001. [10] Ubbink, Onno, Numerical Prediction of Two Fluid Systems with Sharp Interfaces, Doctoral Thesis, Graduate Department of Mechanical Engineering, University of London, London, UK, 1997. [11] Bailey, D. A., P. K. Sweby, and P. Glaister, “A Ghost Fluid, Moving Finite Volume Plus Remap Method for Compressible Euler Flows”, ICFD conference on Numerical Methods for Fluid Dynamics, Oxford, UK, 2004. [12] Kambouchev, N., L. Noels, and R. Radovitzky, “Nonlinear Compressibility Effects in Fluid-Structure Interaction and their Implications on the Air-Blast Loading of Structures”, Journal of Applied Physics, Vol. 100, No. 6, pp. 063519.1-063519.11, 2006.

DRDC Atlantic CR 2011-343

45

[13] Wren, G. P., S. E. Ray, S. K. Aliabadi, and T. E. Tezduyar, “Simulation of Flow Problems with Moving Components, Fluid-Strucutre Interactions and Two-Fluid Interfaces”, International Journal for Numerical Methods in Fluids, Vol. 24, pp. 1433-1448, 1997. [14] Hirt, C.W., A.A. Amsden, and J.L. Cook, “An Arbitrary Lagrangian–Eulerian Computing Method for All Flow Speeds”, Journal of Computational Physics, Vol. 14, No. 1, pp. 227253, 1974. Reprinted in Vol. 135, No. 2, pp. 198-216, 1997. [15] Souli, M., A. Ouahsine, and L. Lewin, “ALE Formulation for Fluid-Structure Interaction Problems”, Computational Methods in Applied Mechanical Engineering, Vol. 190, pp. 659675, 2000. [16] Ahn, H. T., and M. Shashkov, “Multi-material Interface Reconstruction on Generalized Polyhedral Meshes”, Journal of Computational Physics, Vol. 226, pp. 2096-2132, 2007. [17] Chertock, A., S. Karni, and A. Kurganov, “Interface Tracking Method for Compressible Multifluids”, Mathematical Modelling and Numerical Analysis, Vol. 42, pp. 991-1019, 2008. [18] Sussman, M., P. Smereka, and S. Osher, “A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow”, Journal of Computational Physics, Vol. 114, No. 1, pp. 146-159, 1994. [19] Farhat, C., A. Rallu, and S. Shankaran, “A Higher-Order Generalized Ghost Fluid Method for the Poor for the Three-Dimensional Two-Phase Flow Computation of Underwater Implosions”, Journal of Computational Physics, Vol. 227, No. 16, pp. 7674-7700, 2008. [20] Puckett, E. G., A. S. Almgren, J. B. Bell, D. L. Marcus, and W. J. Rider, “A High-Order Projection Method for Tracking Fluid Interfaces in Variable Density Incompressible Flows”, Journal of Computational Physics, Vol. 130, pp. 269-282, 1997. [21] Pilliod, J. E. Junior, and E. G. Puckett, “Second-Order Accurate Volume-of-Fluid Algorithms for Tracking Material Interfaces”, Journal of Computational Physics, Vol. 199, No. 2, pp. 465-502, 2004. [22] Braeunig, J. P., B. Desjardins, and J. M. Ghidaglia, “A Totally Eulerian Finite Volume Solver for Multi-Material Fluid Flows”, European Journal of Mechanics B/Fluids, Vol. 28, pp. 475-485, 2009. [23] Hawken, Donald, “Design Document for 2D Flow Interface Trackers in AGS Code”, CDL 1013, Combustion Dynamics Limited (now Martec Limited), 2002. [24] De Zeeuw, D., and K. G. Powell, “An Adaptively Refined Cartesian Mesh Solver for the Euler Equations”, Journal of Computational Physics, Vol. 104, No. 1, pp. 56-68, 1993. [25] Sachdev, Jai Singh, Parallel Solution-Adaptive Method for Predicting Solid Propellant Rocket Motor Core Flows, Doctoral Thesis, Graduate Department of Aerospace Science and Engineering, University of Toronto, Toronto, Canada, 2007.

46

DRDC Atlantic CR 2011-343

[26] Aftosmis, M., “Robust and Efficient Cartesian Mesh Generation for Component-Based Geometry”, presented at the 35th AIAA Aerospace Sciences Meeting, Reno, Nevada, USA, 1997. [27] Batten, P., N. Clarke, C. Lambert, and D. M. Causon, “On the Choice of Wavespeeds for the HLLC Riemann Solver”, SIAM Journal of Scientific Computing, Vol. 18, No. 6, pp. 1553-1570, 1997. [28] Roe, P. L., “Approximate Riemann Solvers, Parameter Vectors and Difference Schemes”, Journal of Computational Physics, Vol. 43, pp. 357-372, 1981. [29] Jaegle, F. and V. Schleper, “Exact and Approximate Riemann Solvers at Phase Boundaries”, Analysis, 2010. [30] Ericson, Christer, Real Time Collision Detection, Elsevier Inc., Morgan Kaufmann Publishers, 500 Sansome St., Suite 400, San Francisco, California, USA, 2005. [31] Cole, R. H., Underwater Explosions, Princeton University Press, Princeton, New Jersey, USA, 1948. [32] Gregson, J., T. Dunbar, J. J. Lee, and R. Link, “Simulation of Structural Failure from Contact Underwater Explosions”, 78th SAVIAC Shock and Vibration Symposium, Philadelphia, Pennsylvania, USA, 2007. [33] Slater, J. E., G. Rude, and G. Paulgaard, “Experimental Study of Air-Backed and WaterBacked Targets During Near-Contact Explosions”, DRDC Suffield TR 2005-152, Defence Research and Development Canada - Suffield, 2005. [34] Lee, J. J., J. E. Slater, and G. Rude, “Close-Proximity Blast Loading and Damage of Cylinders”, Proceedings of the 7th Canadian Marine Hydrodynamics & Structures Conference, Halifax, Nova Scotia, Canada, 2005.

DRDC Atlantic CR 2011-343

47

This page intentionally left blank.

48

DRDC Atlantic CR 2011-343

List of symbols/abbreviations/acronyms/initialisms ALE

Arbitrary Lagrangian-Eulerian

CFD

Computational Fluid Dynamics

CFL

Courant–Friedrichs–Lewy

DND

Department of National Defence

DRDC

Defence Research & Development Canada

DRDKIM

Director Research and Development Knowledge and Information Management

FE

Finite Element

FSI

Fluid Structure Interaction

R&D

Research & Development

UNDEX

Underwater explosion

DRDC Atlantic CR 2011-343

49

This page intentionally left blank.

50

DRDC Atlantic CR 2011-343

Distribution list Document No.: DRDC Atlantic CR 2011-343 LIST PART 1: Internal Distribution by Centre 1 1 3

Contract Scientific Authority (1 CD) Liam Gannon, DS, NPSS (1 CD) DRDC Atlantic Library

5

TOTAL LIST PART 1

1

LIST PART 2: External Distribution by DRDKIM Library and Archives Canada, Attn: Military Archivist, Government Records Branch National Defence Headquarters MGen. George R. Pearkes Bldg. 101 Colonel By Drive Ottawa, Ontario K1A OK2

1 1

DGMEPM/DMSS 2 DRDKIM

1

Dr. Julian J. Lee DRDC Suffield PO Box 4000 Stn Main Medicine Hat, AB T1A 8K6

1

J.A.A Vaders Defence Material Organisation Van der Burchlaan 31 PO Box 90822 MPC 58A 2509 LV The Hague, The Netherlands

1

J.E. van Aanhold TNO Built Environment and Geosciences Van Mourik Broekmanweg 6 2628 XE Delft The Netherlands

1

Niklas Alin FOI - Swedish Defence Research Agency Grindsjon Research Centre SE-147 25 Tumba, Sweden

7

TOTAL LIST PART 2

12

TOTAL COPIES REQUIRED

DRDC Atlantic CR 2011-343

51

This page intentionally left blank.

52

DRDC Atlantic CR 2011-343

DOCUMENT CONTROL DATA

(Security classification of title, body of abstract and indexing annotation must be entered when the overall document is classified) 1.

ORIGINATOR (The name and address of the organization preparing the document. Organizations for whom the document was prepared, e.g. Centre sponsoring a contractor's report, or tasking agency, are entered in section 8.)

2.

UNCLASSIFIED (NON-CONTROLLED GOODS) DMC A REVIEW: GCEC JUNE 20110

Martec Limited 1888 Brunswick St, Suite 400 Halifax, Nova Scotia, Canada B3J 3J8 3.

SECURITY CLASSIFICATION (Overall security classification of the document including special warning terms if applicable.)

TITLE (The complete document title as indicated on the title page. Its classification should be indicated by the appropriate abbreviation (S, C or U) in parentheses after the title.)

Large Deformation Fluid-Structure Interaction and Interface Tracking for Chinook: 4.

AUTHORS (last name, followed by initials – ranks, titles, etc. not to be used)

Derrick Alexander 5.

DATE OF PUBLICATION (Month and year of publication of document.)

February 2012 7.

6a. NO. OF PAGES 6b. NO. OF REFS (Total containing information, (Total cited in document.) including Annexes, Appendices, etc.)

64

34

DESCRIPTIVE NOTES (The category of the document, e.g. technical report, technical note or memorandum. If appropriate, enter the type of report, e.g. interim, progress, summary, annual or final. Give the inclusive dates when a specific reporting period is covered.)

Contract Report 8.

SPONSORING ACTIVITY (The name of the department project office or laboratory sponsoring the research and development – include address.)

Defence R&D Canada – Atlantic 9 Grove Street P.O. Box 1012 Dartmouth, Nova Scotia B2Y 3Z7 9a. PROJECT OR GRANT NO. (If appropriate, the applicable research and development project or grant number under which the document was written. Please specify whether project or grant.)

11gq 10a. ORIGINATOR'S DOCUMENT NUMBER (The official document number by which the document is identified by the originating activity. This number must be unique to this document.)

9b. CONTRACT NO. (If appropriate, the applicable number under which the document was written.)

W7707-115146 10b. OTHER DOCUMENT NO(s). (Any other numbers which may be assigned this document either by the originator or by the sponsor.)

DRDC Atlantic CR 2011-343 11. DOCUMENT AVAILABILITY (Any limitations on further dissemination of the document, other than those imposed by security classification.)

Unlimited 12. DOCUMENT ANNOUNCEMENT (Any limitation to the bibliographic announcement of this document. This will normally correspond to the Document Availability (11). However, where further distribution (beyond the audience specified in (11) is possible, a wider announcement audience may be selected.))

Unlimited

13. ABSTRACT (A brief and factual summary of the document. It may also appear elsewhere in the body of the document itself. It is highly desirable that the abstract of classified documents be unclassified. Each paragraph of the abstract shall begin with an indication of the security classification of the information in the paragraph (unless the document itself is unclassified) represented as (S), (C), (R), or (U). It is not necessary to include here abstracts in both official languages unless the text is bilingual.)

The ability to assess the loading, response and survivability of naval structures and equipment to withstand the effects of underwater weapons containing non-ideal explosives is of fundamental importance; particularly for determining the threat damage potential of new weapon technologies, and developing effective counter-measures. Further advancement of the Martec Limited’s Chinook computational fluid dynamics software is detailed, with more advanced methods for material interface tracking and structural large deformation capability for fluid structure interaction (FSI). The local mesh adaptation technique modifies the CFD mesh in the cells adjacent to the structure, such that the cell faces correspond to the location of the structure, allowing for precise determination of pressure loads and interface velocities. Details of the numerical implementation and theory are provided. A series of validation cases were performed in order to verify the FSI and material interface tracking functionality and validate the results. The large deformation capability, initially detailed here, can be applied for any FSI interaction, such as applications in underwater explosions, particle break-up, aerodynamic ablation, and structural failure. The large deformation capability allows for assessment of failure properties and the behaviour of structures under high-rate impulse loading.

14. KEYWORDS, DESCRIPTORS or IDENTIFIERS (Technically meaningful terms or short phrases that characterize a document and could be helpful in cataloguing the document. They should be selected so that no security classification is required. Identifiers, such as equipment model designation, trade name, military project code name, geographic location may also be included. If possible keywords should be selected from a published thesaurus, e.g. Thesaurus of Engineering and Scientific Terms (TEST) and that thesaurus identified. If it is not possible to select indexing terms which are Unclassified, the classification of each should be indicated as with the title.)

Chinook, Computational Fluid Dynamics, CFD, underwater explosion, UNDEX, blast, shock, gas bubble dynamics, fluid-structure interaction

This page intentionally left blank.