Hugo Verhelst

My work focusses on numerical methods for solid and fluid mechanics, using the mathematical concept of Isogeometric Analysis (IGA). In particular, I contribute to the fields of adaptive meshing for IGA, phase-field fracture and data-driven analysis, and I am a main developer of the Geometry+Simulation Modules (G+Smo).

Wrinkles are ubiquitous in the world around us. In our daily lives, we encounter wrinkles in various forms, whether in our clothes or on our skin. Wrinkles emerge as a result of a delicate interplay between bending, membrane, and foundation stiffness contributions within membranes. While experimental investigations provide insights into the physics underlying wrinkling, numerical investigations find their purpose in the design, analysis, and optimisation of membranes subjected to wrinkling. Nevertheless, the numerical simulation of membrane wrinkling presents several challenges. Firstly, wrinkling constitutes a buckling phenomenon in membranes with low bending stiffness. Wrinkles have the potential to evolve into folds, creases, or other wrinkling patterns as loads or displacements increase. Secondly, the wavelengths of wrinkling can be orders of magnitude smaller than the overall geometry, requiring a small resolution of the numerical simulation and hence increasing computational costs. Overall, the question arises of how to design robust and accurate numerical models for the analysis of wrinkled membranes. This dissertation is subdivided into four parts and aims to provide answers to this question. The first theme considers hyperelastic material modelling, with a focus on developing wrinkling models under large strains. The shell model employed in this dissertation is based on the isogeometric analysis paradigm. Specifically, the Kirchhoff–Love shell model is used, which leverages the higher-order continuity of underlying spline spaces. Chapter 3 extends hyperelastic material formulations to stretch-based materials, enabling the use of the isogeometric analysis paradigm for rubber-like shells. Since the modelling of wrinkling patterns imposes physical scales limiting element mesh sizes, chapter 4 introduces a hyperelastic isogeometric membrane element that incorporates an implicit wrinkling model, thus avoiding explicit modelling of wrinkling amplitudes. The second theme addresses adaptive methods. On the one hand, spatial adaptivity enhances the local detail in a numerical simulation. Chapter 5 presents an adaptive isogeometric analysis framework based on intuitive goal functions, such as wrinkling amplitudes, to guide adaptive meshing routines. On the other hand, temporal or quasi-temporal adaptivity serves to enhance the efficiency of dynamic or quasi-static simulations. Chapter 6 introduces an adaptive parallel arc-length method. The method’s adaptivity arises as a by-product of parallelisation efforts aimed at reducing computational times for quasi-static simulations. The advantage of the smoothness inherent in the spline spaces used in isogeometric analysis is limited to simple topologies. To benefit from this smoothness in complex geometries, the third theme of this dissertation focuses on complex domain modelling. Chapter 7 presents a qualitative and quantitative comparison of unstructured spline constructions for multi-patch modelling using isogeometric analysis. This chapter offers insights and suggestions for future developments related to unstructured spline constructions. The final theme of this dissertation concerns the reproducibility of the developed methods. In this section, design considerations are presented for an open-source software library, along with small examples, aimed at ensuring easy reproducibility and supporting future research in the three themes mentioned earlier. In summary, this dissertation offers a wide range of methods for the isogeometric analysis of structural instabilities in thin-walled structures, including the modelling of wrinkling. The concepts developed in terms of hyperelasticity expand the applicability of wrinkling models to encompass large strains. The concepts developed in terms of adaptivity provide intuitive error estimators that drive local refinement in space, as well as a novel continuation method that eliminates the inherently serial arc-length methods. Through the use of unstructured splines, complex domains become accessible for the analysis of structural stabilities. By creating an open-source, forward-compatible software library, these concepts are made available for future developments in the field of isogeometric analysis of wrinkling.

In recent decades, the study of fracture propagation in solids has increasingly relied on phase-field models. Several recent contributions have highlighted the potential of this approach in both static and dynamic frameworks. However, a major limitation remains the high computational cost. Two main strategies have been identified to mitigate this issue: the use of locally refined meshes and the adoption of higher-order models. In this work, leveraging Truncated Hierarchical B-splines (THB-splines), we introduce adaptive simulations of higher-order phase-field formulations (AT1 and AT2), focusing primarily on two-dimensional fracture problems.

Isogeometric Analysis (IGA) bridges Computer-Aided Design (CAD) and Finite Element Analysis (FEA) by employing splines as a common basis for geometry and analysis. One of the advantages of IGA is in the realm of thin shell analysis: due to the arbitrary continuity of the spline basis, Kirchhoff–Love shells can be modeled without the need to introduce unknowns for the mid-plane rotations, leading to a reduction in the number of unknowns. In this paper, we provide the background of an implementation of Isogeometric Kirchhoff–Love shells within the Geometry + Simulation Modules (G+Smo). This paper accompanies multiple previous publications and elaborates on the design of the software used in these papers, rather than the novelty of the methods presented therein. The presented implementation provides patch coupling via penalty methods and unstructured splines, goal-oriented error estimators, several algorithms for structural analysis and advanced algorithms for the modeling of wrinkling in hyperelastic membranes. These methods are all contained in three new modules in G+Smo: a module for Kirchhoff–Love shells, a module for structural analysis, and a module for unstructured spline constructions. As motivated in this paper, the modules are implemented to be compatible with future developments. For example, by providing base implementations of material laws, by using black-box functions for the structural analysis module, or by providing a standardized approach for the implementation of unstructured spline constructions. Overall, this paper demonstrates that the new modules contribute to a versatile ecosystem for the modeling of multi-patch shell problems through fast off-the-shelf solvers with a simple interface, designed to be extended in future research.

This paper presents spline-based coupling methods for partitioned multiphysics simulations, specifically designed for isogeometric analysis (IGA) based solvers. Traditional vertex-based coupling approaches face significant challenges when applied to IGA solvers, including geometric accuracy issues, interpolation errors, and substantial communication overhead. The methodology draws on the IGA mathematical framework to deliver coupling solutions that preserve high-order continuity and exact geometric representation of splines. We develop two complementary strategies: (1) a spline-vertex coupling method enabling efficient interaction between IGA and conventional solvers, and (2) a fully isogeometric coupling approach maximizing accuracy for IGA-to-IGA communication. Both theoretical analysis and extensive numerical experiments demonstrate that our spline-based methods significantly reduce communication overhead compared to traditional approaches while enhancing geometric accuracy through exact boundary representation and maintaining higher-order solution continuity across coupled interfaces. We quantitatively confirm communication efficiency benefits through systematic measurements of transfer times and data volumes across various mesh refinement levels. Our benchmark studies demonstrate geometric fidelity advantages while highlighting how splines naturally preserve solution derivatives across interfaces without requiring additional computation. This work provides efficient coupling strategies tailored to IGA-based solvers and establishes a practical bridge between IGA and traditional discretization methods, enabling broader adoption of IGA in established simulation workflows.

Phase separation leads to evolving interfaces that require sufficient spatial resolution to accurately capture their dynamics. We present a multi-dimensional higher-order adaptive isogeometric analysis framework for phase-separation problems, based on a phase-field formulation of the Cahn–Hilliard equation. As basis functions, we employ Truncated Hierarchical B-splines, which form a partition of unity and enable local mesh refinement and coarsening. The adaptive meshing scheme refines the mesh at interfaces and coarsens it in the bulk, with the mesh resolution evolving alongside the solution. Element marking is guided by the solution field, which identifies interface locations, and solution transfer between successive meshes is performed via a quasi-interpolation operator that is naturally parallelizable and efficiently reduces computational cost. The performance of the framework is demonstrated through a spatial convergence study and a series of 2D and 3D numerical examples, showing that locally adaptive meshes accurately track evolving interfaces while reducing the computational effort per iteration compared to uniform tensor-product discretizations.

Wrinkling is the phenomenon of out-of-plane deformation patterns in thin walled structures, as a result of a local compressive (internal) loads in combination with a large membrane stiffness and a small but non-zero bending stiffness. Numerical modelling typically involves thin shell formulations. As the mesh resolution depends on the wrinkle wave lengths, the analysis can become computationally expensive for shorter ones. Implicitly modelling the wrinkles using a modified kinematic or constitutive relationship based on a taut, slack or wrinkled state derived from a so-called tension field, a simplification is introduced in order to reduce computational efforts. However, this model was restricted to linear elastic material models in previous works. Aiming to develop an implicit isogeometric wrinkling model for large strain and hyperelastic material applications, a modified deformation gradient has been assumed, which can be used for any strain energy density formulation. The model is an extension of a previously published model for linear elastic material behaviour and is generalised to other types of discretisation as well. The extension for hyperelastic materials requires the derivative of the material tensor, which can be computed numerically or derived analytically. The presented model relies on a combination of dynamic relaxation and a Newton–Raphson solver, because of divergence in early Newton–Raphson iterations as a result of a changing tension field, which is not included in the stress tensor variation. Using four benchmarks, the model performance is evaluated. Convergence with the expected order for Newton–Raphson iterations has been observed, provided a fixed tension field. The model accurately approximates the mean surface of a wrinkled membrane with a reduced number of degrees of freedom in comparison to a shell solution.

In order to perform isogeometric analysis with increased smoothness on complex domains, trimming, variational coupling or unstructured spline methods can be used. The latter two classes of methods require a multi-patch segmentation of the domain, and provide continuous bases along patch interfaces. In the context of shell modelling, variational methods are widely used, whereas the application of unstructured spline methods on shell problems is rather scarce. In this paper, we therefore provide a qualitative and a quantitative comparison of a selection of unstructured spline constructions, in particular the D-Patch, Almost-C¹, Analysis-Suitable G¹ and the Approximate C¹ constructions. Using this comparison, we aim to provide insight into the selection of methods for practical problems, as well as directions for future research. In the qualitative comparison, the properties of each method are evaluated and compared. In the quantitative comparison, a selection of numerical examples is used to highlight different advantages and disadvantages of each method. In the latter, comparison with weak coupling methods such as Nitsche’s method or penalty methods is made as well. In brief, it is concluded that the Approximate C¹ and Analysis-Suitable G¹ converge optimally in the analysis of a bi-harmonic problem, without the need of special refinement procedures. Furthermore, these methods provide accurate stress fields. On the other hand, the Almost-C¹ and D-Patch provide relatively easy construction on complex geometries. The Almost-C¹ method does not have limitations on the valence of boundary vertices, unlike the D-Patch, but is only applicable to biquadratic local bases. Following from these conclusions, future research directions are proposed, for example towards making the Approximate C¹ and Analysis-Suitable G¹ applicable to more complex geometries.

Parallel computing is omnipresent in today’s scientific computer landscape, starting at multicore processors in desktop computers up to massively parallel clusters. While domain decomposition methods have a long tradition in computational mechanics to decompose spatial problems into multiple subproblems that can be solved in parallel, advancing solution schemes for dynamics or quasi-statics are inherently serial processes. For quasi-static simulations, however, there is no accumulating ’time’ discretization error, hence an alternative approach is required. In this paper, we present an Adaptive Parallel Arc-Length Method (APALM). By using a domain parametrization of the arc-length instead of time, the multi-level error for the arc-length parametrization is formed by the load parameter and the solution norm. Given coarse approximations of arc-length intervals, finer corrections enable the parallelization of the presented method. This results in an arc-length method that is parallel within a branch and inherently adaptive. This concept is easily extended for bifurcation problems. The performance of the method is demonstrated using isogeometric Kirchhoff-Love shells on problems with snap-through and pitch-fork instabilities and applied to the problem of a snapping meta-material. These results show that parallel corrections are performed in a fraction of the time of the serial initialization, achievable on desktop scale.

Mesh adaptivity is a technique to provide detail in numerical solutions without the need to refine the mesh over the whole domain. Mesh adaptivity in isogeometric analysis can be driven by Truncated Hierarchical B-splines (THB-splines) which add degrees of freedom locally based on finer B-spline bases. Labeling of elements for refinement is typically done using residual-based error estimators. In this paper, an adaptive meshing workflow for isogeometric Kirchhoff–Love shell analysis is developed. This framework includes THB-splines, mesh admissibility for combined refinement and coarsening and the Dual-Weighted Residual (DWR) method for computing element-wise error contributions. The DWR can be used in several structural analysis problems, allowing the user to specify a goal quantity of interest which is used to mark elements and refine the mesh. This goal functional can involve, for example, displacements, stresses, eigenfrequencies etc. The proposed framework is evaluated through a set of different benchmark problems, including modal analysis, buckling analysis and non-linear snap-through and bifurcation problems, showing high accuracy of the DWR estimator and efficient allocation of degrees of freedom for advanced shell computations.