A high-fidelity geometric contact method for thin fabrics using shell finite elements in explicit and implicit time scheme

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content A high-fidelity geometric contact method for thin fabrics using shell finite elements in explicit and implicit time scheme by Qingquan Wang A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (AEROSPACE ENGINEERING) August 2024 Copyright 2024 Qingquan Wang Acknowledgments I would like to extend my deepest gratitude to my supervisor, Dr. Carlos Pantano-Rubino, for his unwavering support, insightful guidance, and invaluable feedback throughout this journey. I could not have reached this point without his endless help and support. I am also thankful to my colleagues, Hang Yu and Aditya Heru Prathama, whose discussions have greatly enriched this research and made the process enjoyable. My heartfelt thanks go to my parents and family, who instilled in me the values of perseverance, resilience, and adaptability, and, most importantly, the importance of choosing the right path over the easy one. These qualities have been crucial to the success of my PhD research. Your unconditional love, patience, and encouragement have been a constant source of strength and motivation. I am deeply appreciative of my USC friends Anguo Hu and Weiye Wang, who always generously shared their snacks after long days of work. Additionally, I am grateful for the support and camaraderie of my friends Huairen Zhou, Yao Du, Lucy Chen, Cathy Lan, Melodia Cui, Xiaoyan Lin, and Roger Chen, whose companionship and understanding during gatherings have provided much-needed relief and inspiration. This work is dedicated to all those who have supported me thus far. Achieving a PhD is a collective effort that relies on the dedication and encouragement of a supportive network. I am profoundly thankful for your contributions and belief in me. ii Table of Contents Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Chapter 2: Thin shell formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Chapter 3: Contact mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Non-smooth Lagrange mechanics . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Augmented Lagrange multiplier method . . . . . . . . . . . . . . . . . . . . 24 Chapter 4: Geometric contact method in explicit time integration . . . . . . . . . . . . 30 4.1 Two phase contact detection algorithm . . . . . . . . . . . . . . . . . . . . . 30 4.1.1 Time marching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.1.2 Parallel implementation . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2 Preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2.1 Hemispherical shell collision . . . . . . . . . . . . . . . . . . . . . . . 40 4.2.2 Deployment protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2.3 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2.4 Spherical shell bounces off a rigid wall . . . . . . . . . . . . . . . . . 50 Chapter 5: Geometric contact method in implicit time integration . . . . . . . . . . . 66 5.1 Signed gap function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 iii 5.2 Time marching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2.1 Master iteration with augmented Lagrange multipliers . . . . . . . . 71 5.2.2 Iteration Accelerator with classical Lagrange multipliers . . . . . . . 74 5.3 Contact detection algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.4 Preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.4.1 Triangle-triangle impact . . . . . . . . . . . . . . . . . . . . . . . . . 80 Chapter 6: Programming techniques of contact algorithms . . . . . . . . . . . . . . . . 82 6.1 Data structures and algorithms of geometric contact entities . . . . . . . . . 83 6.2 Template library and user-centric approach . . . . . . . . . . . . . . . . . . . 86 Chapter 7: Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.1 Simulation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.1.1 Fabric folding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 7.1.2 Simple parachute folding . . . . . . . . . . . . . . . . . . . . . . . . . 98 7.1.3 Parachute deployment . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Chapter 8: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 iv List of Tables v List of Figures 1.1 BVH search demonstration. (a) Primitives with dashed bounding boxes are grouped by proximity. The sphere and equilateral triangle share a solid bounding box, which, along with others, is enclosed by a larger solid bounding box for the entire scene. (b) The hierarchy: the root node covers the entire scene, with one child bounding the sphere and triangle, and another bounding the skinny triangle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 An example of 2D spatial subdivision of four objects. . . . . . . . . . . . . . 4 1.3 An Example of the sort and sweep algorithm of three objects. . . . . . . . . 4 2.1 Thin shell geometry in the reference (undeformed) and deformed configurations in curvilinear coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2 A finite element mesh of a spherical shell. . . . . . . . . . . . . . . . . . . . 20 4.1 Mesh distance to a element in the center (marked as 0). . . . . . . . . . . . . 32 4.2 Culling CABB pairs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.3 Sketch of vertex/triangle (left) and edge-edge (right) intersection problems. . 33 4.4 3 MPI task for FEM and contact search partition. . . . . . . . . . . . . . . . 37 4.5 Contact detection in recursive time advancement scheme in parallel. . . . . . 39 4.6 Two hemispheres simulations at t = 0 in head-on collisions (left) and rimcollisions (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.7 Simulation results over time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.8 Main parameters of drogue parachute constructed geometry. . . . . . . . . . 43 vi 4.9 General inflation pressure timeline evolution. . . . . . . . . . . . . . . . . . . 45 4.10 Results for 1-band loading scenario. From top to bottom: at the end of stage 1, at the end of stage 2, after disreefing, where t = t∗ +t1, t∗ +t1 +t2, t∗ +t1 + t2 + 0.00703 s, respectively. From left to right: first no overloading; then full overloading and burst overloading at po = 10, 20, and 40 kPa, respectively. . 46 4.11 Results for 2-band loading scenario; arrangement as in Figure 7.14. . . . . . 46 4.12 Results for 3-band loading scenario; arrangement as in Figure 7.14. . . . . . 47 4.13 Histograms of von Mises stress (σv) for 1-band loading in linear-log coordinates: lower and higher values of σv are 0.1 MPa and 104 MPa, respectively; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.14 Histograms of von Mises stress (σv) for 2-band loading in linear-log coordinates; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . 49 4.15 Histograms of von Mises stress (σv) for 3-band loading in linear-log coordinates; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . 49 4.16 Spherical shell-wall collision simulation: front view (left); top view (right). . 51 4.17 Energy partition over time during impact for various material stiffness. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.18 Impulse histogram for different materials. . . . . . . . . . . . . . . . . . . . . 53 4.19 Energy conservation for various γ of material stiffness E = 2.1 × 106 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.20 Energy conservation for various γ of material stiffness E = 2.1 × 107 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.21 Energy conservation for various γ of material stiffness E = 2.1 × 108 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.22 Energy conservation in one bounce off for various mesh resolutions. . . . . . 57 4.23 Energy conservation in one bounce off for various mesh resolutions. . . . . . 57 4.24 Energy conservation for a medium fine mesh in one bounce off. . . . . . . . . 58 4.26 Mesh configurations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 vii 4.27 Energy conservation for various mesh configuration of material stiffness E = 2.1 × 106 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.28 Energy conservation for various mesh configuration of material stiffness E = 2.1 × 107 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.29 Energy conservation for various mesh configuration of material stiffness E = 2.1 × 108 Pa. (see in color print) . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.1 vertex-triangle (left) and edge-edge (right) intersections. The gap function in both cases are modeled as intersected tetrahedron volumes, whose sides are highlighted in yellow dashed lines, formed by intersected vertices. . . . . . . 69 5.2 Possible intersection cases. Only cases IV and V are considered true intersections. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.3 Initial setup for head-on impacting (left) and sliding (right) respectively. . . 80 5.4 Simulation results of impacting triangle. . . . . . . . . . . . . . . . . . . . . 81 5.5 Simulation results of sliding triangle. . . . . . . . . . . . . . . . . . . . . . . 81 6.1 UML diagram of contact solver. . . . . . . . . . . . . . . . . . . . . . . . . . 86 6.2 UML diagram of driver class inheritance. . . . . . . . . . . . . . . . . . . . . 87 7.1 fabric-folding simulation at t = 0 with front view on the left, and side view on the right. The bend has an inner radius of 3.45 mm. . . . . . . . . . . . . 90 7.2 fabric mesh partition for a 32-processor MPI task. . . . . . . . . . . . . . . . 90 7.3 Results of fabric-folding simulation. . . . . . . . . . . . . . . . . . . . . . . . 92 7.4 Immersive contact happens when two sides of material touches its inner layer simultaneously at 0.016 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.5 Number of contact pairs (logscale) for global and local phase looking at different time windows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 7.6 CPU wall time (logscale). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 viii 7.7 Cross-processor distribution of global contacting pairs for case a (averaged over 10 steps). (see in color print) . . . . . . . . . . . . . . . . . . . . . . . 96 7.8 Cross-processor distribution of local contacting pairs for case a (averaged over 10 steps). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 7.9 The configuration of a simple parachute with cables. . . . . . . . . . . . . . . 98 7.10 Temporal development of cable force. . . . . . . . . . . . . . . . . . . . . . . 99 7.11 Simulation results of simple parachute folding driven by cable force. . . . . . 100 7.12 Main parameters of drogue parachute constructed geometry. . . . . . . . . . 102 7.13 General inflation pressure timeline evolution. . . . . . . . . . . . . . . . . . . 104 7.14 Results for 1-band loading scenario. From top to bottom: at the end of stage 1, at the end of stage 2, after disreefing, where t = t∗ +t1, t∗ +t1 +t2, t∗ +t1 + t2 + 0.00703 s, respectively. From left to right: first no overloading; then full overloading and burst overloading at po = 10, 20, and 40 kPa, respectively. . 105 7.15 Results for 2-band loading scenario; arrangement as in Figure 7.14. . . . . . 105 7.16 Results for 3-band loading scenario; arrangement as in Figure 7.14. . . . . . 106 7.17 Histograms of von Mises stress (σv) for 1-band loading in linear-log coordinates: lower and higher values of σv are 0.1 MPa and 104 MPa, respectively; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.18 Histograms of von Mises stress (σv) for 2-band loading in linear-log coordinates; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . 108 7.19 Histograms of von Mises stress (σv) for 3-band loading in linear-log coordinates; arrangement as in Figure 7.14. . . . . . . . . . . . . . . . . . . . . . . 108 ix Abstract While numerical models for contact mechanics are increasingly common, achieving efficient geometric contact resolution with a robust search strategy remains challenging. This research addresses this issue by introducing a comprehensive solution that includes an exact geometric contact mechanics algorithm for thin shell finite elements, applicable to both explicit and implicit time schemes. In the explicit time scheme, the method offers precise geometric resolution of self-contact interactions through a sub-time-step marching technique, employs adaptive data structures to minimize computational overhead, and features a specialized parallelization approach with load-balancing capabilities. To manage the polynomial time complexity, an efficient detection algorithm divides the problem into global and local contact detection phases. Impact equations are then used to resolve contact events while preserving kinematic energy and momentum. This algorithm is integrated with an MPI-based parallelization framework for thin-shell finite element solvers, ensuring balanced load distribution. The research also investigates the integration of the two-phase contact resolution approach within the implicit time scheme. Here, the use of a signed gap function provides precise contact constraints, offering an alternative solution for contact-intensive cases where explicit methods may be time-consuming due to sub-time step marching. The robustness and accuracy of the algorithm are validated through three numerical studies, and strong scaling analysis demonstrates the scalability of the parallelization approach. x Chapter 1 Introduction 1.1 Background Large deformations in structural mechanics often lead to scenarios where a material may come into contact with itself, a phenomenon known as self-contact. This presents significant challenges in accurately modeling and simulating such behavior, particularly when the material in question is thin and flexible, like fabric or other thin-shell structures. The inherent complexity of self-contact, combined with the nonlinearity of large deformations, requires sophisticated computational techniques to ensure accurate and efficient simulations. Addressing these challenges is crucial for applications in engineering structures to biomechanics, where understanding the material behavior under large deformations is essential. With this in mind, a contact dynamics formulation, in its simplest and most comprehensive form, involves three key problems: 1) identifying potential contact points between colliding bodies in a multibody system, 2) evaluating the contact-impact forces, and 3) managing the transitions between contact and non-contact scenarios to reach the equilibrium states. In the finite element frameworks, these steps involve searching for contacting elements over N distinct elements, and handling the contact responses by applying contact forces on these elements. The contact forces further deform the shell. As a result of the deformation, some parts of the previous contact pairs may separate, or other parts of the shell may come 1 into contact. Neither contacting elements nor contact surfaces are unknown at the beginning of the simulation. They are typically determined in an iterative approach, where the stable finite element solution is found by repeatedly searching for contacting elements and applying contact forces. Extensive research has been devoted to exploring solutions to contact problems, by proposing efficient, reliable, and accurate solutions or modeling methods in all the aforementioned areas [1, 2, 3, 4, 5]. For contact detection, the most intuitive and expensive search method is brute force, which enumerates all possible element-element pairs in the mesh and checks whether each pair is in contact. Some early explorations in the 80s to 90s on this topic attempted to reduce this natural polynomial complexity, O(N2 ), by decomposing the problem into a two-phase contact detection procedure [2, 1, 3, 4, 5]. This is accomplished by identifying elements in the mesh that are potentially in contact with each other (global contact detection), followed by a detailed contact check to determine which candidate elements are actually in contact (local contact detection). The global contact detection was often referred to as neighborhood identification [2, 1]. Numerous collision detection algorithms have been implemented and employed. Here, we give a brief summary of some of the most fundamental approaches, namely spatial subdivision, coordinate sorting, and bounding volume hierarchy (BVH). Other existing methods can be largely regarded as variants and derived versions of one or more of these fundamental methods [4]. BVH and coordinate sorting methods require constructing a bounding box of the occupied space of a single element or grouped elements from time t to t + dt. 2 Figure 1.1: BVH search demonstration. (a) Primitives with dashed bounding boxes are grouped by proximity. The sphere and equilateral triangle share a solid bounding box, which, along with others, is enclosed by a larger solid bounding box for the entire scene. (b) The hierarchy: the root node covers the entire scene, with one child bounding the sphere and triangle, and another bounding the skinny triangle. The spatial subdivision technique partitions the 3D domain occupied by the considered mesh into a uniform grid, such that each grid cell is at least as large as the biggest element in the mesh. At the beginning of the simulation, a hash table is established, which maps each cell to the elements whose centroids are contained within the cell. Collision tests are thus only performed among elements that belong to the same cell or adjacent cells. In the case of four objects contact search in 2D space (see Figure 1.2), a collision test is performed between objects o1 and o2 as both have centroids in cell 1, and between o2 and o3 as both appear in cell 5 with o3’s centroid there. No test is performed between o2 and o4, despite both appearing in cell 2, as neither has its centroid in that cell. 3 Figure 1.2: An example of 2D spatial subdivision of four objects. Coordinate sorting, often referred to as sort and sweep, projects the bounding volume of each object onto the x, y, and z-axis and forms a list of the one-dimensional interval [bi , ei ] along each axis, where bi marks the beginning of the bounding box interval, and ei marks its end. The objects can only possibly contact each other if their intervals of bounding boxes on all three axes overlap. In this method, a sorted list of all the beginning markers bi and end marker ei for each axis is kept in ascending order. To find out the intersecting pairs of objects, the list is traversed from beginning to end, and whenever a bi marker is found in the list, object i is added to a separate set of active objects; whenever a ei is found in the list, object i is removed from the set (see a demonstration in Figure 1.3). This is followed by a contact check among all possible pairs of objects in the set whenever a new object is added. Figure 1.3: An Example of the sort and sweep algorithm of three objects. 4 BVH method, as its name implies, builds up a hieratical structure, such as an octree or quadtree, from the bounding volumes of objects. The root node usually contains the entire domain occupied by all the objects. Each of its children nodes takes part in the domain. Then the domain is further subdivided into its grandchildren until the leaf nodes contain the domain on the order of the size of the averaged bounding volume. The merit of this method is that it allows fast pruning to quickly filter out the subtrees whose corresponding bounding volumes do not intersect with the object of interest. Similar to BVH method, spatial subdivision and coordinate sorting can also employ some query acceleration strategies to achieve an average O(NlogN) time complexity by using a particular tree structure [6, 7]. At the end of the broad phase contact detection, each element obtains a neighbor list of elements of potential collisions. In practice, a shortcut is sometimes used to simplify the process, which assumes that only small changes associated with the mesh movements take place within one or few time steps for a properly chosen time step [2]. Consequently, this list stays unchanged until the next global search is needed. Following the global contact detection, the exact contact point/location usually needs to be determined to compute impact forces. This process is referred to as local contact detection. In this phase, where a detailed contact check is performed on each element against members in its neighbor list to select the true contacting pairs. Here, the two most common approaches are introduced: one is the penalty-based method [4, 8, 9], and the other involves solving the co-planarity equations from exact geometrical considerations [2]. The penalty-based method first examines the actual intersection among objects. For an unstructured finite element mesh, this becomes an element-element penetration test by typically checking if one of the vertices of one element lies within the interpenetration side of the other element. Then the interpenetrated vertex is projected back to its nearest surface, applied with the proper repulsion force, whose magnitude is proportional to the intersected volume. This method works well with rigid-body problems, like robot motion-planning, particles collisions, structure modeling, etc, but encounters failures in simulating large de5 formations of thin shell materials. The critical reason for failures here is that a unique side of thin shell material, modeled as a surface mesh with zero thickness, cannot be defined. Furthermore, the back-projection point for an interpenetrated vertex is arbitrary when multiple instantaneous contacts happen around that point [4, 9]. Alternatively, the geometric method does not admit any interpretation by predicting the exact contact time and location. In this process, a set of co-planarity nonlinear polynomial equations are solved iteratively by Newton’s method. Using the gap volume or distance function, the repulsion force is determined and applied to the contacting triangles whose contact time is within the small step [2]. Clearly, this method applies not only to rigid body collision problems but also to thin shell contact mechanics. However, its capability to solve problems of multiple instantaneous contacts, like fabric folding, has not been clearly revealed yet in the field of computational mechanics. The study of contact-impact force evaluation began gaining traction in the 1980s and 1990s, spurred by advances in computational mechanics and the increasing need to model complex deformations in engineering and biomedical applications. Early work focused primarily on the theoretical foundations and basic numerical methods for handling contact problems. One of the seminal works in this area is K.L. Johnson’s ”Contact Mechanics” (1985), which provides a comprehensive theoretical framework for understanding the mechanics of contact and impact between elastic bodies [10]. Johnson’s work laid the groundwork for subsequent studies by detailing the behavior of materials under contact conditions and introducing fundamental concepts such as Hertz contact theory. Here, contact formulations can be broadly categorized into two types: regularized formulations and non-smooth formulations. Regularized methods model the transition from non-contact to contact situations using continuous functions, providing a smooth representation of the system’s kinematic quantities. Simplest examples here are Hooke’s law and Hertz contact theory. In regularized approaches, the exact contact point between contacting bodies does not coincide precisely, leading to the existence of numerous potential or candidate contact points. The actual con6 tact point is determined by the one with the maximum indentation. In contrast, non-smooth formulations exhibit abrupt discontinuities in these quantities when collisions occur since the contacting parts are considered rigid. The unilateral constraints are utilized to prevent the local interpenetration and to determine the impulses. At the same time, once the contact force is determined, incorporating it into the system and solving the resulting set of equations is another significant topic that attracts considerable interests. In high level, the resulting equations can be regarded as constrained problems. For regularization methods, the most common approaches to solve these systems include the Lagrange Multiplier method, the penalty method, and the augmented Lagrange Multiplier method. These methods are crucial in enforcing contact constraints within finite element analyses. Classical Lagrange multiplier method typically includes a vector λ as surface contact forces. This method proceeds by treating λ as unknowns and solving the equations of motions as well as constraint functions simultaneously. The resulting equations for a small displacement problem where internal forces are independent of strain rate and proportional to displacement can then be either directly by time integration or in an iterative approach. Early contributions include the work in [11], where a contact resolution with the Lagrange Multiplier method that effectively incorporated these techniques into finite element models to handle contact problems is developed. Later, a forward increment Lagrange multiplier formulation is integrated into explicit time integrator, where displacement constraints at tn+1 are related with Lagrange multipliers at time tn [12]. The penalty method offers an alternative approach to enforcing contact constraints by incorporating a penalty term into the system’s potential energy. This term penalizes any violation of the contact constraints, effectively discouraging interpenetration between bodies. The penalty method does not require additional variables, simplifying the system of equations compared to the Lagrange Multiplier method. However, the method’s effectiveness largely depends on the choice of the penalty parameter; a high value ensures strict enforcement of 7 the contact constraints but can lead to numerical stiffness, while a low value may result in significant interpenetration. Despite its advantages, the penalty method can struggle with stability and accuracy in highly dynamic or complex contact problems. An extensive discussion of this method can be found in [13]. The augmented Lagrange multiplier method combines the strengths of the Lagrange multiplier and penalty methods, addressing their respective weaknesses. This method introduces a penalty term and Lagrange multipliers, augmenting the potential energy of the system. The penalty term ensures approximate enforcement of the contact constraints, while the Lagrange multipliers refine the solution to meet the constraints precisely. This approach balances the computational efficiency of the penalty method with the accuracy of the Lagrange Multiplier method. By iteratively updating the multipliers and penalty parameters, the augmented Lagrange Multiplier method achieves stable and accurate enforcement of contact constraints. This method is particularly effective in handling large deformations and complex contact interactions, making it suitable for various applications in engineering and biomechanics. Although multiple techniques have been proposed to solve these equations, the issue of non-linearity presents a unique set of challenges that must be tackled. In [14], a method to handle elastoplasticity with contact was introduced, emphasizing the complexities introduced by nonlinear material behavior and large deformations. Similarly, [12] explored the challenges associated with modeling contact mechanics in transient nonlinear finite element problems, including issues of numerical stability and convergence in grouped equations in the Lagrange multiplier method involving large deformations and contact interactions. The turn of the millennium saw significant advancements in computational power and numerical methods, enabling more sophisticated and efficient modeling of self-contact phenomena. Enhanced numerical techniques became a focal point of research. A comprehensive review of contact modeling methods is summarized in [15], highlighting the evolution of numerical techniques and their application to various engineering problems. Afterward, a 8 decomposition-based contact handling method, which improved the efficiency and robustness of simulations involving complex contact scenarios, was developed in [9]. In parallel, finite element methods and reduced-order models saw considerable development. For example, meshfree methods for contact problems provide an alternative to traditional finite element methods and address some of their limitations [16]. Reduced-order models, in other cases, simplify simulations without losing the degree of freedom at contact interfaces, particularly for large-scale problems involving many deformable bodies [17]. More recent research has focused on addressing the remaining challenges in material selfcontact, such as computational efficiency, promoting advanced parallel computing methods, and incorporating advanced material models. Here, advanced contact algorithms and data structures have been a major area of focus. To name a few, a semi-explicit algorithm for multibody contact systems reduced computational costs by combining explicit and implicit methods [18]. An implicit material point method eliminated the need for conformal mesh construction, simplifying the process of contact pair identification and handling [19]. Parallelization and computational efficiency have been critical areas of progress. Addressing load balancing concerns in computational mechanics, particularly in the context of contact detection in transient dynamics, high-level dynamic load balancing strategies are proposed to tackle these challenges [20]. More recent implementations of the parallelization advancements include asynchronous contact mechanics [21] and [22], utilization of GPU for multithreading features [23] and [24]. Finite element contact mechanics is a critical field within computational mechanics, focusing on the study of deformations, stresses, and contact forces between interacting surfaces. This domain has broad applications, ranging from engineering to biomedical sciences, and has evolved significantly with advances in computational power and numerical techniques. This section reviews notable works and advancements in finite element contact mechanics, discussing foundational theories, numerical methods, algorithmic developments, and practical applications. 9 The foundational theories of contact mechanics were established by Hertzian contact mechanics in the late 19th century, providing solutions for elastic bodies with curved surfaces in contact. This classical theory, though simplistic, laid the groundwork for more complex models that consider various factors like plastic deformation, surface roughness, and adhesion. Later, the Johnson-Kendall-Roberts (JKR) contact model expanded upon Hertzian contact mechanics by incorporating adhesive forces, which is crucial for understanding microscale interactions. The Derjaguin-Muller-Toporov (DMT) model, introduced later, further refined the understanding of adhesive contact mechanics by considering different deformation regimes. A detailed review of Hertzian contact mechanics, JKR model, and DMT model can be found in [25]. The finite element method (FEM) has become a dominant computational technique for solving contact mechanics problems due to its versatility and robustness. FEM discretizes a continuous domain into smaller finite elements, making it possible to solve complex problems involving large deformations, non-linear material properties, and intricate geometries. In [26], a seminal work on the finite element method provided a comprehensive foundation for its application in various fields, including contact mechanics. In the early developments of FEM for contact mechanics, researchers faced challenges in accurately modeling contact interactions. One significant issue was the detection of contact, where determining whether two bodies are in contact or not is critical. Algorithms for contact detection, such as the node-to-surface and surface-to-surface approaches, were developed to address this. The node-to-surface method, introduced in [27], was one of the first to be implemented, but it faced limitations in handling complex contact interactions and multiple contact points. The surface-to-surface contact algorithm, proposed by [28], improved upon the nodeto-surface method by allowing more accurate and stable contact force calculations. This method considers the entire surface of contact bodies, reducing the likelihood of numerical instabilities and improving convergence rates. Furthermore, penalty methods and Lagrange 10 multipliers have been employed to enforce contact constraints. The penalty method introduces a penalty parameter to prevent penetration between contact bodies, while the Lagrange multiplier method enforces contact conditions exactly but requires solving additional equations, as discussed by [29]. One of the major advancements in finite element contact mechanics is the development of efficient and robust algorithms for large-scale problems. The augmented Lagrangian method combines the advantages of both penalty and Lagrange multiplier methods. This approach iteratively adjusts the Lagrange multipliers and penalty parameters to enforce contact constraints more accurately. The method has been widely adopted in commercial finite element software due to its robustness and efficiency. Another critical area of research is the treatment of frictional contact. Coulomb friction, which is the simplest and most widely used friction model, poses significant challenges for numerical implementation due to its non-linearity and discontinuity. Researchers like [30] have developed algorithms to handle frictional contact within the finite element framework, ensuring convergence and stability. These algorithms often involve iterative schemes that adjust the frictional forces until equilibrium conditions are satisfied. Advances in computational power and parallel computing have enabled the solution of increasingly complex contact mechanics problems. The development of parallel algorithms for FEM has significantly reduced computation times, making it feasible to analyze largescale problems with millions of degrees of freedom. Work in [31] on domain decomposition methods has been pivotal in this regard. These methods divide the computational domain into smaller subdomains, which can be solved independently and in parallel, thus enhancing computational efficiency. The integration of FEM with other numerical methods has also been explored to tackle specific challenges in contact mechanics. For instance, coupling FEM with the boundary element method (BEM) has been shown to effectively handle problems involving infinite or semi-infinite domains. Researchers such as [32] demonstrated the effectiveness of this 11 approach in elastohydrodynamic lubrication problems, where the BEM handles the infinite domain while the FEM models the contact region. Experimental validation remains a crucial aspect of finite element contact mechanics. Techniques such as digital image correlation (DIC) and scanning electron microscopy (SEM) are employed to measure contact deformations and stresses experimentally. These experimental data are essential for validating and calibrating finite element models. Studies in [33] have demonstrated the effectiveness of DIC in capturing detailed deformation fields, providing a robust means for validating numerical models. Finite element contact mechanics has found applications in various engineering fields. In the automotive industry, it is used to analyze the contact interactions in components such as tires, brakes, and seals. The work in [34] on tire-road interactions has been instrumental in understanding the contact mechanics of tires, leading to improved tire designs and better performance. In the aerospace industry, finite element analysis is employed to study the contact interactions in components like turbine blades and landing gear, where accurate modeling of contact stresses is crucial for ensuring structural integrity and reliability. In the biomedical field, finite element contact mechanics is used to study the interactions between biological tissues and medical devices. For example, the design of joint replacements and dental implants heavily relies on accurate modeling of contact stresses to ensure their longevity and performance. Studies in [35] on joint contact mechanics have provided valuable insights into the biomechanical behavior of joints, informing the design of more effective and durable joint prostheses. Microelectromechanical systems (MEMS) represent another area where finite element contact mechanics is crucial. The performance and reliability of MEMS devices, such as sensors and actuators, depend on the accurate modeling of contact interactions at micro and nanoscale levels. Research in [36] has addressed issues such as stiction, adhesion, and wear in MEMS, leading to improved device designs and materials that enhance their performance and reliability. 12 Emerging technologies continue to push the boundaries of finite element contact mechanics. One such area is soft robotics, where the contact mechanics of soft materials is crucial for developing compliant and adaptable robots. Research in [37] has explored the behavior of soft materials under contact, paving the way for innovative robotic designs that can safely interact with humans and delicate objects. Nanotechnology is another frontier where finite element contact mechanics plays a vital role. The manipulation and assembly of nanoscale components rely on precise control of contact forces and deformations. Studies by [38] have demonstrated the importance of understanding contact mechanics in the development of nanoscale devices and systems. Finite element contact mechanics is a dynamic and evolving field with significant theoretical, computational, and practical advancements. From the foundational theories of Hertz and the development of robust numerical algorithms to the application in diverse engineering and biomedical domains, this field continues to address new challenges and enable technological innovations. The integration of advanced computational techniques, experimental validation, and interdisciplinary applications ensures that finite element contact mechanics will remain a cornerstone of modern engineering and scientific research. As computational capabilities expand and new materials and technologies emerge, the scope and impact of finite element contact mechanics are poised to grow even further, driving progress and innovation across various industries and scientific disciplines. With the progress in artificial intelligence and high-performance computing, ongoing research has started focusing on exploring machine learning for predictive modeling and optimizing contact algorithms. This integration aims to lower computational costs and enhance accuracy in simulations involving highly nonlinear and dynamic contact interactions. Utilization of an artificial neural network to predict contact states, locations, and penetration depths, improves computational efficiency by separating the local contact search from dynamic contact analysis [39]. [40] investigates the use of machine learning techniques to predict structural mechanical behaviors in contact scenarios, leveraging the ability of neural 13 networks to map complex functions with minimal assumptions. Finally, continuing to improve geometric contact detection and resolution methods, particularly for materials with highly complex and nonlinear geometries, is essential. These advancements will ensure accurate and stable simulations even in challenging scenarios, paving the way for more sophisticated and reliable modeling of material self-contact. From the foundational work in the 1980s and 1990s to the advanced computational techniques of today, the field of material self-contact has seen substantial progress. While significant challenges remain, particularly in handling complex geometries and nonlinear behaviors, ongoing research promises to further enhance our ability to accurately model and simulate these phenomena, with wide-ranging applications across structural mechanics, bio-mechanics, material science, robotics and beyond. 1.2 Motivations Describing contact accurately poses significant challenges due to the inherent geometrical nonlinearity of deformations. This nonlinearity can greatly increase the computational cost of simulations, especially in large-scale scenarios involving numerous deformable bodies. Achieving accurate modeling of contact nonlinearity often necessitates high-resolution meshes, leading to a high number of degrees of freedom in the simulation. Moreover, discontinuous contact forces can lead to numerical instabilities and convergence issues in simulations [12, 41]. In addition, proper enforcement of contact constraints are crucial in accurately capturing the non-smooth contacting surfaces and computing the contact forces between them. Here, an associated challenge is to determine the exact contact point when it comes to a detailed contact detection, i.e. closest point projection. Existing approaches typically focus on interaction among slave nodes and master surfaces, which requires a two-pass search due to the symmetry nature of contact [2]. Theoretical studies have been focused on the existence 14 and uniqueness to the solution of this problem and corresponding solution strategies are discussed to construct a valid projection space for C 0 and C 1 continuous surfaces in general [27, 42]. Other heuristic research, dealing with certain contact cases when C 1 continuity is broken, has also been presented in [43, 44, 28]. In these studies, the gap function, or distance function, is usually defined as the normal distance between the master surface and the slave node, and therefore, the contact constraints are modeled as extremal problems to minimize the gap function. These methods are suitable and efficient for bounded volumes since the sign of gap function is well-defined, which is crucial for enforcement determination in most state-of-art contact algorithms. Or in other cases, regardless of the sign of the gap function, a back projection is applied at the end of time step to enforce the geometrical non-interpenetration conditions [9]. However, thin materials, such as fabrics, pose a more challenging problem due to the possibility of three-dimensional open surfaces, which lack a pre-defined interior and exterior [45]. Consequently, determining whether interpenetration actually occurs becomes challenging. In situations involving highly nonlinear geometries, distinguishing between ”interpenetration” and ”non-interpenetration” sides can be ambiguous, particularly in cases where multiple layers may become stuck together. It is typically impossible to reverse an incorrect deformation once the surface interpenetrates. As a result, the back-projection approach is unsuitable here because the projection point is not clear. To circumvent some of the difficulties mentioned above, a high-fidelity geometric contact method is presented to simulate thin fabric dynamics represented by shell finite elements both in an explicit and implicit time integration scheme. The primary focus of this paper is predominantly centered on the computational and algorithmic intricacies associated with geometric contact detection and resolution. In this study, we focus exclusively on analyzing elastic impacts without frictional effects. This deliberate choice stems from our primary objective, which centers on investigating computational algorithms for implementing geometric contact detection and resolution. It’s important to note that the presence or absence of 15 friction does not alter the underlying contact detection algorithms being studied. While there exists a plethora of reputable sources delving into the physics of contact mechanics, our contribution does not aim to introduce novel insights in this domain. Instead, this research seeks to address what we perceive as a notable lack of detailed information concerning the implementation of precise geometric contact methods. Numerous excellent references, as mentioned above, elucidate various solution strategies, yet they often lack comprehensive insights into the practical implementation of these algorithms. At this point, the novelties and contributions of the new method are summarized as follows, i) precise geometrical resolution of self-contact interactions enabled by a sub-time-step marching method that recursively monitors possible contact events, ii) use of advanced data structures to minimize computational overhead, and iii) a dedicated parallelization implementation with dynamic load-balancing capability. 1.3 Outline The rest of this thesis is organized as follows: a brief overview of the fundamental work from [46, 47] on thin shell mechanics is provided in Chapter 2. This is followed by an in-depth exploration of contact mechanics applicable to both explicit and implicit time schemes in Chapter 3. The core contribution of this research, detailing a contact algorithm to ensure non-interpenetration for explicit time integration, is presented in Chapter 4, with corresponding methods for implicit time integration discussed in Chapter 5. Chapter 6 covers advanced data structures for programming efficiency and high-level object-oriented programming (OOP) design to facilitate user customization. The dissertation then validates the correctness and performance of the contact resolution via comprehensive simulation results demonstrating the validity and accuracy of the contact algorithm across both time integration schemes in Chapter 7. Finally, the dissertation concludes with a brief conclusion and discussion in Chapter 8. 16 Chapter 2 Thin shell formulation The contact resolution algorithm is integrated in a finite element method (FEM) for linear elastic thin-shell mechanics, reviewed below. First, the assumptions associated with linear thin shell theory are disscussed and governing equations in equilibrium following the formulation of Mindlin–Reissner plate theory are briefly reviewed and summarized [48]. Second, the discretized equilibrium energy equation is given within a unique surface subdivision paradigm to ensure C 1 displacement field interpolation developed in previous work [46, 47]. The foundation thin shell theory of this contact method assumes i) small thickness h in comparison with the least radius of curvature of the mid surface R, ii) large displacements and small strains so that the associated higher-order terms are neglected in comparison with first-order terms, and iii) small stress component in the direction normal to the mid surface compared with other components of stress [49]. Furthermore, the shell deformation can be described only by stretching and bending of its mid surface. For an arbitrary shell body in 3D, its undeformed geometry is defined by a mid surface of the domain Ω with boundary ¯ δΩ. ¯ Upon application of external loads, the shell body deforms characterized by a mid surface of the domain Ω with boundary δΩ. The position vector ψ¯ and ψ of a material point, associated with the curvilinear coordinates (θ1, θ2, θ3), is readily obtained for reference (undeformed) and current (deformed) configuration (see Figure 2.1), respectively: 17 x �! x ∘ �!!" Figure 2.1: Thin shell geometry in the reference (undeformed) and deformed configurations in curvilinear coordinates. ψ¯(θ1, θ2, θ3) = x¯(θ1, θ2) + θ3a¯3(θ1, θ2), (2.1a) ψ(θ1, θ2, θ3) = x(θ1, θ2) + θ3a3(θ1, θ2), (2.1b) where θ3 ∈ [− h 2 , h 2 ]; x¯ and x presents of the undeformed and deformed middle surfaces of shell parametrized by the surface curvilinear coordinates (θ1, θ2), respectively. They are directly related via the displacement field u at the middle surface under linearization: x(θ1, θ2) = x¯(θ1, θ2) + u(θ1, θ2). (2.2) Hereby, the corresponding covariant base vectors in the reference and current frame are given as follows: g¯α = ∂ψ¯ ∂θα = ∂x¯ ∂θα + θ3 ∂a¯3 ∂θα , g¯3 = ∂ψ¯ ∂θ3 = a¯3, (2.3a) gα = ∂ψ ∂θα = ∂x ∂θα + θ3 ∂a3 ∂θα , g3 = ∂ψ ∂θ3 = a3, (2.3b) 18 here α = 1, 2. The Green-Lagrange strain tensor is therefore computed as: Eij = 1 2 (gi · gj − g¯i · g¯j ). (2.4) We can expand equation Eq. (2.4) and group the first-order terms into membrane strain ϵij and bending strain κij as: Eij = ϵij + θ3κij , (2.5) where ϵij and κij are functions of displacement field u. According to the linear elastic theory, the strain energy density takes the form [46]: W(ϵij ) = 1 2 Eh 1 − ν 2 Hijklϵij ϵkl, (2.6a) W(κij ) = 1 24 Eh3 1 − ν 2 Hijklκijκkl, (2.6b) where E is Young’s modulus, ν is the Poisson ratio, h is the shell thickness and H is expressed as: Hijkl = ν ∂x¯ ∂θi · ∂x¯ ∂θj ∂x¯ ∂θk · ∂x¯ ∂θl + 1 2 (1 − ν)( ∂x¯ ∂θi · ∂x¯ ∂θk ∂x¯ ∂θj · ∂x¯ ∂θl + ∂x¯ ∂θi · ∂x¯ ∂θl ∂x¯ ∂θj · ∂x¯ ∂θk ). (2.7) In general, the shell experiences the combined effects of internal and external loads, and the corresponding energy equation accounts for the summation of these terms along with kinetic energy: Π(u) = Πint(u) + Πext(u) + Π kin(u), (2.8) and from Eq. (2.6), we can see that internal energy takes the integration form of the energy density as: Π int(u) = Z Ω W(ϵij ) + W(κij ) dΩ. (2.9) In finite element discretization, the domain Ω is partitioned into a finite element mesh in the domain ΩNEL, where NEL denotes the number of elements. The thin-shell objects are 19 vertex element Figure 2.2: A finite element mesh of a spherical shell. modeled in an unstructured triangular mesh (see Figure 2.2 for a spherical shell example). Object motion is thus entirely defined by nodal motions. Taking the weak form following the minimum principle of Eq. (2.8), one can get the first variation in the direction of virtual displacement δu in the form: ⟨DΠ(u), δu⟩ = 0. (2.10) For the present study, we assume rotation free at each node, and thus the nodal displacement is a result of the sum weighted on shape functions. By taking this into account, the discretized formulation of Eq. (2.10) yields the dynamic equation of motion as follows: Mx¨ + f int(x, t) = f ext(x, t), (2.11) colorblackwhere M represents mass matrix for the system, and f int(x, t) and f ext(x, t) are the internal and external force vector, respectively with the definitions below: f int(x, t) = ∂Πint ∂x , (2.12a) f ext(x, t) = ∂Πext ∂x . (2.12b) 20 The FEM uses subdivision-surface paradigm to create smooth shape functions on the discretization of unstructured mesh. In subdivision schemes, an initial relatively coarse mesh, often referred to as control mesh, is refined until a satisfying surface smoothness is achieved. The preceding study leverages Loop’s scheme, which approximates the newly refined mesh by allowing the old nodal positions to adjust as needed [46]. The resulting limit surfaces are globally C 2 , and details of refinement rules and the convergence analysis can be found in [46]. In the mesh region where each triangle has exactly 6 incident edges, denoted as a regular patch, the resulting limiting surface is interpolated by Box-spline shape functions [46]. Consequently, the corresponding nodal positions parameterized by local barycentric coordinates (t1, t2) are expressed as: x(t1, t2) = X 12 i=1 Ni(t1, t2)xi , (2.13) where Ni are shape functions and 12 nodes are the neighbors of the vertex (see Figure 9 in [46]). For the irregular patch, where an arbitrary number of incident edges of one or more vertex is present, the mesh needs to be subdivided in sufficient times so that each element is contained within a regular patch upon which the Box-spline shape functions are applicable. 21 Chapter 3 Contact mechanics As introduced in Chapter ??, contact-impact forces typically originate from either a regularized or discontinuous framework. This study explores two fundamental methods based on discontinuous framework for resolving contact interfaces: the non-smooth variational method and the Augmented Lagrange multiplier method. The former, aligned with the explicit time integration scheme, builds upon foundational work from [9]. The latter employs a similar contact detection strategy but formulates problems within augmented systems, making it suitable for implicit time integration schemes. These two methodologies are presented in section 3.1 and 3.2 with their discretized impact equations respectively. 3.1 Non-smooth Lagrange mechanics Non-smooth variational methods provide a robust framework in finite element analysis (FEA) for addressing contact mechanics, particularly in scenarios where traditional smooth optimization techniques may not suffice. These methods handle the inherent non-smoothness of contact conditions, such as gaps and frictional forces, effectively. The geometric contact resolution in explicit time integration focuses on a contact-enforcement method, termed as decomposition contact response (DCR), to solve for post-impact dynamics. At the same time, the non-interpenetration constraints are enforced by a geometric method separately [9]. Impenetrability is enforced through techniques like closest point 22 projections, while momentum transfer occurs via self-equilibrating impulses applied to the impacted nodes. Notably, the DCR method diverges from conventional Lagrange multiplier or penalty methods in self-equilibrating impulses. This characteristic eliminates the necessity of selecting a penalty parameter crucial for convergence, as the method is derived directly from variational nonsmooth Lagrangian mechanics. Upon a high-level summary of Lagrange mechanics in non-smooth settings, it is derived based on the foundations of smooth settings to ensure energy conservation and particularly conserved discrete symplectic form. The nonsmooth situation considers a configuration manifold Q with the tangent space TqQ consisting of all velocity vectors based at the configuration q ∈ Q, and a submanifold with boundary C ∈ Q which represent the subset of admissible configurations. Here, the boundary ∂C is the admissible contact set. Therefore, impacts happen at q ∈ ∂C at spatial location x ∈ R3 in a continuum system in R3 between two material points. A comprehensive review regarding variational nonsmooth Lagrangian mechanics can be found in [50]. We start with the action term of the system S(x,x, t ˙ c), which is the integral of Lagrangian L(x,x˙) by taking contact time tc into account: S(x,x, t ˙ c) = Z tc 0 L(x,x˙)dt + Z t tc L(x,x˙)dt. (3.1) Following the stationary principle of action, the first-order change in action δS is zero, which is given by: δS(x,x, t ˙ c) = δ( Z tc 0 L(x,x˙)dt + Z t tc L(x,x˙)dt). (3.2) Expanding the right-hand side of Eq. (3.2), we can finally reach the derivative at tc: Lδtc + ∂L ∂x˙ · δx t + c t − c = 0. (3.3) The geometrical impenetrability constraints to second-order accuracy are satisfied with a properly defined gap function g(x). Here g(x) is defined as the tetrahedron volume formed 23 by vertices participating in contact, and its variation term is: δg(x(tc)) = ∇g(x(tc)) · (δx(tc) + xδt ˙ c) = 0. (3.4) Finally, for δtc = 0, replacing Eq. (3.4) in Eq. (3.3) leads to the governing equations for contact: ∂L ∂x˙ · δx t + c t − c = 0, (3.5a) ∂L ∂x˙ t + c t − c = λ∇g(x(tc)), (3.5b) where λ ∈ R is a scalar, which means that the only admissible momentum change aligns with the direction of ∇g during the impact. Furthermore, energy conservation of contact is also encapsulated. The above equations are used to compute the momentum right after impact if the momentum immediately prior to the contact is known. Here, we denote a symbol C as a contact operator to summarize such process in the following equation: x˙ + (tc) = C(x˙ − (tc)) = x˙ − (tc) − 2 ∇g Tx˙ − (tc) ∇g TM−1∇g M−1∇g, ∇g = ∂g ∂x tc , (3.6) where x˙ − (tc) and x˙ + (tc) are pre- and post-collision velocities respectively, and M is the lumped mass matrix of the vertices in motion, which is defined and constructed by the shell FEM discretization. At this point, we apply impact mechanics conservation equations. 3.2 Augmented Lagrange multiplier method For the implicit time integration scheme, an Augmented Lagrange multiplier method is considered due to its benefits in convergence and stability. To derive the equations of motion 24 within this framework, a brief review of the Lagrange multiplier method and the penalty method is provided, as these form the foundations of the Augmented Lagrange multiplier method. Any continuum problem involving hyperelastic material can be formulated as an optimization problem. This is based on the principle that the total potential energy reaches a minimum at the solution point, as formulated in the following: Minimize Π(x), (3.7) subject to g(x) ≤ 0, (3.8) where Π represents the total energy of the system, x is the displacement vector, and g is the gap function vector respectively. In Lagrange multiplier method within finite element analysis, it introduces additional variables, Λ, known as Lagrange multipliers, to transform a constrained optimization problem into an unconstrained one, ensuring that contact conditions are satisfied accurately: Π LM (x,Λ) = Π(x) + ΛT g(x), (3.9) where Π(x) is the potential energy of the original system as described in Eq. (2.8), and together with the Kuhn-Tucker-Kraush (KKT) conditions: g(x) ≤ 0, Λ ≥ 0, g(x)Λ = 0. (3.10) At the equilibrium state, the variation form of Eq. (3.9) yields: δΠ LM(x,Λ) = δΠ(x) + δΛ T g(x) + δxT C(x) TΛ = 0, (3.11) where C(x) T δx is the matrix form of δg(x). Finally, the resulting system of equations can 25 be written as:    K CT C 0       x Λ    =    f g    , (3.12) where K is the stiffness matrix, C represents the constraint matrix, and f is the force vector respectively. Solve the above equations to obtain the displacements x and the Lagrange multipliers Λ. This system of equations is typically nonlinear due to the nature of contact constraints, necessitating the use of iterative solution techniques such as the Newton-Raphson method. The Lagrange multiplier method provides a rigorous and systematic approach to handling contact constraints in finite element analysis. It ensures that the constraints are satisfied exactly, resulting in accurate modeling of contact interactions. However, this method can be computationally intensive, particularly due to the need to calculate the contact gap function g(x) and the subsequent sorting of active and inactive nodes. Additionally, solving the resulting indefinite system of equations is generally more computationally demanding than solving positive-definite systems. Therefore, efficient numerical techniques and robust algorithms are essential for addressing these challenges in large-scale problems with complex contact conditions. The penalty method is a another prevalent technique in FEA for enforcing contact constraints between interacting bodies. This approach introduces a penalty term into the objective function, penalizing constraint violations and thereby driving the solution towards satisfying these constraints. Here, contact conditions are formulated as constraints that must be fulfilled at the interface between bodies to prevent penetration and ensure proper transmission of contact forces. The penalized objective function Π(x) is constructed by adding a penalty term to the original objective function. This term imposes a cost proportional to the violation of the contact constraints: Π P (x) = Π(x) + µ 2 g(x) T g(x). (3.13) 26 Here, µ is the penalty parameter, a large positive number controlling the severity of the penalty. Minimizing the penalized objective function Π(x) with respect to the nodal displacements x drives the solution to satisfy the contact constraints by making violations energetically unfavorable. The variation form is given by: δΠ P (x) = δΠ(x) + µC(x) T g(x) = 0. (3.14) Therefore, the finite element matrix formulation, incorporating the penalty term, can be expressed as: Kx + µgT gx = f. (3.15) Again, solving this system yields the displacements x. The nonlinearity introduced by the penalty term and contact constraints typically necessitates iterative solution techniques, such as the Newton-Raphson method. During each iteration, contact conditions are updated, and the penalty term is adjusted to ensure convergence. This iterative process continues until the solution converges, indicating that the contact constraints are satisfied and the system’s potential energy is minimized. The penalty method is straightforward and effective for managing contact constraints in FEA. It approximately satisfies the constraints by imposing high costs for violations, resulting in accurate modeling of contact interactions. However, selecting an appropriate penalty parameter µ is crucial. If µ is too small, the constraints may not be sufficiently enforced; if too large, it can cause numerical instability and ill-conditioning of the equations. To address the numerical challenges associated with the traditional Lagrange multiplier and penalty methods, the Augmented Lagrange multiplier method is considered for use in implicit time integration during contact resolution. The corresponding potential energy equation is formulated as: Π ALM (x) = Π(x) + ΛT g(x) + µ 2 g(x) T g(x) − 1 2µ λ Tλ. (3.16) 27 Here, λ represents the Lagrange multipliers associated with active constraints. The mixtures of Λ and λ ensure C 1 differentiability of Eq. (3.16). Starting with the Eq. (3.16), the system of equations to be integrated is: Mx¨ = f (x, t), (3.17) where f (x, t) = ∂ΠALM (x) ∂x , (3.18) and it is subject to admissible contact condition: g(x) ≤ 0. (3.19) First, the form of a normal Lagrangian without any augmentation is given in Eq. (3.20) after expanding the governing equation, Eq. (3.9): L(x, z) = Π(x) +X k λkgk(x)+, (3.20) where λk are the Lagrange multipliers for each individual contact event and gk(x) represent the active contact constraints where gk(x) = max (0, gk(x)). Next, the formulation of augmented Lagrangian is given by: L ALM (x, z) = F(x) + λg(x) + µ 2 g 2 (x), (3.21) with µ being the penalty parameter, which is usually a large constant. The constraints term g(x) is a sum of individual active constraint term in Eq. (5.15) as in: g(x) = X k gk(x)+. (3.22) The dual use of the Lagrangian functionals arises from a hybrid solution approach: L ALM is 28 used to solve for deformations away from contact regions, accounting for contact constraints, while L accelerates the iteration process with individual contact constraints. 29 Chapter 4 Geometric contact method in explicit time integration This chapter provides a detailed discussion on two phase contact detection and time marching methodologies within the geometric contact method for explicit time integration. A preliminary result on energy conservation is included to verify the accuracy of the contact method, demonstrating its effectiveness in preventing interpenetration. 4.1 Two phase contact detection algorithm In the global phase, the motion envelope for all triangles from tn to tn+1 are computed, and oriented coordinate aligned bounding boxes (CABB) are generated; see one CABB example in Figure 4.2b. Two triangles can potentially contact if their bounding boxes intersect. Intersection of CABB is performed by a modified BVH method with a quadtree structure from the computational geometry library (CGAL) [51]. To enhance computational efficiency, we incorporate a user customizable variable exclusion size EXCL as a threshold into the implementation, which filters out element pairs with mesh distances below. Here, mesh distance in a fully connected mesh between any two elements is defined recursively as follows: 1) mesh distance of any element to itself is zero. 2) elements that share a single edge are considered to have one mesh distance. 3) moving forward, the distance of any element to a 30 target element is one plus the minimum distance of its immediate neighbours (see Figure 4.1 for example). In practice, the bending energy required to make the two adjacent elements have contact with each other is tremendous. As a result, immediately connected triangles are removed from global contact consideration in applications since they usually do not participate in contact, where the default exclusion size is set to 2. A culling result for a hemispherical mesh in motion is demonstrated in Figure 4.2b where a default exclusion size is considered. For global phase contact detection, we leverage the high-fidelity open-source computational geometry library (CGAL) to search for intersections of the triangles corresponding to the valid CABB pairs. CGAL algorithm for intersection search utilizes a modified BVH method in a quadtree structure (details of the CGAL algorithm are described in [51]). Center to the CGAL global contact search is the function box intersection d() to compute the pairwise intersecting boxes between two sequences of iso-oriented boxes in arbitrary dimension. It utilizes segment trees, structured as balanced binary search trees, to efficiently manage and query intervals and intersecting multidimensional boxes [51]. At the end of global-phase pass, a list of element pairs, denoted as LGP , whose bounding volumes are intersected, which indicate the elements potentially in contact, is obtained for further detailed contact examination in the local phase. Putting everything together, the procedure to perform global contact check among any two mesh domains D1 and D2 is summarized in Algorithm 1, where a self contact check has D1 = D2. Algorithm 1 Global phase contact detection of time advancement tn → tn+1 1: procedure ContactBetween(D1, D2, tn, tn+1, EXCL) 2: Lbox1 ← ConstructBoundingBoxListFrom(D1, tn, tn+1) 3: Lbox2 ← ConstructBoundingBoxListFrom(D2, tn, tn+1) 4: LGP ← CGALIntersectionBetween(Lbox1 , Lbox2 , EXCL) 5: return LGP 6: end procedure 31 0 1 2 Figure 4.1: Mesh distance to a element in the center (marked as 0). (a) CABB example for one time step. (b) Cull out the non-intersected CABB pairs. The right shows the resulting CABBs saved for future use. Figure 4.2: Culling CABB pairs. In the local phase contact detection, for each pair of elements in LGP , two elementary geometrical situations are considered in the local phase: i) vertex/triangle and ii) edge/edge intersections. Vertex/triangle intersection occurs when a vertex from one triangle is about to contact with another triangle. In the example provided in Figure 4.3 (left), vertex 0 is on the verge of contacting the triangle formed by vertices 1, 2, and 3. Edge/edge intersection 32 occurs when an edge from one triangle is in contact with an edge from another triangle. In the illustrated example in Figure 4.3 (right), edge 0-1 is about to contact with edge 2-3. The contact time, tc, for each potential intersection pair is determined by solving the co-planarity condition for tc, given by (x3(tc) − x1(tc)) × (x2(tc) − x1(tc)) · (x1(tc) − x0(tc)) = 0, (4.1) Figure 4.3: Sketch of vertex/triangle (left) and edge-edge (right) intersection problems. where x0 to x3 denote nodal coordinates of vertices being examined. This equation applies to triangle-triangle and edge-edge intersections; see Figure 4.3. Eq. (4.1) is solved iteratively for tc since it is a high order polynomial of tc, with Newton method to high precision. The possible contact times, tc, are calculated for all triangle pairs having survived the global phase sorting. There are fifteen co-planarity conditions to test: six for vertex/triangle and nine for edge/edge contact (for each pair of triangles). Only admissible contact times within the time interval are retained, i.e., tn < tc ≤ tn+1, sorted and the earliest contact time is acted upon. The corresponding contact vertices are also saved in a list LLP for collision handling. This contact time computation procedure is summarized in Algorithm 2. 4.1.1 Time marching The discretized governing equations, Eq. (2.11), are advanced in time using the explicit Newmark time integration scheme, commonly used for second-order systems, written in 33 the predictor/corrector approach [52]. The nodal coordinates x(tn), velocities x˙(tn), and accelerations x¨(tn) at time tn are updated to t∗ = tn + ∆t∗ according to x(t∗) = x(tn) + ∆t∗ x˙(tn) + 1 2 ∆t 2 ∗ x¨(tn), (4.2a) x˜˙ −(t∗) = x˙(tn) + (1 − γ)∆t∗x¨(tn), (4.2b) x˜˙ +(t∗) = C x˜˙ −(t∗) , (4.2c) x¨(t∗) = f ext(x(t∗)) − ∂Πint ∂x [x(t∗)] m , (4.2d) x˙(t∗) = x˜˙ +(t∗) + γ∆t∗x¨(t∗), (4.2e) where γ is the Newmark time integration parameter for temporal acceleration interpolation, 0 ≤ γ ≤ 1, which we usually set to 0.5 as it is equivalent to explicit central difference scheme. Here, we consider undamped system. The time step ∆t employed in simulation is usually set to the stable time step, which is the maximum Courant stability time-step in the explicit time scheme, given by ∆t = 2 ωGmax . Here, ω G max is the maximum eigen frequency of the entire finite element assembly, estimated from the smallest finite element in the mesh. If there is contact at time tc between t = tn and tn+1 = tn +∆t, where ∆t is the time step in simulation, then ∆t∗ = tc − tn and t∗ = tc. Otherwise, one sets t∗ = tn+1 and ∆t∗ = ∆t when no contact is predicted in the time step interval. The pre-contact intermediate velocity x˜˙ −(t∗) of elements involved in contact are generally modified accordingly by the contact operator C to the post-contact value x˜˙ +(t∗), as described in the previous section. After the mesh configuration is updated to t∗, we check again for more instances of contact and repeat the process until no further contact is detected within t ≤ tn+1. At this point, we complete the time step to bring the mesh exactly to tn+1, and repeat until the simulation reaches its end time. The corresponding local phase contact detection procedure incorporated with recursive time scheme is summarized in Algorithm 3. 34 Algorithm 2 Contact time t∗ computation 1: procedure ComputeContactTime(LGP ) 2: t∗ ← Double Max 3: LLP ← V ector() 4: for each (s1, s2) in LGP do 5: for each t-t and e-e couple do 6: Compute contact time t ′ ∗ ▷ Eq. (4.1) 7: if t ′ c < t∗ − ϵ then 8: t∗ ← t ′ ∗ 9: LLP clear 10: LLP add (x0, x1, x2, x3) 11: else if |t ′ ∗ − t∗| ≤ ϵ then 12: LLP add (x0, x1, x2, x3) 13: end if 14: end for 15: end for 16: return t∗, LLP 17: end procedure Algorithm 3 Local phase contact detection of time advancement tn → tn+1 1: procedure LocalContact(D, LGP , tn, tn+1) 2: tinit ← tn 3: while True do 4: t∗, LLP ← ComputeContactT ime(LGP ) 5: if not tinit < t∗ < tn+1 6: break 7: Predict D from tinit to t∗ ▷ Eq. (4.2) (a) & (b) 8: for each (x0, x1, x2, x3) ∈ LLP do 9: x˙ + (t∗) ← C(x˙ − (t∗) ▷ Eq. (4.2) (c) 10: end for 11: Correct D from tinit to t∗ ▷ Eq. (4.2) (d) & (e) 12: tinit ← t∗ 13: t∗ ← Double Max 14: end while 15: Time advance D from tinit to tn+1 without contact treatment 16: end procedure In highly geometrically nonlinear cases, contacts may occur at various time scales, with some sequential contacts having much smaller time intervals between them than others. Consequently, when performing numerical integration over the overall time step ∆t composed of small and large time steps, the susceptibility of the mesh to interpenetrations is increased. 35 To address this issue, we devised a recursively time step subdivision scheme to account for the various contact time-scale in some applications. The idea is to examine interpenetrations among the triangles that survived the global phase at the end of each full time step. If interpenetration is found, as there should be none, the time step is reduced by half and the algorithm restarts at the beginning of the original time step and proceeds with the reduced time step. Otherwise, if the contact events in time span ∆t reach to a maximum threshold, the program follows the same procedure. This maximum threshold is a user defined parameter, which is by default set to 1000. Here we assume that the intersected CABBs will remain unchanged throughout the entire time step, therefore we perform the global contact check only once for each full time step. The procedure repeats itself until no more interpenetration is found at the end of the full time step (see Figure 4.5 for demonstration). 4.1.2 Parallel implementation The task of contact parallelization in distributed-memory computers is to break the global contact search into local search, which then becomes a self-contained problem that can be resolved on each individual process. In this context, we aim to establish a separate parallelization strategy characterized by uniform actual computational load distribution, in contrast to the conventional domain decomposition approach commonly employed in thin shell mechanics without considering contact. Our contact algorithm is parallelized for simulations using the standard message passing interface (MPI). The FEM shell solver employs a domain decomposition strategy with load balancing (Zoltan) and an optimized point-to-point communication pattern [53] such that each processor Pi is assigned with a local mesh Di (see Figure 4.4a as an example). Without considering contact, each rank solely retains its local mesh, and communication between processors only occurs among the elements at the boundaries of processors. 36 �! �" �# (a) FEM shell solver mesh partition. (b) Contact search domain Dc 0 , (c) Dc 1 , (d) Dc 2 . Figure 4.4: 3 MPI task for FEM and contact search partition. In contact parallelization, each MPI rank maintains a contact search domain Dc in addition to its local mesh throughout the simulation. This Dc defines the region for which the rank is responsible for conducting inter-process contact checks against its own FEM mesh domain. It is intuitive that each processor should be responsible for the triangle intersection tests within its own mesh domain. For inter-processor contact detection where two elements are from different processors, due to the symmetry nature of contact, it is clearly a waste for requiring two processes to resolve contact detection for the same element pair. To eliminate the duplicate pass, we devise a separate communication map for the parallel global search phase in contact detection to check for the same pair only once in either processor. For example, in a three processes application, three inter-processor contact checks are required: 37 rank 0 needs to check contact against rank 1; rank 1 is responsible for checking contact against rank 2; rank 2 needs to check against rank 0. Overall, the contact search domain is D0∪D1, D1∪D2 and D2∪D0 for processor 0, 1, and 2 respectively (see Figure 4.4a (b)-(d) for demonstration). In general, the communication takes a round-robing pattern to distribute the load of contact search evenly among processes, according to the following equation: D c i = (∪ min(j+k i ,n−1) j=i Dj ) ∪ (∪ i+k i−n j=0 Dj ), (4.3) where Dc i is the contact search domain for processor i (0-indexed), k i is the number of cross-processor contact checks that this processor needs to perform, n is the number of total processors. The total number of cross-processor contact checks is n(n−1) 2 , therefore each processor needs to consider at least n−1 2 contact checks, and additional processors need to take the remainders as following: k i = n − 1 2 + (i < mod(n(n − 1) 2 , n)) (4.4) which ensures no duplicate contact check on two processors and a balanced load distribution since the difference in k for each MPI rank cannot exceed one. To make the global contact problem self-contained, it is essential to account for multiple interprocessor communications and synchronizations of nodal kinematics during both phases of contact detection in the process. Before the global phase, the contact search domain is synchronized regarding nodal kinematics. Then it performs the CABB intersection tests among elements in Dc . The resulting list of intersected CABB pairs LGP is then obtained for a detailed contact check in the local phase. In the local phase parallelization, the local mesh D is predicted to the global earliest admissible contact time via Newmark predictor in Eq. (4.2) (a) & (b) for all ranks if any contact in the local phase was detected. The nodes that are involved in this contact event must be synchronized w.r.t. their nodal kinematics again before the repulsive force computation in 38 Eq. (3.6). After applying the contact force and updating the post-impact velocities, another synchronization occurs to send the updated post-impact velocities back to the contact triangles from other processors (synchronization regarding LLP ). To proceed with subsequent local contact detections, it is necessary to factor in synchronization related to the LGP . The abovementioned parallelization algorithm is summarized in Figure 4.5 for an arbitrary processor Pi with FEM mesh domain Di in a full time step integration: Start Time Step End Time Step Broad Phase Narrow Phase Narrow Phase Narrow Phase Max Contact Events Penetration Time Integrator yes yes halve Δ� 1 3 3 Compute Contact Time Predict to �∗ 2 a Resolve Collision 2 b Correct to �∗ 2 b Narrow Phase 1. Gather & sync �" 3. Sync nodal kinematics in �#$ 2 a. Allreduce to get global min �" 2 b. Sync nodal kinematics in �%$ 3 repeat to completion of full Δ� Figure 4.5: Contact detection in recursive time advancement scheme in parallel. The algorithm’s efficiency in parallelization and its ability to resolve complex contact scenarios are demonstrated through simulation examples in Chapter 7. These results confirm the algorithm’s success in generating realistic shell deformation configurations while preventing interpenetration. 4.2 Preliminary results In this section, preliminary results to demonstrate the correctness of the contact algorithm in explicit time integration are showcased in two simulations: 1)contact between two hemispherical shells in head-on and rim-rim alignments and 2) impact between a spherical shell 39 and a rigid wall. The first example simulates a process in which two hemispherical shell objects make contact and then subsequently separate. The latter example is aimed to demonstrate the energy conservation properties of the contact algorithm. As discussed in [9], a notable drawback of back-projection, necessary for geometric non-interpenetration conditions, is the need to adjust the artificial virtual work associated with it. However, with the geometric contact resolution presented in this chapter, this is no longer a concern, as the enforcement of geometric constraints is inherently embedded in the algorithm itself within each time step. This major advantage in automatic energy conservation is demonstrated through a preliminary example of a spherical shell impacting a rigid wall, where the energy is nearly perfectly conserved. 4.2.1 Hemispherical shell collision Two simulations are considered in this example: head-on collision, and rim-rim collision (see Figure 4.6 left and right respectively). In both simulations, we use two identical hemispherical shells, marked as red and blue in Figure 4.6, with radius 5 cm. At the beginning of the simulation, the two hemispheres are assigned initial velocities to move them toward each other: the red and blue hemispheres for both cases are initialized with velocities of 10 m/s, and 5 m/s downward and upward, respectively. The corresponding material is from a family of Nylon type B, with the following properties: density 1.15 × 103 kg/m3 ; Young’s modulus 0.83 × 109 kg/m s2 ; Poisson ratio 0.33; thickness 6.86 × 10−5m. The finite-element mesh of each hemispherical shell has 6,456 triangular elements. 40 Figure 4.6: Two hemispheres simulations at t = 0 in head-on collisions (left) and rimcollisions (right). Figure 4.7a and 4.7b present the results from the onset of contact to the separation of two hemispheres for head-n and rim-rim respectively. The head-on simulation results demonstrate typical blunt soft materials deformation behaviors due to elastic collision. The rim-rim simulation sees a longer time of contact. Notably, in this case, the objects have rotated ∼ 3 4 π about the contacting edges before separation due to the angular momentum imposed upon contact. 41 0.004s 0.006s 0.007s 0.009s (a) Head-on collision. 0.004s 0.008s 0.011s 0.015s 0.019s 0.020s (b) Rim-rim collision. Figure 4.7: Simulation results over time. Parachute deployment Unlike the other two simulations, the parachute deformation is driven by a specified deployment protocol to simulate the parachute behavior with payload. The protocol includes a few parts unique to parachute application, such as cable lines and reefing lines, to control the parachute’s open area. The simulations use shell mechanics for the fabric, cable mechanics for the suspension and reefing lines, and self-contact algorithms for contact detection and handling. The canopy has been fitted with a reefing line that can be expanded, contracted, and cut as needed, see Figure 7.12. We assume there is a single riser attached to a fixed point (where the payload would be) with 24 suspension lines extending from the riser to the rim of the canopy. The riser line has the properties of the aggregated 24 suspension lines. The finite-element model of the parachute has 7056 triangular shell 42 elements, 1511 finite cable segments in the shell, and 506 external finite cable segments in the suspension lines, riser and reefing lines. The material properties of the fabric are the same as hemispheres (see section 4.2.1). The suspension lines, reefing line, and riser line have linear density 0.0175 kg/m, 0.0248 kg/m, and 0.42; kg/m, respectively, and stiffness is 6.31 × 105 kg m/s 2 , 7.81 × 105 kg m/s 2 and 1.51 × 107 kg m/s 2 , respectively. We also use a special material type for the cables along the canopy skirts, which has linear density 0.0301 kg/m and stiffness 5.52 × 104 kg m/s 2 . Figure 4.8: Main parameters of drogue parachute constructed geometry. 4.2.2 Deployment protocol The strategy devised to partially collapse the original constructed geometry follows a loading sequence described schematically in Figure 7.13, and consists of the following processes: • A prescribed load pressure, pi , is applied on the inner surface of the chosen panels of the canopy. In the present simulations we apply a constant load either on: a) the end disk, b) the end disk and next band of panels, or (c) the end disk and subsequent two bands of panels of the parachute. These are denoted by 1-, 2- or 3-band loading, respectively, throughout this section. • The initial length of the reefing line is equal to the perimeter of the canopy rim, of diameter D = 4.48 m, in its constructed geometry specification. This length is substantially larger than the actual length of the reefing line in the first and second stages of reefing which have diameters of D1 = 0.656D and D2 = 0.9D, respectively. 43 To achieve a geometrical configuration with the actual diameter of the reefing line we proceed by reducing its diameter, in four steps from t = 0 and up to t ∗ , by 90% each step (for a total of 0.9 4 ≈ 0.656). In between steps, we allow the parachute to deform dynamically to accommodate this change in the reefing line diameter. Each shortening step of the reefing line pulls the skirt inwards (radially). Therefore, we freeze the canopy by setting all velocities to zero shortly after shrinking the reefing line to prevent the parachute from collapsing owing to the inertia induced by the abrupt shortening of the reefing line. We have experimented with two protocols during this stage. In one case we impose the uniform loading pressure pi during the whole duration of the simulation while in a second strategy we imposed an over-pressure po during the first step of shortening from 0 < t ≤ to; this is referred to as the over-pressured startup. • After the reefing line has achieved its nominal first-stage reefing diameter D1 at t = t∗, the canopy is allowed to deform dynamically up to t = t∗ + t1. The stresses and loads on the suspension lines are allowed to follow the normal dynamics without any changes to the geometry or topology of the parachute. • Immediately after t = t∗+t1, the reefing line is extended instantaneously to the secondstage reefing diameter D2, and the simulation is carried out up to t = t∗+t1+t2 without further modification. • At t = t∗ + t1 + t2, the reefing line is removed and the parachute is allowed to expand to its fully deployed configuration. The values of the applied load and timings can be varied. The inflation pressure pi used here is approximately the transmural pressure that a parachute deployed at 4000 m altitude. We conducted our experiments in three overloading settings for po = 10, 20, and 40 kPa, respectively. The other parameters we used are pi = 4 kPa, t∗ = 0.0141 s, t1 = 0.00703 s, and t2 = 0.00703 s. Two overpressure loadings were used, one applied over t0 = 0.00351 s 44 and one over the full initial shrinking of the reefing line and including t1, i.e., t0 = 0.0211 s. These are denoted as burst and full overloading, respectively. Figure 4.9: General inflation pressure timeline evolution. 4.2.3 Simulation results The simulations show typical behavior where the pressurized canopy induces axial elongation and subsequent partial collapse of the parachute, see figures 7.14, 7.15 and 7.16. Shortly after simulation startup, there is global self-contact of the structure and some panel folding, which needs to be resolved correctly for the parachute to be in a physical state, free of interpenetration. The suspension lines display their own dynamics, which are not fully displayed in the figures, but evidently couple with the geometry displayed by the canopy. The figures show that the chosen load protocol has a large influence on the state of the parachute, which is to be expected. Evidently, different load protocols can be used to generate even more varied geometries, and the possibilities are very large. Here, we focus on the loading protocol sketched in Figure 7.13. To aid in the description of the results, we will refer to each simulation by the row-column ordering in the figures, e.g., R2-C1 denotes the second row and the first column. The first row and column correspond to the top row and leftmost column, respectively. 45 Figure 4.10: Results for 1-band loading scenario. From top to bottom: at the end of stage 1, at the end of stage 2, after disreefing, where t = t∗ + t1, t∗ + t1 + t2, t∗ + t1 + t2 + 0.00703 s, respectively. From left to right: first no overloading; then full overloading and burst overloading at po = 10, 20, and 40 kPa, respectively. Figure 4.11: Results for 2-band loading scenario; arrangement as in Figure 7.14. 46 Figure 4.12: Results for 3-band loading scenario; arrangement as in Figure 7.14. Loading only the back disk (1-band) of the parachute generates a more elongated geometry, as can be seen by comparing the first row of all three figures 7.14, 7.15 and 7.16. Loading 3 bands tend to produce more hemispherical geometries, while loading 2 bands is more cylindrical in nature, although there is variability depending on the overloading. The largest overloading pressures result, invariably, in the most collapsed geometries, see R3-C6 results. Some patterns are interesting, like 1-band R3-C4 and C6, and 3-band R2-C6 and R3-C6, showing a conically expanded skirt and overinflated disk (first band). Figures 7.17, 7.18 and 7.19 show histograms of von Mises stresses at the corresponding times and for all cases discussed before. The stress is essentially bounded between 0.1 MPa and 104 MPa for all cases. The distribution is more or less Gaussian with a mode about 100 MPa and a few exceptions that show bimodal distributions, such as in 1-band loading R2-C5 and C7 and 3-band loading R1-C1. The distribution of stresses is informative but obviously not as specific as that coming from an fluid-solid interface simulation. Furthermore, it is well known that the details of the constitutive law of the fabric have a large influence on the calculated stresses, but not as much on the geometry of the parachutes since that is 47 regulated by the boundary conditions with the flow. Figure 4.13: Histograms of von Mises stress (σv) for 1-band loading in linear-log coordinates: lower and higher values of σv are 0.1 MPa and 104 MPa, respectively; arrangement as in Figure 7.14. 48 Figure 4.14: Histograms of von Mises stress (σv) for 2-band loading in linear-log coordinates; arrangement as in Figure 7.14. Figure 4.15: Histograms of von Mises stress (σv) for 3-band loading in linear-log coordinates; arrangement as in Figure 7.14. 49 4.2.4 Spherical shell bounces off a rigid wall The study of the impact between an elastic object and a rigid wall has a long history of research and has made significant contributions to the field of mechanics. Numerous studies in the past have investigated this classical problem, providing valuable insights into the fundamental physics of impact and collision. Some critical aspects explored in these studies include i) exploring the restitution properties w.r.t. different incident angles and velocities [54, 55], ii) evaluating the extent of validity of Hertz theory by comparing contact pressure and forces under various conditions [56, 57], and iii) demonstrating the efficacy of proposed contact algorithms through conservation of energy analyses [58, 59, 60, 61]. This study focuses on the energy conservation behaviors of a gravity-driven thin spherical shell bouncing off a rigid wall for various material stiffness. We set a spherical shell 10 cm from a rigid wall so that it initially falls toward the wall under the gravitational acceleration of 9.81 m/s 2 (see Figure 4.16). For proper contact handling, both the spherical shell and the wall are modeled in FE mesh of similar mesh size, with 6,362 and 11,861 triangular elements, respectively. The spherical shell has outer radius 10 cm, thickness 1 mm, density 4.0 ×103 kg/m3 , Poisson ratio 0.33, and Young’s modulus varying between 2.1×106kg/m s2 , 2.1×107kg/m s2 , and 2.1×108kg/m s2 . The wall mesh is of a squared-shape with side length 40 cm. 50 g Figure 4.16: Spherical shell-wall collision simulation: front view (left); top view (right). To appropriately simulate the rigid wall, the corresponding wall mesh is considered to have infinite nodal mass and zero velocities, which results in zero contribution to the computation of reaction momentum and velocities. In practice, we only update the post-impact velocities of the triangles from the spherical shell mesh whenever there is a sphere-wall contact, leaving the wall mesh to remain stationary. As expected, our simulation exhibits the total energy conservation for multiple bounces (kinetic energy, gravitational potential energy, and elastic internal energy) (see Figure 4.17a-4.17c). In all cases, the system begins with an initial total energy of 0.981 J, solely attributed to potential energy. As the spherical shell moves, this potential energy progressively transforms into kinetic energy until the moment of impact with the wall. Upon collision, the spherical shell experiences an increase in internal elastic and kinetic energy due to significant deformation caused by the impact. The potential energy ultimately asymptotes because there is an increasing transfer of kinetic energy to the mesh since there is no source of dissipation in this elastic model. 51 0.0 0.2 0.4 0.6 0.8 1.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 energy (J) kinetic energy potential energy internal energy total energy (a) E = 2.1 × 106 Pa. 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 energy (J) (b) E = 2.1 × 107 Pa. 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 energy (J) (c) E = 2.1 × 108 Pa. Figure 4.17: Energy partition over time during impact for various material stiffness. (see in color print) Materials with higher stiffness demand additional time to attain a stable equilibrium in terms of potential, kinetic, and internal energy. For instance, the least stiff material, as depicted in Figure 4.17a, quickly splashes onto the wall immediately upon impact, resulting in the shortest time required for stability. On the other hand, materials with medium stiffness (Figure 4.17b), achieve energy stability after undergoing approximately 4-5 cycles of rebound. Furthermore, the stiffest material, shown in Figure 4.17c, takes around 30 cycles to approach an asymptotic state. Here, we have also included a histogram of nodal contact impulse throughout the simulation, denoted as m(x˜˙ +(t∗) − x˜˙ −(t∗)), for the sphere vertices 52 contacting with the wall for varying material stiffness in Figure 4.18. All cases exhibit very similar behavoir despite the elastic modulus changing by two orders of magnitude. The distributions are approximately log-normal and peak at about 1.8 N s. The similarity of the statistics if a consequence of the nature of the impulse, which is related to the mass and velocity of the mesh, and less to the elastic properties of the material. 0 2 4 6 impulse (Ns) ×10−4 (a) E = 2.1 × 106 Pa. 0 2 4 6 impulse (Ns) ×10−4 (b) E = 2.1 × 107 Pa. 0 2 4 6 impulse (Ns) ×10−4 (c) E = 2.1 × 108 Pa. Figure 4.18: Impulse histogram for different materials. We further investigate the impact of different choices of the algorithmic damping parameter, γ, on the energy exchange behavior. In addition to the case where γ = 0.5, we consider two other scenarios where γ takes the values of 1 and 10, respectively. The corresponding results are depicted in Figure 4.19 for varying material stiffness. In damped cases 53 where γ > 0.5, the asymptotic state potential energy closely resembles that of undamped systems (γ = 0.5). The introduction of artificial damping dissipates kinetic energy, which would otherwise be exchanged with internal energy in undamped systems at high frequencies. Consequently, the internal energy converges to significantly lower asymptotic values compared to undamped cases and the total energy is merely the sum of potential energy and internal energy at the asymptotic state. Further, The values of γ in damped cases influence the convergence speed to the asymptotic state, with higher values of γ introducing more damping to the kinetic energy initially, leading to quicker convergence. 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy γ=0.5 γ=1 γ=10 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.19: Energy conservation for various γ of material stiffness E = 2.1 × 106 Pa. (see in color print) 54 0.0 0.5 1.0 1.5 2.0 2.5 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy γ=0.5 γ=1 γ=10 0.0 0.5 1.0 1.5 2.0 2.5 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.0 0.5 1.0 1.5 2.0 2.5 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.0 0.5 1.0 1.5 2.0 2.5 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.20: Energy conservation for various γ of material stiffness E = 2.1 × 107 Pa. (see in color print) 55 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy γ=0.5 γ=1 γ=10 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.21: Energy conservation for various γ of material stiffness E = 2.1 × 108 Pa. (see in color print) 56 0.00 0.05 0.10 0.15 0.20 time(s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy coarse medium fine fine 0.00 0.05 0.10 0.15 0.20 time(s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.00 0.05 0.10 0.15 0.20 time(s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.00 0.05 0.10 0.15 0.20 time(s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.22: Energy conservation in one bounce off for various mesh resolutions. 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 energy (J) (a) Coarse mesh (2329 spherical shell elements). 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 energy (J) (b) Fine mesh (36,649 spherical shell elements). Figure 4.23: Energy conservation in one bounce off for various mesh resolutions. 57 Figure 4.24: Energy conservation for a medium fine mesh in one bounce off. Energy conservation remains well-preserved throughout all cases, as the variation in total energy for each case is kept within 0.8%. This small variation is simply a consequence of truncation error in the spatial and temporal discretization. To conclusively establish the preservation of energy, we further examined this conservation attribute across a range of mesh resolutions under different mesh configurations during a single bouncing event.Two mesh configurations are considered in this study: one with uniform meshes and the other with skewed meshes. The uniform meshes maintain consistent element shapes, while the skewed meshes feature triangular elements in the spherical shell with significantly varied aspect ratios, leading to reduced uniformity. For both configurations, three mesh resolutions are analyzed with respect to energy conservation: coarse, medium, and fine (see Figure 4.25a for demonstration). In the uniform mesh setting, the coarse, medium, and fine meshes (which is the one in Figure 4.16) consist of 192, 2,329, and 6,362 elements, respectively. In the skewed mesh setting, these resolutions consist of 458, 1,080, and 2,904 elements, respectively. For 58 both configurations, the wall mesh uses the same as in Figure 4.16. This mesh convergence test is conducted on a material with a stiffness of 2.1 × 106 Pa. As expected, for each mesh configuration, the finer the mesh, the closer it gets to exact values across all energy partitions as well as total energy (see Figure 4.25b). 59 Uniform mesh Skew mesh (a) Mesh configuration for uniform and skew meshes. From left to right: coarse, medium and fine mesh. 0.00 0.05 0.10 0.15 0.20 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy coarse (u) medium (u) fine (u) coarse (s) medium (s) fine (s) 0.00 0.05 0.10 0.15 0.20 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.00 0.05 0.10 0.15 0.20 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.00 0.05 0.10 0.15 0.20 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy (b) Energy conservation during a single bounce for different mesh resolutions for both uniform and skew mesh configurations. Red curves labeled (u) correspond to the uniform mesh setting, while blue curves labeled (s) correspond to the skewed mesh setting. 60 To delve deeper into the robustness of the contact search strategy, we extend our analysis to evaluate the conservation attribute across different mesh configurations. Here, we utilize the notation ’a-a’ to differentiate between the treatment of spherical and rigid wall shell meshes: the first letter designates the selection of the spherical shell mesh, while the latter letter denotes the choice of the rigid wall mesh. The letter can be either ’r’ or ’c’. ’r’ represents the reference mesh used in Figure 4.17. ’c’ indicates the uniform medium mesh utilized to assess mesh resolution convergence in Figure 4.25a. Therefore, ’r-r’ denotes the mesh configuration that utilizes the reference mesh for both the sphere and rigid wall, mirroring the mesh setup employed in Figure 4.17. This case serves as the baseline configuration and is later used for comparison with other mesh configurations. In this context, we explore two other mesh configurations: ’r-c’ and ’c-r’ (see Figure 4.26a and 4.26b respectively). Additionally, we consider a skew mesh, wherein the spherical shell is relatively coarser than the rigid wall, yet all mesh elements exhibit significantly higher skewness and less uniformity. The corresponding results are displayed in Figure 4.27, 4.28 and 4.29 for materials of different stiffness. In all cases, energy conservation is maintained. The outcomes of the diverse mesh investigations are predominantly shaped by how the sphere is meshed. The structure of the rigid wall is unlikely to undergo significant changes with a coarser mesh compared with a finer mesh; however, there may be slight changes in the contacts. Specifically, in the ’r-c’ case, consistent behavior akin to ’r-r’ is observed across all material stiffness levels. This similarity arises from the fact that both cases employ identical meshes for the spherical shell, and the rigid wall does not engage in momentum exchange with the sphere materials. Conversely, the ’c-r’ and skew cases swiftly approach their asymptotic states, primarily attributed to the relatively coarse meshes used for the spherical shell compared to the other cases. The coarse meshes are insufficiently fine to display behavior within the asymptotic convergence range, as the values’ rates are not very close to their exact values. Notably, the skewed mesh case exhibits significantly lower initial potential energy and total energy compared to other cases. This is due to its smaller surface area of 0.1168 m2 , in contrast to the reference spherical 61 mesh (uniform fine mesh in Figure 4.25a) with a surface area of 0.1254 m2 , and the coarse spherical mesh (uniform medium mesh in Figure 4.25a) with a surface area of 0.1246 m2 . Since all meshes have the same material density and thickness, the reduced surface area in the skewed mesh results in a lower overall mass and, consequently, lower potential energy. As demonstrated in Figure 4.25b, increasing the resolution of this skewed mesh ultimately brings it closer to the exact values. (a) r-c (top view). (b) c-r (top views). (c) skew mesh. From left to right: top view, the impact location on the rigid wall marked with a blue circle, and side view depicting the sphere impacting on the wall respectively. Figure 4.26: Mesh configurations. 62 0.0 0.2 0.4 0.6 0.8 1.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy r-r r-c c-r skew 0.0 0.2 0.4 0.6 0.8 1.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.0 0.2 0.4 0.6 0.8 1.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.0 0.2 0.4 0.6 0.8 1.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.27: Energy conservation for various mesh configuration of material stiffness E = 2.1 × 106 Pa. (see in color print) 63 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy r-r r-c c-r skew 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0.0 0.5 1.0 1.5 2.0 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.28: Energy conservation for various mesh configuration of material stiffness E = 2.1 × 107 Pa. (see in color print) 64 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 kinetic energy r-r r-c c-r skew 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 potential energy 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 internal energy 0 2 4 6 time (s) 0.0 0.2 0.4 0.6 0.8 1.0 total energy Figure 4.29: Energy conservation for various mesh configuration of material stiffness E = 2.1 × 108 Pa. (see in color print) 65 Chapter 5 Geometric contact method in implicit time integration In the previous chapter, a geometrically accurate contact algorithm incorporating an explicit time scheme, has been proven to deliver physically accurate solutions that preserve energy conservation. This is achieved by significantly reducing contact search efforts through a coupled global-local phase searching strategy followed by exact geometric constraint enforcements. However, several simulation difficulties have been observed in the presented study: (i) Recursively computing contact time is inefficient in parallel environments because the earliest contact time must be broadcast to all processors, requiring synchronization (serialization). (ii) The contact search can be highly unbalanced during parallelization due to its dependence on the initial mesh domain decomposition. (iii) The inherent nonlinearity in the contact problems makes the detailed contact behaviors highly sensitive to the time step size [62]. Some efforts have been made to address the issue mentioned in point (ii), where a separate mesh domain is proposed for handling contact detection [63]. However, there is a clear tradeoff between achieving overall load balancing and minimizing communication overheads, and the frequency of contact detection decomposition remains a problem-dependent parameter. Given all the trade-offs and limitations, an implicit approach, though traditional, offers a classical resolution in contact mechanics to avoid many recursive contact checks. This 66 implicit approach is based on the foundation of the Augmented Lagrange multiplier method introduced in section 3.2, which involves augmenting the equations of motion and solving them simultaneously. The implicit method does not require as intensive consideration of material contact within a single time step as the explicit approach does due to geometric non-linearity, thereby significantly reducing the effort needed for recursive contact searches. In the implicit time scheme, the contact detection follows the similar two phase contact detection procedures in the explicit time scheme described in section 4.1, but they are fundamentally different in the local phase. Here, the algorithm first detects triangle intersections to identify the exact contacting triangle pairs whilst in an explicit time scheme, such interpenetration is not allowed. Then, a signed gap function is computed to apply the corresponding constraint contribution to the augmented system according to the formulation in section 3. After establishing the constraint system, a GMRES approach is employed to solve the linearized equations during Newton iterations. This chapter begins with a discussion on the calculation of the gap function, which is crucial for determining the contacting triangle pairs in the local phase. It then gives a comprehensive analysis of the contact algorithm, detailing the contact detection methods and the procedures used to solve large augmented systems. Following this, the chapter describes time-marching procedures for solving the augmented constraint system. Finally, the chapter concludes with a simple example demonstrating the validity of the contact algorithm in implicit time integration schemes. 5.1 Signed gap function As discussed in the formulation of the augmented Lagrange Multiplier method, the scalar gap function g(x) is crucial for implementing the implicit constraint. Similar to the gap function evaluation in the explicit time scheme, the same procedures are applied here, considering both vertex-triangle and edge-edge contact in the local phase. For instance, consider two triangles A and B with vertex coordinates xA and xB, respectively. The tetrahedral volume 67 defined by a given 4-tuple of vertices T(x j B ; xA), where j = 1, 2, 3, is used to measure the gap function value when the vertices of triangle A are ordered. This operator is applied to all 4-tuples of vertices in a given pair of triangles, as each triangle is tested against multiple vertex events. There are a total of six possible tetrahedral volumes for vertex-triangle pairs and nine possible tetrahedral volumes for edge-edge pairs. The tetrahedral volume (without the factor 1/6) is defined by: T(x 0 ; x 1 , x2 , x3 ) = v a · v b × v c = ϵilmv a i v b l v c m, (5.1) where v a = x 1 − x 0 , vb = x 2 − x 0 , vc = x 3 − x 0 , (5.2) and ϵilm is Levi-Civita tensor. At some time t n , assume there is no interpenetration, so the gap function g n < 0 by definition in section 3.2. If contact occurs between t n and t n+1, the triangles will intersect at t n+1 and the sign of g n+1 will differ from that of g n . In the illustrated example shown in Figure 5.1, during the vertex-triangle intersection, the bottom vertex from A intersects with B from t n to t n+1, resulting in a sign change in its corresponding gap function value (the gap function at t n+1 is highlighted in yellow dashed lines.). Similarly, in the edge-edge case, the gap function changes its sign from negative to positive as the bottom edge from A intersects B. 68 n g n < 0 g n+1 > 0 x n x n+1 TB TA n g n < 0 g n+1 > 0 Figure 5.1: vertex-triangle (left) and edge-edge (right) intersections. The gap function in both cases are modeled as intersected tetrahedron volumes, whose sides are highlighted in yellow dashed lines, formed by intersected vertices. Furthermore, the gradient of this quantity is ∇x1T = v b×v c , ∇x2T = v c×v a , ∇x3T = v a×v b , ∇x0T = −∇x1T−∇x2T−∇x3T, (5.3) and the Hessian is ∇x0∇x0T = ∇x1∇x1T = ∇x2∇x2T = ∇x3∇x3T = 0, (5.4) ∇x1∇x2T = −∇x2∇x1T = v c×, (5.5) ∇x3∇x1T = −∇x1∇x3T = v b×, (5.6) ∇x2∇x3T = −∇x3∇x2T = v a×, (5.7) ∇x0∇x1T = −∇x1∇x0T = (v c − v b )×, (5.8) 69 ∇x0∇x2T = −∇x2∇x0T = (v a − v c )×, (5.9) ∇x0∇x3T = −∇x3∇x0T = (v b − v a ) × . (5.10) Once the actual contact is determined through contact detection, the augmented formulation is then to sum of all contacting triangle-triangle pairs identified by the contact search algorithm. The implementation challenge stems from the high dimensionality of the gap function g(x), which has 3N arguments, where N represents the number of vertices in the FEM mesh. This function resides in a high-dimensional space without a clear geometric interpretation. 5.2 Time marching A generalized Newmark-β method is employed to integrate the implicit time scheme. For the unconstrained problem where contact is not considered, the time integration is given by x n+1 = x n + ∆t vn + ∆t 2 2 (1 − 2β)a n + 2βan+1 , (5.11) v n+1 = v n + ∆t (1 − γ)a n + γan+1 , (5.12) where v = x˙ and a = x¨, with β and γ as constants. The acceleration a(x, t) is defined by a(x, t) = f ext(x, t) − f int(x, t) m . (5.13) Here, n represents the discrete time steps with t n+1 = t n + ∆t, where ∆t is the time step size. The implicit midpoint method corresponds to β = 1/4 and γ = 1/2, while the most dissipative method has β = 1 and γ = 3/2. These parameters are related to the amplification parameter ρ through β = 1 (1 + ρ) 2 , γ = (3 − ρ) 2(1 + ρ) . (5.14) For stability, it may be preferable to use an implicit scheme (β ̸= 0). The nonlinear system then becomes: f (x) = x − x ∗ − β∆t 2 a(x, tn+1) = 0, (5.15) where x = x n+1. The predictor is given by x ∗ = x n + ∆t vn + ∆t 2 2 (1 − 2β)a n . (5.16) The velocity equation in Eq. (5.12) is decoupled and updated after solving for x. 5.2.1 Master iteration with augmented Lagrange multipliers The deformation of the shell is primarily due to the displacement of vertices that are not in contact, even when substantial contact occurs. Classical Lagrange multipliers produce algorithms that converge less effectively than unconstrained formulations. However, enforcing the contact constraints is still necessary. Therefore, all constraints are condensed into a single one by summing them, as shown in Eq. (3.22), and an augmented Lagrange multiplier technique is used. The unconstrained minimization of the augmented Lagrange functional in Eq. (3.21) and Eq. (3.17) leads to: fi(x) + µ ν g(x) + λ ν µν X k ∇xi gk(x)+ = 0, (5.17) g(x) + λ ν µν = 0, (5.18) where there is only a single Lagrange multiplier λ ν at each outer iteration. The latter equation can be explicitly simplified by noting that contact requires z = 0, thus leading to: ¯fi(x) = fi(x) + ˜fi(x) = 0, (5.19) 71 where ˜fi(x) = (µ ν g(x) + λ ν ) P k Gk,i(x) and Gk,i(x) is the corresponding gap function gradient vector, which is given by Gk,i(x) = ∇xi gk(x)+. (5.20) if there is no contact, the constraints can be disregarded (as the inequality is inactive), and the solution can be obtained for fi(x) = 0, (5.21) with λ ν = 0, which indicates none of the constraints is active. In the augmented Lagrange multiplier method, µ ν is initially set to a moderate value µ 0 . The problem is then solved using an appropriate method, such as Newton’s method, which forms the basis of an outer iteration. Upon convergence of this outer iteration, the Lagrange multiplier is updated by λ ν+1 = λ ν + µ ν g(x). (5.22) The penalty parameter µ ν is updated by a factor ζ according to µ ν+1 = ζ µν , (5.23) where ζ > 1 (typically ζ = 2). These outer iterations proceed until the solution is deemed accurate enough. The iterative solution of the system of equations Eq. (5.19) by the augmented Lagrange multiplier method using Newton’s method is given by ∆xi = [Jij (x m)]−1 ¯fj (x m), (5.24) x m+1 i = x m i − ∆xi , (5.25) where x m denotes the configurations at m-th Newton iteration and J is the Jacobian of the 72 system of equations. These equations are implemented using a matrix-free Newton method with the GMRES Krylov space linear solver. The approach requires only the product of the Jacobian of the dynamic equations (excluding the contact constraints) with a general vector y, which can be calculated using the first-order difference formula via ∇xfi(x)y = fi(x + ϵy) − fi(x) ϵ , (5.26) or the second-order central difference formula: ∇xfi(x)y = fi(x + ϵy) − fi(x − ϵy) 2ϵ , (5.27) where ϵ is updated within Newton’s iterations based on a very small number ϵo according to Eq. (5.28), which by default is set to 10−9 . ϵ = ϵo ||x||2 ||y||2 + 1 . (5.28) The Jacobian of Eq. (5.19) associated with the constraints, which can be calculated from the explicit formula, is given by J˜ ij = ∇xj ˜fi(x), (5.29) and ∇xj ˜fi(x) = µ ν Gk,j (x)Gk,i(x) + (µ ν gk(x) + λk) Hk,ij (x) (5.30) where the Hk,ij (x) is the Hessian matrix, which is given by Hk,ij (x) = ∇xi∇xj gk(x) (5.31) Since there is no need to assemble the Jacobian of the constraint explicitly, the part of 73 Eq. (5.29) concerned with the constraint only needs to be evaluated as J˜ k,ijyj = µ ν Gk,i(x)w + (µ ν gk(x) + λk) Qk,i, (5.32) where w = Gk,j (x)yj and Qk,i = Hk,ij (x)yj . Observation of the algorithm in practice shows that the above iteration usually terminates without problems for a given error tolerance, but it cannot assure that gk(x) ≤ 0. The numerical convergence process may result in some contacting vertices having a small positive value of the gap function. This is unsatisfactory since the calculation of the gap function depends on the ability to ensure that the mesh is not interpenetrated at the start of a time step. This prompted the creation of the acceleration step described below to ensure gk(x) ≤ −gϵ for some small gϵ > 0. 5.2.2 Iteration Accelerator with classical Lagrange multipliers Convergence of the augmented Lagrange multiplier method above needs to be followed by an accelerator to ensure the vertices are projected back to the admissible space of solutions (non-interpenetrating structure). This is done here by projecting the solution back to the admissible space using a classical Lagrange multiplier method. The main difficulty is to ensure that gk(x) ≤ −gϵ because only the positive branch of gk(x) is defined. After simplification to restrict to the active constraints, the constrained minimization of the classical Lagrange functional gives, fi(x˜) + λkGk,i(x˜) = 0, (5.33) gk(x˜) = 0, (5.34) The vertices’ coordinates are segregated between those involved in the calculation of gk(x), denoted by x˜, and the rest, which are not modified by this accelerator. Although the notation 74 gk(±˜x)+ is used to indicate that only active constraints are considered, the primary issue is that gk is only defined when positive. Therefore, the definition of the gap functions is extended by extrapolating from the latest known value and modifying the constraint to ensure that the vertices are always just above the contact surface by requiring that gk(x m) + Gk,i(xi − x m i ) + gϵ = 0, (5.35) where x m is the latest solution of the previous global iteration. No distinction is made here between x˜ and x since gk(x) = gk(x˜) by definition. The system of equations to be solved by Newton’s method is therefore formulated as:    ∇xj fi + λkHk,ij Gk,i Gk,i 0       ∆xi ∆λk    = −    fi(x m) + λkGk,i(x m) gk(x m) + Gk,i(x − x m i ) + gϵ    , (5.36) with x m+1 i = x m + ∆xi , (5.37) λ m+1 k = λ m k + ∆λk. (5.38) 5.3 Contact detection algorithm At the beginning of each time step, the contact algorithm assumes a geometric admissible mesh to start searching for contact. It begins by constructing coordinate-aligned bounding boxes (CABB) for all triangles based on their positions at xn and assumed displacements for xn+1 inferred from the velocity at xn. Since the final state at xn+1 is not known, to capture as many as possible contacting pairs, a rather conservative approach is utilized to compute the swept bounding boxes as it introduces a user-defined bounding box padding in x,y and 75 z direction. This padding length, denoted as lpad, is set as in the following: lpad = p|v|max(tn − tn+1), (5.39) where |v|max stands for the maximum velocity magnitude of all nodal velocities, and p is a user defined padding parameter, which is set to 1.0 by default. A global search list, LGP , is then generated to identify potential intersections between these CABBs, which is the same procedure as described in algorithm 1. Again similar to the approach taken in the explicit time scheme, the global contacting pairs are assumed to remain unchanged from state xn to xn+1. In the most cautious approach, the list LGP can be recalculated after each update in nodal displacements during every Newton iteration. However, this is generally avoided in most cases due to its potential to significantly increase computational costs. The next step involves determining the intersected triangles and evaluating the appropriate gap functions. The key distinction in the contact algorithm used here compared to the one in explicit time schemes is the use of a signed gap function to precisely identify the contacting triangle pairs. Clearly, it is not sufficient to determine actual interpenetration only from sign shift from the gap function of each triangle-triangle pair. In the example illustrated in Figure 5.2, when A intersects B, only cases IV and V show true intersections. However, a single signed function check cannot distinguish cases II and III from case IV and V if A is initially on the non-interpenetration side of B in all cases (In non-interpenetration side, all gap functions are positive). Therefore, a straightforward triangle-triangle intersection test can be used to identify triangle pairs that are in actual contact, thus eliminating cases I to III. 76 2.4. Algorithm The algorithm to compute the gap function has several steps. It is designed to not rely on previous time step information since this can be problematic. The steps in the algorithm below are here described in more detail. 1) First, we construct coordinate-aligned bounding boxes (CABB) for all triangles. We use the state at time step x n and assumed displacements (inferred from the velocity) to estimate where the triangles would go at the end of the time step, x n+1 . Obviously, we do not know what is the final state x n+1 when there is contact, and so this must be chosen conservatively. It is possible to recalculate CABB at every step of the Newton iteration, using that iterate x m, but this will necessitate the re-processing of the sliding nodes (see below) at every Newton iteration. It will be more costly, but it is a possible alternative. 2) Second, generate the list of CABB intersections, which I call broad search list (BSL). This is the broad phase of the algorithm, the same as in the explicit contact method. 3-6) Third, process all triangle vertices in the BSL to determine their relative position with respect to the corresponding triangle pair at x n . We use the sign of the tet volume as a measure of distance to a triangle (using the definition above). This does not entirely cover all possibilities because triangles can have tet volumes that appear to be of alternating signs while the triangles are not really intersecting. 7-16) Handle any sliding nodes that, while not really interpenetrated, may be left from the previous Newton iteration just slightly touching. This is not to mean that a triangle has interpenetrated but numerically appears like so. Figure 1 shows most possibilities of triangle B with respect to a triangle A (shown flat perpendicular to the figure). First we need to check there is interpenetration between A and B (check for triangle intersection). This will exclude all cases I, II, and III. The remaining cases IV and V (and VI not drawn where all vertices are below A) can be detected by the ϵ-displacement test. Let’s calculate the variation of the gap function of a vertex j with respect to a small increment in the position of the node in the direction of the normal, ˆn, of triangle A. The gap variation is equal to ∆g n j = ϵ G n j nˆ for some small ϵ. For any vertex j below triangle A, the gap function will be positive (in the orientation of the figure) and therefore g n j + ∆g n j will become negative. We must then check if at least one of the gap function of the other vertices of this triangle B is negative and did not change sign when performing the ϵ-displacement test and use that reference to decide that we should displace the vertex to the correct side of the space. Change the gap and position of the vertices accordingly. Also consider the reverse situation where the triangle B is below A. The value of ϵ is connected to the tolerance of the Newton method. Its value will need a few trials to see what works well. Warning: if the gap function of all three vertices of the triangle change sign when doing the ϵ-displacement test, then we should abort because the algorithm cannot tell where we are. Figure 1: Example of TB arrangement against TA at x n . 5 Figure 5.2: Possible intersection cases. Only cases IV and V are considered true intersections. Thus, in the local contact detection phase, a pairwise triangle intersection is first performed to identify the intersection triangle pairs in the local phase with the help of CGAL libraries. The resulting pairs will further be passed to a more detailed check to determine the valid 4-vertex tuple either in vertex-triangle or edge-edge where the proper gap function is calculated upon. All vertex-triangle and edge-edge pairs in the LGP are examined using the techniques described in section 5.1 to check for intersections at the x n+1 coordinates. As a result, all the intersecting vertex-triangle and edge-edge pairs, both formed in a 4-vertices tuple, are obtained and it is ready to compute individual gap function as well as its gradients and Hessian values described in the section 5.1. The Newton method is then introduced, beginning with the initialization of the cumulative gap and gradient to zero. An augmented system is therefore formulated and solved using iterative Newton’s method. In this process, the GMRES algorithm is employed to solve the linear system in each Newton iteration. When this Newton loop sees convergence, a solution of displacements and corresponding Λ vector is obtained. However, at this point, the geometric interpenetration has not been fully prevented yet due to the possibility that some nodes might be left just slightly touching after Newton iterations. This is not to mean that a triangle has interpenetrated but numerically 77 appears like so. To resolve this, a back-projection approach, denoted as ϵ-displacement, is then used to adjust gap functions accordingly. The cumulative gap and gradient is set to zero at the beginning of each local contact search. The variation of the gap function of a vertex j is calculated with respect to a small increment in the position of the node in the direction of the normal, ˆn, of triangle A. The gap variation is expressed as ∆g n j = ϵ Gn j nˆ for some small ϵ. For any intersecting vertex j from A , the gap function will be positive and therefore g n j + ∆g n j will become negative. It is necessary to check if at least one of the gap functions of the other vertices of this triangle B is negative and did not change sign when performing the ϵ-displacement test. This reference is used to decide that the vertex should be displaced to the correct side of the space. The gap and position of the vertices are then adjusted accordingly. The reverse situation where triangle B is below A is also considered. The value of ϵ is connected to the tolerance of the Newton method and may require several trials to determine an appropriate value. Notably, if the gap function of all three vertices of the triangle changes sign during the ϵ-displacement test, the algorithm is unable to determine the correct position and the process is aborted. This ϵ-projection algorithm is summarized in Algorithm 4. The overall algorithm for contact detection is summarized in Algorithm 5. 5.4 Preliminary results In this chapter, a simple triangle-triangle contact in a setting of head-on collision and sliding is considered to verify the validity of the implicit algorithm (see Figure 5.3 for corresponding scenarios). In all simulations, the acceleration iteration method as well as the ϵ-projection is adapted. 78 Algorithm 4 ϵ-projection 1: loop over all triangle pairs A − B 2: loop over all vertices j = {A1, A2, A3, B1, B2, B3} at x n ▷ this is the reference state 3: Calculate g n j = T n j and Gn j = ∇xT n j 4: end loop 5: if triangles intersect at state x n then 6: Calculate ∆g n j = ϵ Gn j nˆJ for some small ϵ ▷ J denotes the complementary triangle to vertex j 7: if g n j > 0 and g n j + ∆g n j < 0 and at least one g n k < 0 for vertices of same triangle then 8: Move vertex to other side, x n j ← x n j + ϵnˆJ and g n j ← g n j + ∆g n j 9: end if 10: if g n j < 0 and g n j − ∆g n j > 0 and at least one g n k > 0 for vertices of same triangle then 11: Move vertex to other side, x n j ← x n j − ϵnˆJ and g n j ← g n j − ∆g n j 12: end if 13: end if 14: end loop ▷ ensured non-interpenetrated geometry Algorithm 5 Contact resolution using signed gap function 1: for each Newton iteration do 2: Initialize g(x), the gradient vector G(x) = ∇xg(x), and hessian H(x) = ∇x∇xg(x) to zero 3: for each triangle-triangle pair in LGP do 4: if triangles intersect at state x n+1 then 5: for every vertex j do 6: Calculate g n+1 j = T n+1 j and G n+1 j = ∇xT n+1 j ▷ Eq. (5.1) and 5.3 7: if g n+1 j g n j < 0 then 8: Add this g n+1 j to global g(x), G n+1 j to G(x) and H n+1 j to H(x) 9: Calculate H n+1 j = ∇x∇xT n+1 j ▷ Eq. (5.4) 10: end if 11: end for 12: end if 13: end for 14: Formulate augmented linear system in Eq. (5.36) 15: Solve the augmented linear system in GMRES loop 16: if Converge then 17: Update x, λ, µ, and ϵ ▷ Eq. (5.28) and 5.37 18: end if 19: end for 79 5.4.1 Triangle-triangle impact The triangle mesh configurations are exactly the same for both scenarios in Figure 5.3: each triangle has 4 triangular elements, thickness 0.1 mm, density 500 kg/m3 , Young’s modulus 5 × 107 kg/ms2 , Poisson ratio 0.33. In a head-surface scenario (left in the Figure 5.3, a initial velocity of 10 m/s, perpendicular to the surface of blue triangle on bottom, is set to the red triangle on the top while the bottom blue triangle remains stationary. Similarly, in sliding triangle scenario (right in the Figure 5.3), the red triangle on the top is set with an initial velocity of the same magnitude, parallel to the surface of blue triangle on the bottom, as well as an downward external acceleration of -1000 m/s2 to the vertices at the bottom rim. 10 m/s 10 m/s Figure 5.3: Initial setup for head-on impacting (left) and sliding (right) respectively. The simulation result of the head-surface contact is illustrated in Figure 5.4. The results display the contact response of the red triangle, where the blue triangle is dented due to the impact. Subsequently, the blue triangle bounces back, causing the red triangle to move upwards. After several instances of bouncing on and off, the two triangles eventually separate. 80 0.0037s 0.0074s 0.011s 0.015s 0.018s 0.022s 0.026s 0.030s Figure 5.4: Simulation results of impacting triangle. In the triangle sliding case, the contact response exhibits similar behavior as the blue triangle becomes slightly folded when the red triangle slides to its center position. Due to the external accelerations applied to the red triangle, it continues to move downward after sliding out of the blue triangle (see Figure 5.5). 0.0028s 0.0056s 0.0085s 0.011s 0.014s 0.017s 0.020s 0.023s Figure 5.5: Simulation results of sliding triangle. 81 Chapter 6 Programming techniques of contact algorithms This chapter delves into the implementation of contact algorithms utilizing advanced data structures, concurrent inter-processor data transmission, and OOP in low-level software design incorporating templated libraries and user-centric functionalities in C++. The sophisticated data structures are pivotal in managing the complex geometric entities involved in contact mechanics, ensuring efficient calculations. The concurrent inter-processor data transmission mechanism enhances the performance and scalability of these algorithms, facilitating real-time computations across multiple processors. Additionally, the templated library approach and user-centric functions highlight their roles in creating flexible, reusable, and easily maintainable codes that can be tailored to meet specific user needs. These advanced methodologies aim to address the inherent challenges in contact algorithm implementation, providing robust solutions that leverage modern computational techniques. The chapter is organized to introduce programming techniques from low-level concepts to high-level details. At the low level, the design leverages modern C++ STL, appropriate data structures and algorithms (DSA), and smart pointers to enhance sorting efficiency and memory usage. In the high-level design, the focus is on extending the basic finite element shell solver, as detailed in [47, 53], using polymorphism (see section 6.2). 82 6.1 Data structures and algorithms of geometric contact entities This section discusses two key strategies that enhance the performance and maintainability of such algorithms: the integration of shell and contact solvers on a unified dataset and the use of hashmaps for rapid data mapping in contact search. The execution of both the shell solver and the contact solver on the same dataset. This approach eliminates the need to maintain separate datasets for each solver, reducing redundancy and minimizing the potential for errors. By sharing a unified dataset, the solvers can operate more cohesively, ensuring consistency in data handling and reducing the overhead associated with synchronizing multiple datasets. This integration simplifies the overall software architecture, making it easier to manage and update the data structures as the algorithms evolve. Another critical aspect of efficient data management in contact algorithms is the use of hashmaps for fast data mapping. Hashmaps offer a highly efficient way to store and retrieve data, allowing for constant time complexity in average cases. This speed is particularly beneficial when dealing with large sets of geometric entities, where quick access to data is essential for maintaining the performance of the solvers. By employing hashmaps, the implementation can ensure that the mapping of geometric entities and their associated attributes is both fast and reliable, significantly enhancing the responsiveness and scalability of the contact algorithms. The basic shell solver ShellSolver, developed in [53], includes methods for Finite Element Method (FEM) operations and geometric entities. Specifically, geometric entities such as vertices, edges, and elements are stored in C++ STL containers: vector<vertex*>, vector<edge*>, and vector<element*>, respectively. These containers are designated as protected fields to provide direct access to derived classes. When dealing with contact, a specialized contact solver, ContactSolver, is derived from ShellSolver, encapsulating data and algorithms specific to contact interactions. A 83 key feature of this solver is the contact elements container, elementsContact , which is a vector<element*> holding all elements relevant for contact search. Each element in this container corresponds to an element allocated in the parent class. In the current contact algorithm, elementsContact is initialized at the start and remains unchanged throughout the simulation, as dynamic load balancing has not yet been implemented. Global contact detection is performed on elementsContact, and the results are stored in vector<pair<element*, element*>>. To avoid redundant contact checks, this list is further decomposed into a list of element-vertices elemVtxList (unordered map<element*, SvertexSet>) for vertex-triangle test, and a list of edge-edge pairs edgeEdgeList (set<pair<edge*, edge*>, edgePairComp>) for edge-edge tests. Here, SvertexSet is a set<vertex*> used for vertex-triangle checks, and edgePairComp is a custom comparator to prevent duplicate edge-edge checks (see 6.1). Listing 6.1: customized comparator for edge pairs. bool edgePairComp ( co n s t pai r<edge ∗ , edge∗> &ep1 , co n s t pai r<edge ∗ , edge∗> &ep2 ) co n s t { /∗ compare by a d d r e s s ∗/ co n s t edge ∗ e1 = min ( ep1 . f i r s t , ep1 . second ) ; co n s t edge ∗ e2 = min ( ep2 . f i r s t , ep2 . second ) ; /∗ r e t u r n ∗/ r e t u r n ( e1 < e2 ) | | ( ( e1 == e2 ) && max( e1 . f i r s t , e1 . second ) < max( e2 . f i r s t , e2 . second ) ) ; } In local contact detection, a specialized data structure, contactPair, encapsulates the contact pairs. From a mathematical perspective, edge-edge and vertex-pair contacts are treated equivalently (refer to the data structure declaration for contactPair in 6.2). Valid contact pairs identified during the local contact check are stored in a container set<contactPair, contactPairComp>. A custom comparator, contactPairComp, is employed to maintain the sorted order of contact pairs, which is crucial for efficient point-to-point communication and data exchange during local phase parallelization (see 6.3). Finally, the class ContactSolver is summarized in Figure 6.1. 84 Listing 6.2: contact pair data structure declaration. s t r u c t c o n t a c tP ai r ( v e r t e x ∗v0 , v e r t e x ∗v1 , v e r t e x ∗v2 , v e r t e x ∗v3 ) { co n s t v e c to r<v e r t e x∗> v e r t e x = {v0 , v1 , v2 , v3 } ; /∗ v e r t i c e s c o n t ai n e r ∗/ double gap ; /∗ gap f u n c ti o n val u e ∗/ double g r a d i e n t [ 3 ] ; /∗ g r a di e n t v e c t o r ∗/ double h e s s i a n [ 9 ] ; /∗ h e s si a n ma t rix ∗/ double t ∗ ; /∗ co n ta c t time ∗/ } Listing 6.3: customized comparator for contact pairs. bool contactPairComp ( co n s t c o n t a c tP ai r &cp1 , co n s t c o n t a c tP ai r &cp2 ) co n s t { auto i t 1 = cp1 . b egin ( ) ; auto i t 2 = cp2 . b egin ( ) ; /∗ compare two co n ta c t p ai r s , r e t u r n t r u e i f the s m a l l e s t v e r t e x i n cp1 i s l e s s than t ha t i n cp2 ∗/ i f (∗ i t 1 == ∗ i t 2 ) { i f ( ∗ ( i t 1 + 1 ) == ∗( i t 2 + 1 ) ) i f ( ∗ ( i t 1 + 2 ) == ∗( i t 2 + 2 ) ) r e t u r n ∗( i t 1 + 3 ) < ∗( i t 2 + 3 ) ; e l s e r e t u r n ∗( i t 1 + 2 ) < ∗( i t 2 + 2 ) ; e l s e r e t u r n ∗( i t 1 + 1 ) < ∗( i t 2 + 1 ) ; } e l s e r e t u r n ∗ i t 1 < ∗ i t 2 ; } 85 contact_shell_inh - Untitled ShellSolver # vertices_: vector<vertex*> # edges_: vector<edge*> # elements_: vector<element*> - readMesh(File meshFile) - computeNodalMass() - internalForce() - externalForce() ...... ContactSolver - elementsContact_: vector<element*> - globalList_: vector<pair<element*, element*>> - elemVtxList_: unordered_map <element*, SvertexSet > - edgeEdgeList_: set<pair<edge*, edge*>, edgePairComp>> - narrowList_: set<contactPair, contactPairComp> - computeGap() - computeGapGradientHessian() + setContactElements() + globalContactDetection() + localContactDetection() + resolveCollision() ...... 1 0 Figure 6.1: UML diagram of contact solver. 6.2 Template library and user-centric approach To enhance the existing shell software, an wrapper class, ShellManager, is then developed which wraps key functions from ShellSolver. From user end-point view, this class provides interfaces to invoke key functions such as initialize, advanceTime, output, and checkPoint, for controlling the simulation flow. Further, to easily maintain parallel shell solver, a parallel shell handler class ShellParallelManager is inherited from ShellManager, where data exchange and synchronization are perform within each time advance step. And similarly to the design paradigm in shell solver, an additional manager class ContactManager, inherited from ShellManager, is employed in contact-included simulations. Here, different polymorphism techniques, overloading and overriding, are widely used through virtual functions and inheritance. The contact parallelization is extended from the parallel shell solver ShellParallelManager and ContactManager in class ContactParallelManager, because it needs essential functionalities in both classes. These inheritance relations are described in the Figure 6.2. 86 New Diagram - Untitled ShellManager + timeStamp_ + outputFile_ + initialize(): void + advance(): void + visualize(): ofstream ContactManager ShellParallelManager ContactParallelManager 1 0 1 1 1 0 0 0 ContactSolver ShellSolver Figure 6.2: UML diagram of driver class inheritance. In user stand point, one can directly create his or her own customized manager class derived from any of the four manager classes described above. On a high level, the class to derive from is determined by C++ macro combinations, PARALLEL and CONTACT, where the possible outcomes are listed as below: Listing 6.4: user defined driver class. #i f d e f CONTACT #i f d e f PARALLEL t y p e d ef Con tac tPa rallelManage r<Con tac tSol ve r> Manager ; #e l s e t y p e d ef ContactManager<Con tac tSol ve r> Manager ; #e l s e #i f d e f PARALLEL t y p e d ef S h ellPa rall elMa nag e r Manager ; #e l s e 87 t y p e d ef ShellManager Manager ; #end c l a s s MyManager : p u bli c Manager { // cu s tomized f u n c t i o n s } ; 88 Chapter 7 Results 7.1 Simulation results This chapter presents results for geometric contact resolution both in explicit time scheme as well as in implicit time scheme through three distinct examples. The first two examples utilize the explicit time scheme, while the third employs the implicit time scheme. The first example evaluates the performance of the parallelization algorithm. The parallel contact algorithm constructs a specialized communication map (see Section 4.1.2) to ensure load balancing among processes during the broad phase. The performance of this parallel approach is compared to the system’s performance without contact detection, as detailed in section 7.1.1. The second example, presented in section 7.1.2, demonstrates a simple parachute folding scenario, validating the geometric contact resolution implementation within the implicit time integration scheme. Finally, the last example illustrates a complex parachute deployment mechanism driven by cable motion, featuring significant self-contact with many folded structures within the material. 89 10.7cm 4.9cm 10.0cm Figure 7.1: fabric-folding simulation at t = 0 with front view on the left, and side view on the right. The bend has an inner radius of 3.45 mm. Figure 7.2: fabric mesh partition for a 32-processor MPI task. 7.1.1 Fabric folding The paper in [63] explores the simulation of a fabric suspended from one corner, employing two distinct domain decompositions to handle material simulations and contact detection separately. Notably, the simulation in their study utilizes a relatively small mesh size, comprising only around 1000 mesh elements, and lacks a scaling study. To expand upon this research, particularly in exploring the scalability of our algorithm in managing a significantly larger number of simultaneous contacts, we introduce a fabric folding simulation in this section: a model of multi-layer fabric folding (see geometry in Figure 7.1) is considered, where an equal and opposite velocity, 0.1 m/s, is applied on the outer-most layers of two sides to fold the structure. The material properties are: density 103 kg/m3 ; Young’s modulus 105 kg/m s2 ; Poisson ratio 0.33; thickness 10−4m. The corresponding finite-element mesh has 499,392 triangular elements. The simulation results are shown here in Figure 7.3 with 32 MPI processors (see corresponding mesh partition in Figure 7.2). As time evolves, the fabric continues to be squeezed from both sides symmetrically, as the initial velocities are of the same magnitude, until the inner-most layers stick together (0.016-1.07 s in Figure 7.3). Then one after another 90 the layers separate (2.14-4.28 s in Figure 7.3). Eventually, after separating the two innermost layers, initial curvatures are largely flattened and the fabric is left with wrinkles from deformations (2.68 s in Figure 7.3). To analyze the scales of contact occurring in the simulation, we monitor the number of contacting pairs in both global and local phases throughout the entire duration (refer to the full simulation in Figure 7.5). In the case of global contacting pairs, also known as intersecting CABB pairs, their count remains within the range of 105 to 106 under the default EXCL setting. Conversely, the local contacting pairs show a fluctuation between 0 and 100. The onset of contact sees a sudden surge in both global and local pairs, occurring almost immediately at the start (at 0.016 s, as seen in Figure 7.3). This intensive contact involves both sides of the material simultaneously, as indicated by the green box in Figure 7.4. Subsequently, the material continues to compress, leading to a substantial increase in the number of global contacting pairs and frequent local contact events until around 0.8 s. Following this, as the material begins to expand and separate layers, the number of contacts significantly decreases. At this point, sporadic local contact events occur from time to time. Finally, the material experiences zero instances of local contact, indicating that it has been fully flattened. To examine the parallelization performance of our contact algorithms, we concentrate on three specific cases, denoted as case a, b, and c, each characterized by significantly different contact scales at different points in time, which follow in the order of decreasing local contacts. Case a to c were sampled at approximately 0.016 s, 0.141 s, and 0.001 s, each spanning time windows of 10, 500, and 1000 steps, with average local contacting pairs of 178, 5, and 0, respectively, as depicted in Figure 7.5 for cases a-c. 91 s 0.016s 0.14s 0.53s 0.80s 1.07s 1.34s 1.61s 1.87s 2.14s 2.41s 2.68s Figure 7.3: Results of fabric-folding simulation. 92 Figure 7.4: Immersive contact happens when two sides of material touches its inner layer simultaneously at 0.016 s. The parallelization of contact problems can easily lead to load imbalance and thus reduce parallel efficiency. A key metric for evaluating the effectiveness of a parallel program is strong scaling, which measures how increasing computing tasks helps reduce the execution time of a given problem, and thus has been widely used in various applications [64, 65]. For this reason, we chose this fabric folding case to conduct a scaling test in strong form since it has a large mesh size and has the most contact-intensive occurrences where the layers stick together from both sides. The scaling test, employing up to 256 processors and including a pure FE shell solver with contact disabled and cases a to c as previously mentioned, is summarized in Figure 7.6 with the corresponding normalized CPU wall time results. This normalized CPU wall time is computed as the standard grind time g: g = Texe N ∗ cycles, (7.1) where Texe is the execution wall time for the sampled window, N is the total number of mesh elements, cycles is the number of sampled time steps. In all scenarios, the wall time decreases as more MPI tasks are employed, and notably, the greater the contact intensity, the more time it takes to complete the calculation. Further93 more, the time scales for these cases exhibit substantial variation: the FE solver completes a step without contact considerations in just a few nanoseconds, while in contrast, case a, the most contact-intensive, requires more than 20 minutes in a serial application and around 30 seconds in a 256-processor application. It is evident that the local contact handling is the most time-consuming aspect of the algorithm, as increasing the count of local pairs by 10x leads to a 106x increase in execution time to complete one step when transitioning from case b to c, given that their number of global contacting pairs remain similar. The decrease in execution time associated with increased computing tasks can be attributed to the distribution of workload among processors, as illustrated in Figure 7.7 and 7.8 for global and local cases, respectively. As anticipated, case c exhibits the most evenly balanced load distribution. In contrast, cases a and b display less balanced load distributions, with a significant concentration of tasks on the latter half of processors, where contact points are concentrated. step 25000 50000 75000 100000 0 100 101 102 103 104 105 106 full simulation broad phase narrow phase 1121 1124 1127 1130 0 100 101 102 103 104 105 106 case a 10000 10200 10400 0 100 101 102 103 104 105 106 case b 0 250 500 750 1000 0 100 101 102 103 104 105 106 case c time (s) 0.695 1.389 2.084 0.016 0.141 0.001 Figure 7.5: Number of contact pairs (logscale) for global and local phase looking at different time windows. 94 1 4 8 16 32 64 128 256 NCPU 10−1 101 103 105 107 109 grind time ( µs) FE shell case a case b case c Figure 7.6: CPU wall time (logscale). 95 0 rank 500000 600000 serial case a case b case c 0 1 2 3 150000 200000 4p 0 2 4 6 75000 100000 8p 0 5 10 15 40000 60000 16p 0 10 20 30 20000 30000 32p 0 20 40 60 10000 20000 64p 0 20 40 60 80 100 120 5000 10000 15000 128p 0 50 100 150 200 250 2000 4000 6000 8000 256p 0 rank 500000 600000 serial case a case b case c 0 1 2 3 150000 200000 4p 0 2 4 6 75000 100000 8p 0 5 10 15 40000 60000 16p 0 10 20 30 20000 30000 32p 0 20 40 60 10000 20000 64p 0 20 40 60 80 100 120 5000 10000 15000 128p 0 50 100 150 200 250 2000 4000 6000 8000 256p 0 rank 500000 600000 serial case a case b case c 0 1 2 3 150000 200000 4p 0 2 4 6 75000 100000 8p 0 5 10 15 40000 60000 16p 0 10 20 30 20000 30000 32p 0 20 40 60 10000 20000 64p 0 20 40 60 80 100 120 5000 10000 15000 128p 0 50 100 150 200 250 2000 4000 6000 8000 256p Figure 7.7: Cross-processor distribution of global contacting pairs for case a (averaged over 10 steps). (see in color print) 96 0 rank 0 100 serial 0 1 2 3 0 100 4p 0 2 4 6 0 100 8p 0 5 10 15 0 50 16p 0 10 20 30 0 20 32p 0 20 40 60 0 20 64p 0 20 40 60 80 100 120 0 10 20 30 128p 0 50 100 150 200 250 0 5 10 15 256p 0 rank 0 100 serial 0 1 2 3 0 100 4p 0 2 4 6 0 100 8p 0 5 10 15 0 50 16p 0 10 20 30 0 20 32p 0 20 40 60 0 20 64p 0 20 40 60 80 100 120 0 10 20 30 128p 0 50 100 150 200 250 0 5 10 15 256p 0 rank 0 100 serial 0 1 2 3 0 100 4p 0 2 4 6 0 100 8p 0 5 10 15 0 50 16p 0 10 20 30 0 20 32p 0 20 40 60 0 20 64p 0 20 40 60 80 100 120 0 10 20 30 128p 0 50 100 150 200 250 0 5 10 15 256p Figure 7.8: Cross-processor distribution of local contacting pairs for case a (averaged over 10 steps). 97 7.1.2 Simple parachute folding This section delves into a more intricate scenario than the one discussed in section 5.4, which involves a simple drogue parachute folding driven by external ramping cable forces. The simulation setup is illustrated in Figure 7.9. The parachute has a radius of 3 m and a depth of 1 m, supported by 20 external cables. The suspension lines are 5 m in length, extending from the parachute. Similar to the cable modeling in section 7.1.3, the cable can carry string force. The material properties of the parachute are: density 500 kg/m3 , Young’s modulus 5 × 108 kg/m s2 , Poisson ratio 0.33, and thickness 10−4m. The material properties for the cable embedded in the shell are: density 0.002 kg/m, cable stiffness 3 × 104 kgm/s 2 , and the suspension line has a density of 0.002 kg/m and cable stiffness 6 × 104 kgm/s 2 . The parachute mesh consists of 8860 vertices and 17360 elements. The suspension lines are modeled as a single piece. 5 m 1 m 6 m 0.6 m Figure 7.9: The configuration of a simple parachute with cables. Here, a ramping force mechanism is employed to fold the parachute within the suspension lines, consisting of three phases. In the first phase, the cable force magnitude gradually increases to 5000 N within 0.1 s. In the second phase, the force maintains a plateau at 5000 N until 0.115 s. In the final phase, the force quickly diminishes to 0 within 0.001 s. The 98 temporal force development can be expressed mathematically as: f(t) =    5000 t 0.1 2 for 0 ≤ t ≤ 0.1, 5000 for 0.1 < t ≤ 0.115, 5000 1 − t−0.115 0.001 2 for 0.115 < t, (7.2) which is also summarized in Figure 7.10. Figure 7.10: Temporal development of cable force. As time progresses, the parachute undergoes folding due to the external cable force (see Figure 7.11). This force causes the parachute to fold rapidly, transitioning quickly into a highly collapsed state. As the parachute continues to collapse, the contact events become more frequent, with approximately 30 contact events captured within each time step after 0.1 s. The algorithm effectively captures the non-interpenetration requirement 0.02s 0.04s 0.08s 0.11s 0.115s 0.118s Figure 7.11: Simulation results of simple parachute folding driven by cable force. 7.1.3 Parachute deployment The high-fidelity simulation of parachutes for space exploration is maturing at an accelerated pace, both in the subsonic and supersonic regimes [66, 67, 68]. However, there are still 100 significant problems related to early time dynamics, such as unbagging, deployment, and early inflation. Particularly, it is recognized that the early dynamics might result in excess loads on the fabric, seams, and suspension lines, which one wants to measure to decrease uncertainty and increase safety margins. Here, we concentrate on the early inflation of partially folded parachutes, where there are major geometrical constraints in producing suitable configurations from which coupled fluid-structure interaction (FSI) simulation can be carried out. The geometry of the parachutes is typically characterized in terms of their stressfree construction. However, this nominal shape is inappropriate for deployment simulation, and to facilitate deployment research, the parachute must be partially folded or collapsed. It is necessary to use a deformation and/or loading technique to bring the built geometry to a desired partially inflated condition to generate partially deflated parachutes. Since the likelihood of contact increases significantly during folding, the key computational difficulty, in this case, is the efficient and reliable implementation of a self-contact algorithm. Existing computational contact approaches are quite efficient in working with three-dimensional objects with a finite volume[69, 9, 1, 8, 70, 71, 72]. Thin materials, like fabrics, present a slightly more complex scenario since their geometries are surfaces in three-dimensional space and lack a pre-defined interior and exterior. Since it is typically impossible to reverse an incorrect deformation once the surface interpenetrates, contact must be managed effectively at every instance, which poses significant reliability issues. Unlike the other two simulations, the parachute deformation is driven by a specified deployment protocol to simulate the parachute behavior with payload. The protocol includes a few parts unique to parachute application, such as cable lines and reefing lines, to control the parachute’s open area. The simulations use shell mechanics for the fabric, cable mechanics for the suspension and reefing lines, and self-contact algorithms for contact detection and handling. The canopy has been fitted with a reefing line that can be expanded, contracted, and cut as needed, see Figure 7.12. We assume there is a single riser attached to a fixed point (where the payload would be) with 24 suspension lines extending from the riser to the 101 rim of the canopy. The riser line has the properties of the aggregated 24 suspension lines. The finite-element model of the parachute has 7056 triangular shell elements, 1511 finite cable segments in the shell, and 506 external finite cable segments in the suspension lines, riser and reefing lines. The material properties of the fabric are the same as hemispheres (see section 4.2.1). The suspension lines, reefing line, and riser line have linear density 0.0175 kg/m, 0.0248 kg/m, and 0.42; kg/m, respectively, and stiffness is 6.31 × 105 kg m/s 2 , 7.81 × 105 kg m/s 2 and 1.51 × 107 kg m/s 2 , respectively. We also use a special material type for the cables along the canopy skirts, which has linear density 0.0301 kg/m and stiffness 5.52 × 104 kg m/s 2 . Figure 7.12: Main parameters of drogue parachute constructed geometry. 7.1.3.1 Deployment protocol The strategy devised to partially collapse the original constructed geometry follows a loading sequence described schematically in Figure 7.13, and consists of the following processes: • A prescribed load pressure, pi , is applied on the inner surface of the chosen panels of the canopy. In the present simulations we apply a constant load either on: a) the end disk, b) the end disk and next band of panels, or (c) the end disk and subsequent two bands of panels of the parachute. These are denoted by 1-, 2- or 3-band loading, respectively, throughout this section. • The initial length of the reefing line is equal to the perimeter of the canopy rim, of diameter D = 4.48 m, in its constructed geometry specification. This length is 102 substantially larger than the actual length of the reefing line in the first and second stages of reefing which have diameters of D1 = 0.656D and D2 = 0.9D, respectively. To achieve a geometrical configuration with the actual diameter of the reefing line we proceed by reducing its diameter, in four steps from t = 0 and up to t ∗ , by 90% each step (for a total of 0.9 4 ≈ 0.656). In between steps, we allow the parachute to deform dynamically to accommodate this change in the reefing line diameter. Each shortening step of the reefing line pulls the skirt inwards (radially). Therefore, we freeze the canopy by setting all velocities to zero shortly after shrinking the reefing line to prevent the parachute from collapsing owing to the inertia induced by the abrupt shortening of the reefing line. We have experimented with two protocols during this stage. In one case we impose the uniform loading pressure pi during the whole duration of the simulation while in a second strategy we imposed an over-pressure po during the first step of shortening from 0 < t ≤ to; this is referred to as the over-pressured startup. • After the reefing line has achieved its nominal first-stage reefing diameter D1 at t = t∗, the canopy is allowed to deform dynamically up to t = t∗ + t1. The stresses and loads on the suspension lines are allowed to follow the normal dynamics without any changes to the geometry or topology of the parachute. • Immediately after t = t∗+t1, the reefing line is extended instantaneously to the secondstage reefing diameter D2, and the simulation is carried out up to t = t∗+t1+t2 without further modification. • At t = t∗ + t1 + t2, the reefing line is removed and the parachute is allowed to expand to its fully deployed configuration. The values of the applied load and timings can be varied. The inflation pressure pi used here is approximately the transmural pressure that a parachute deployed at 4000 m altitude. We conducted our experiments in three overloading settings for po = 10, 20, and 40 kPa, respectively. The other parameters we used are pi = 4 kPa, t∗ = 0.0141 s, t1 = 0.00703 s, 103 and t2 = 0.00703 s. Two overpressure loadings were used, one applied over t0 = 0.00351 s and one over the full initial shrinking of the reefing line and including t1, i.e., t0 = 0.0211 s. These are denoted as burst and full overloading, respectively. Figure 7.13: General inflation pressure timeline evolution. 7.1.3.2 Simulation results The simulations show typical behavior where the pressurized canopy induces axial elongation and subsequent partial collapse of the parachute, see figures 7.14, 7.15 and 7.16. Shortly after simulation startup, there is broad self-contact of the structure and some panel folding, which needs to be resolved correctly for the parachute to be in a physical state, free of interpenetration. The suspension lines display their own dynamics, which are not fully displayed in the figures, but evidently couple with the geometry displayed by the canopy. The figures show that the chosen load protocol has a large influence on the state of the parachute, which is to be expected. Evidently, different load protocols can be used to generate even more varied geometries, and the possibilities are very large. Here, we focus on the loading protocol sketched in Figure 7.13. To aid in the description of the results, we will refer to each simulation by the row-column ordering in the figures, e.g., R2-C1 denotes the second row and the first column. The first row and column correspond to the top row and 104 leftmost column, respectively. Figure 7.14: Results for 1-band loading scenario. From top to bottom: at the end of stage 1, at the end of stage 2, after disreefing, where t = t∗ + t1, t∗ + t1 + t2, t∗ + t1 + t2 + 0.00703 s, respectively. From left to right: first no overloading; then full overloading and burst overloading at po = 10, 20, and 40 kPa, respectively. Figure 7.15: Results for 2-band loading scenario; arrangement as in Figure 7.14. 105 Figure 7.16: Results for 3-band loading scenario; arrangement as in Figure 7.14. Loading only the back disk (1-band) of the parachute generates a more elongated geometry, as can be seen by comparing the first row of all three figures 7.14, 7.15 and 7.16. Loading 3 bands tend to produce more hemispherical geometries, while loading 2 bands is more cylindrical in nature, although there is variability depending on the overloading. The largest overloading pressures result, invariably, in the most collapsed geometries, see R3-C6 results. Some patterns are interesting, like 1-band R3-C4 and C6, and 3-band R2-C6 and R3-C6, showing a conically expanded skirt and overinflated disk (first band). Figures 7.17, 7.18 and 7.19 show histograms of von Mises stresses at the corresponding times and for all cases discussed before. The stress is essentially bounded between 0.1 MPa and 104 MPa for all cases. The distribution is more or less Gaussian with a mode about 100 MPa and a few exceptions that show bimodal distributions, such as in 1-band loading R2-C5 and C7 and 3-band loading R1-C1. The distribution of stresses is informative but obviously not as specific as that coming from an fluid-solid interface simulation. Furthermore, it is well known that the details of the constitutive law of the fabric have a large influence on the calculated stresses, but not as much on the geometry of the parachutes since that is 106 regulated by the boundary conditions with the flow. Figure 7.17: Histograms of von Mises stress (σv) for 1-band loading in linear-log coordinates: lower and higher values of σv are 0.1 MPa and 104 MPa, respectively; arrangement as in Figure 7.14. 107 Figure 7.18: Histograms of von Mises stress (σv) for 2-band loading in linear-log coordinates; arrangement as in Figure 7.14. Figure 7.19: Histograms of von Mises stress (σv) for 3-band loading in linear-log coordinates; arrangement as in Figure 7.14. 108 Chapter 8 Conclusion To achieve high-fidelity simulations of thin-shell deformations in computational mechanics, practical solutions are necessary to resolve material self-contact effectively. In this context, a geometrically accurate contact algorithm is proposed, ensuring physically accurate simulations of material folding and collapsing by significantly reducing contact search efforts through a coupled global-local phase searching strategy. This contact algorithm has been integrated into both explicit and implicit time schemes. In the explicit time scheme, the algorithm, based on the DCR approach outlined in [9], nearly perfectly conserves energy and demonstrates robustness in capturing various contact structures, as well as scalability in parallelization. The simulation results highlight the algorithm’s success in providing realistic shell deformation configurations by preventing interpenetration and ensuring good load balancing in parallel computations. However, the contact resolution method in the explicit time scheme faces challenges, particularly in intensive contact scenarios, where the recursive contact time advancement scheme becomes inefficient. To address these issues, the augmented Lagrange multiplier method is proposed as a promising solution, allowing the augmented system to be solved in a single step without needing to consider internal contact events. In the implicit time scheme, the contact algorithm retains the global-local phase contact search strategy, with significant modifications in the local phase to accurately compute the active contact constraints required by the Lagrangian formulation. The algorithm’s effectiveness is demonstrated through a simple 109 parachute folding example, where it successfully captures multiple self-contact structures. The research is ongoing, with several areas for improvement: 1) optimizing the parallel contact resolution in the explicit time scheme to improve load balancing and efficiency, 2) future work could focus on parallelizing contact resolution in the implicit time scheme, and 3) more rigorous testing is required to verify the correctness of contact resolution in the implicit time scheme, particularly concerning mesh convergence and energy conservation properties. 110 Bibliography [1] MW Heinstein, SW Attaway JW Swegle, and FJ Mello. “A general contact detection algorithm for finite element analysis”. In: WIT Transactions on Engineering Sciences 1 (1970). [2] MW Heinstein et al. A general-purpose contact detection algorithm for nonlinear structural analysis codes. Tech. rep. Sandia National Labs., Albuquerque, NM (United States), 1993. [3] LM Taylor and DP Flanagan. PRONTO 3D: A three-dimensional transient solid dynamics program. Tech. rep. Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 1989. [4] David J Benson and John O Hallquist. “A single surface contact algorithm for the post-buckling analysis of shell structures”. In: Computer methods in applied mechanics and engineering 78.2 (1990), pp. 141–163. [5] Ted Belytschko and Mark O Neal. “Contact-impact by the pinball algorithm with penalty and Lagrangian methods”. In: International Journal for Numerical Methods in Engineering 31.3 (1991), pp. 547–572. [6] Yan Gu et al. “Efficient BVH construction via approximate agglomerative clustering”. In: Proceedings of the 5th High-Performance Graphics Conference. 2013, pp. 81–88. [7] Matthias Teschner et al. “Optimized spatial hashing for collision detection of deformable objects.” In: Vmv. Vol. 3. 2003, pp. 47–54. [8] Jonathan Boustani et al. “A Numerical Investigation of Parachute Deployment in Supersonic Flow”. In: AIAA Scitech 2020 Forum. 2020, p. 1050. [9] Fehmi Cirak and Matthew West. “Decomposition contact response (DCR) for explicit finite element dynamics”. In: International Journal for Numerical Methods in Engineering 64.8 (2005), pp. 1078–1110. [10] Kenneth Langstreth Johnson. Contact mechanics. Cambridge university press, 1987. [11] Thomas JR Hughes et al. “A finite element method for a class of contact-impact problems”. In: Computer methods in applied mechanics and engineering 8.3 (1976), pp. 249–276. 111 [12] Nicholas J Carpenter, Robert L Taylor, and Michael G Katona. “Lagrange constraints for transient finite element surface contact”. In: International journal for numerical methods in engineering 32.1 (1991), pp. 103–128. [13] D Boukari and AV Fiacco. “Survey of penalty, exact-penalty and multiplier methods from 1968 to 1993”. In: Optimization 32.4 (1995), pp. 301–334. [14] Juan C Simo. “A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part II: computational aspects”. In: Computer methods in applied mechanics and engineering 68.1 (1988), pp. 1–31. [15] Gianni Gilardi and Inna Sharf. “Literature survey of contact dynamics modelling”. In: Mechanism and machine theory 37.10 (2002), pp. 1213–1239. [16] Shaofan Li et al. “A meshfree contact-detection algorithm”. In: Computer methods in applied mechanics and engineering 190.24-25 (2001), pp. 3271–3292. [17] G Battiato and CM Firrone. “Reduced order modeling of large contact interfaces to calculate the non-linear response of frictionally damped structures”. In: Procedia Structural Integrity 24 (2019), pp. 837–851. [18] Lei Peng, Zhi-Qiang Feng, and Pierre Joli. “A semi-explicit algorithm for solving multibody contact dynamics with large deformation”. In: International Journal of NonLinear Mechanics 103 (2018), pp. 82–92. [19] Chuanqi Liu and WaiChing Sun. “ILS-MPM: An implicit level-set-based material point method for frictional particulate contact mechanics of deformable particles”. In: Computer Methods in Applied Mechanics and Engineering 369 (2020), p. 113168. [20] Bruce Hendrickson and Karen Devine. “Dynamic load balancing in computational mechanics”. In: Computer methods in applied mechanics and engineering 184.2-4 (2000), pp. 485–500. [21] David Harmon et al. “Asynchronous contact mechanics”. In: ACM SIGGRAPH 2009 papers. 2009, pp. 1–12. [22] Samantha Ainsley et al. “Speculative parallel asynchronous contact mechanics”. In: ACM Transactions on Graphics (TOG) 31.6 (2012), pp. 1–8. [23] He Liu et al. “A novel GPGPU-parallelized contact detection algorithm for combined finite-discrete element method”. In: International Journal of Rock Mechanics and Mining Sciences 144 (2021), p. 104782. 112 [24] M Mohammadnejad et al. “GPGPU-parallelized 3D combined finite–discrete element modelling of rock fracture with adaptive contact activation approach”. In: Computational Particle Mechanics 7 (2020), pp. 849–867. [25] Kenneth L Johnson. “One hundred years of Hertz contact”. In: Proceedings of the Institution of Mechanical Engineers 196.1 (1982), pp. 363–378. [26] Olek C Zienkiewicz and Robert L Taylor. The finite element method set. Elsevier, 2005. [27] JO Hallquist, GL Goudreau, and DJ822741 Benson. “Sliding interfaces with contactimpact in large-scale Lagrangian computations”. In: Computer methods in applied mechanics and engineering 51.1-3 (1985), pp. 107–137. [28] Tod A Laursen. Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis. Springer Science & Business Media, 2003. [29] Peter Wriggers and Tod A Laursen. Computational contact mechanics. Vol. 2. Springer, 2006. [30] TA Laursen and JC1242927 Simo. “A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems”. In: International Journal for numerical methods in engineering 36.20 (1993), pp. 3451– 3485. [31] Charbel Farhat and Michael Lesoinne. “Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems”. In: Computer methods in applied mechanics and engineering 182.3-4 (2000), pp. 499–515. [32] AA Lubrecht and E Ioannides. “A fast solution of the dry contact problem and the associated sub-surface stress field, using multilevel techniques”. In: (1991). [33] Richard S Sutton et al. “Fast gradient-descent methods for temporal-difference learning with linear function approximation”. In: Proceedings of the 26th annual international conference on machine learning. 2009, pp. 993–1000. [34] Seyedmeysam Khaleghian, Anahita Emami, and Saied Taheri. “A technical survey on tire-road friction estimation”. In: Friction 5 (2017), pp. 123–146. [35] Van C Mow and Rik Huiskes. Basic orthopaedic biomechanics & mechano-biology. Lippincott Williams & Wilkins, 2005. [36] Alban de Vaucorbeil and Vinh Phu Nguyen. “Modelling contacts with a total Lagrangian material point method”. In: Computer Methods in Applied Mechanics and Engineering 373 (2021), p. 113503. 113 [37] Ryan L Truby and Jennifer A Lewis. “Printing soft matter in three dimensions”. In: Nature 540.7633 (2016), pp. 371–378. [38] Metin Sitti and Hideki Hashimoto. “Teleoperated touch feedback from the surfaces at the nanoscale: modeling and experiments”. In: IEEE/ASME transactions on mechatronics 8.2 (2003), pp. 287–298. [39] Atsuya Oishi and Genki Yagawa. “A surface-to-surface contact search method enhanced by deep learning”. In: Computational Mechanics 65.4 (2020), pp. 1125–1147. [40] Kalle Kalliorinne et al. “Artificial neural network architecture for prediction of contact mechanical response”. In: Frontiers in Mechanical Engineering (2021), p. 105. [41] Daniel A Tortorelli. “Sensitivity analysis for non-linear constrained elastostatic systems”. In: International journal for numerical methods in engineering 33.8 (1992), pp. 1643–1660. [42] Alexander Konyukhov and Karl Schweizerhof. “On the solvability of closest point projection procedures in contact analysis: Analysis and solution strategy for surfaces of arbitrary geometry”. In: Computer Methods in Applied Mechanics and Engineering 197.33-40 (2008), pp. 3045–3056. [43] JC Simo and TJR Hughes. Elastoplasticity and Viscoplasticity: computational aspects. Stanford Univ., Division of Applied Mechanics, 1988. [44] J-H Heegaard and A Curnier. “An augmented Lagrangian method for discrete largeslip contact problems”. In: International Journal for Numerical Methods in Engineering 36.4 (1993), pp. 569–593. [45] James Law et al. “Speech and language therapy interventions for children with primary speech and language delay or disorder”. In: Cochrane Database of Systematic Reviews 2015.5 (1996). [46] Fehmi Cirak, Michael Ortiz, and Peter Schr¨oder. “Subdivision surfaces: a new paradigm for thin-shell finite-element analysis”. In: International Journal for Numerical Methods in Engineering 47.12 (2000), pp. 2039–2072. [47] Fehmi Cirak and Michael Ortiz. “Fully C1-conforming subdivision elements for finite deformation thin-shell analysis”. In: International Journal for Numerical Methods in Engineering 51.7 (2001), pp. 813–833. [48] Juan C Simo and David D Fox. “On a stress resultant geometrically exact shell model. Part I: Formulation and optimal parametrization”. In: Computer Methods in Applied Mechanics and Engineering 72.3 (1989), pp. 267–304. 114 [49] Paul M Naghdi. The Theory of Shells and Plates In Handbuch der Physik, Vol. VIa/2 C. Truesdell. 1972. [50] Razvan C Fetecau et al. “Nonsmooth Lagrangian mechanics and variational collision integrators”. In: SIAM Journal on Applied Dynamical Systems 2.3 (2003), pp. 381– 416. [51] Afra Zomorodian and Herbert Edelsbrunner. “Fast software for box intersections”. In: Proceedings of the sixteenth annual symposium on Computational geometry. 2000, pp. 129–138. [52] Nathan M Newmark. “A method of computation for structural dynamics”. In: Journal of the engineering mechanics division 85.3 (1959), pp. 67–94. [53] Fehmi Cirak and Julian C Cummings. “Generic programming techniques for parallelizing and extending procedural finite element programs”. In: Engineering with Computers 24.1 (2008), pp. 1–16. [54] XW Zhang, Z Tao, and QM Zhang. Dynamic behaviors of visco-elastic thin-walled spherical shells impact onto a rigid plate. 2014. [55] Chuan-yu Wu, Long-yuan Li, and Colin Thornton. “Rebound behaviour of spheres for plastic impacts”. In: International journal of impact engineering 28.9 (2003), pp. 929– 946. [56] JPA Tillett. “A study of the impact of spheres on plates”. In: Proceedings of the Physical Society. Section B 67.9 (1954), p. 677. [57] James Guilkey, Robert Lander, and Linda Bonnell. “A hybrid penalty and grid based contact method for the Material Point Method”. In: Computer Methods in Applied Mechanics and Engineering 379 (2021), p. 113739. [58] Francisco Armero and Eva Pet˝ocz. “Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems”. In: Computer methods in applied mechanics and engineering 158.3-4 (1998), pp. 269–300. [59] RH Bao and TX Yu. “Collision and rebound of ping pong balls on a rigid target”. In: Materials & Design 87 (2015), pp. 278–286. [60] Alfredo Gay Neto, Paulo de Mattos Pimenta, and Peter Wriggers. “Contact between spheres and general surfaces”. In: Computer Methods in Applied Mechanics and Engineering 328 (2018), pp. 686–716. 115 [61] XuHai Tang, Adriana Paluszny, and Robert W Zimmerman. “An impulse-based energy tracking method for collision resolution”. In: Computer Methods in Applied Mechanics and Engineering 278 (2014), pp. 160–185. [62] Jae Hyuk Sung and Byung Man Kwak. “Large displacement dynamic analysis with frictional contact by linear complementarity formulation”. In: Computers & structures 80.11 (2002), pp. 977–988. [63] Jelle Muylle and Barry HV Topping. “Parallel contact detection strategies for cable and membrane structures”. In: Computational Science—ICCS 2002: International Conference Amsterdam, The Netherlands, April 21–24, 2002 Proceedings, Part II 2. Springer. 2002, pp. 787–796. [64] Seung Hoon Paik et al. “Parallel performance of large scale impact simulations on Linux cluster super computer”. In: Computers & structures 84.10-11 (2006), pp. 732– 741. [65] G Guillamet et al. “A parallel algorithm for unilateral contact problems”. In: Computers & Structures 271 (2022), p. 106862. [66] K. Karagiozis et al. “A computational study of supersonic disk-gap-band parachutes using Large-Eddy Simulation coupled to a structural membrane”. In: Journal of Fluids and Structures 27.2 (Feb. 2011), pp. 175–192. [67] Daniel Z Huang et al. “Modeling, simulation and validation of supersonic parachute inflation dynamics during Mars landing”. In: AIAA Scitech 2020 Forum. 2020, p. 0313. [68] John Lingard and Matt Darley. “Simulation of parachute fluid structure interaction in supersonic flow”. In: 18th AIAA aerodynamic decelerator systems technology conference and seminar. 2005, p. 1607. [69] Tod Alan Laursen. “Formulation and treatment of frictional contact problems using finite elements”. PhD thesis. Stanford University, 1992. [70] Sheng Ping Wang and Eiiji Nakamachi. “The inside–outside contact search algorithm for finite element analysis”. In: International Journal for Numerical Methods in Engineering 40.19 (1997), pp. 3665–3685. [71] Qingquan Wang, Hang Yu, and Carlos Pantano. “Structural simulations of conical ribbon parachute deployment with a new detailed fabric contact algorithm, reefing line interactions and suspension line dynamics”. In: 26th AIAA Aerodynamic Decelerator Systems Technology Conference. 2022, p. 2702. 116 [72] Hang Yu, Qingquan Wang, and Carlos Pantano. “Inflation Simulations of a Conical Ribbon Parachutes at Subsonic Conditions”. In: 26th AIAA Aerodynamic Decelerator Systems Technology Conference. 2022, p. 2750. 117

Asset Metadata

Creator Wang, Qingquan (author)

Core Title A high-fidelity geometric contact method for thin fabrics using shell finite elements in explicit and implicit time scheme

Contributor Electronically uploaded by the author (provenance)

School Andrew and Erna Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Degree Conferral Date 2024-08

Publication Date 09/02/2024

Defense Date 08/05/2024

Publisher Los Angeles, California (original), University of Southern California (original), University of Southern California. Libraries (digital)

Tag contact detection,exact contact constraint enforcement,parallel contact,thin surface contact.

Format theses (aat)

Language English

Advisor Pantano-Rubino, Carlos (committee chair), Garikipati, Krishna (committee member), Ghanem, Roger (committee member), Oberai, Assad (committee member), Plucinsky, Paul (committee member)

Creator Email qingquan@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC11399A6H4

Unique identifier UC11399A6H4

Identifier etd-WangQingqu-13467.pdf (filename)

Legacy Identifier etd-WangQingqu-13467

Document Type Dissertation

Format theses (aat)

Rights Wang, Qingquan

Internet Media Type application/pdf

Type texts

Source 20240903-usctheses-batch-1206 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email uscdl@usc.edu

Abstract (if available)

Abstract While numerical models for contact mechanics are increasingly common, achieving efficient geometric contact resolution with a robust search strategy remains challenging. This research addresses this issue by introducing a comprehensive solution that includes an exact geometric contact mechanics algorithm for thin shell finite elements, applicable to both explicit and implicit time schemes.

In the explicit time scheme, the method offers precise geometric resolution of self-contact interactions through a sub-time-step marching technique, employs adaptive data structures to minimize computational overhead, and features a specialized parallelization approach with load-balancing capabilities. To manage the polynomial time complexity, an efficient detection algorithm divides the problem into global and local contact detection phases. Impact equations are then used to resolve contact events while preserving kinematic energy and momentum. This algorithm is integrated with an MPI-based parallelization framework for thin-shell finite element solvers, ensuring balanced load distribution.

The research also investigates the integration of the two-phase contact resolution approach within the implicit time scheme. Here, the use of a signed gap function provides precise contact constraints, offering an alternative solution for contact-intensive cases where explicit methods may be time-consuming due to sub-time step marching.

The robustness and accuracy of the algorithm are validated through three numerical studies, and strong scaling analysis demonstrates the scalability of the parallelization approach.