DC1 CFD techniques for IBRA-type discretizations:

Explore the use of IBRA-type discretizations in the context of CFD.


The Shifted Boundary Method in IGA

    IGA has transformed computational mechanics by bridging the gap between Computer-Aided Design (CAD) and Computer-Aided Engineering (CAE). Originally introduced by Hughes et al. [1,2,3,4,5], IGA enables precise geometric representation and high continuity across element boundaries [6], making it an effective tool for simulating complex geometries [7,8,9]. The foundation of IGA lies in B-Spline and NURBS basis functions, which allow for smooth transitions and local refinement, enhancing the accuracy and robustness of simulations [10,11].
    Despite its advantages, traditional boundary-fitted IGA methods face challenges with complex geometries, such as the need for watertight geometries and high computational costs in handling trimmed models [12,13]. Immersed boundary IGA methods, such as the Finite Cell Method (FCM) [14,15] and Isogeometric Boundary Representation Analysis (IBRA) [16,17,18,19,20], alleviate some of these challenges by using non-boundary-fitted discretizations. However, these approaches are often hindered by issues like the small cut-cell problem [21,22], which leads to poor computational performance and difficulty in solver convergence.

    The SBM [23,24], initially developed within the FEM context, offers a promising solution to these challenges. SBM simplifies integration over complex geometries by shifting boundary conditions to a surrogate boundary and modifying boundary values using Taylor expansions. This approach avoids the small cut-cell problem, maintains optimal accuracy, and enables significant simplification in mesh generation and refinement. Previous applications of SBM in FEM have demonstrated its effectiveness in areas such as elasticity and incompressible fluid dynamics [25,26,27,28]. This study pioneers the integration of SBM with IGA, proposing an alternative to IBRA that eliminates the need for trimming techniques. By treating the boundary as unfitted, the proposed approach simplifies pre-processing and retains the accuracy and continuity advantages of IGA. The structure includes an introduction to the SBM framework, its application to the Poisson problem, and numerical validations.[31].

    The Shifted Boundary Method (SBM) introduces key innovations.

    • Surrogate boundary representation: the true boundary, denoted as Γ, is replaced by a surrogate boundary, which is constructed using the edges of a Cartesian background grid. This simplification avoids complications associated with small cut-cells in boundary-fitted methods (see Figure 1).
    • Taylor expansion for boundary conditions: boundary conditions are imposed on the surrogate boundary using a Taylor expansion of arbitrary order mmm. This approach accurately approximates the conditions that would have been applied on the true boundary, preserving the desired level of accuracy.
    • Weak imposition of boundary conditions: the shifted boundary operator is employed to weakly enforce Dirichlet and Neumann boundary conditions, ensuring numerical stability and compatibility with the unfitted nature of the method.

    These features collectively enable SBM to achieve robust and accurate solutions without the computational overhead associated with traditional boundary-fitted approaches.

    Figure 1: An example of a domain 𝛺 with a ring-like shape treated with the SBM in a Cartesian grid. Blue and red solid lines denote the true and surrogate boundaries respectively. The computational domain is shaded in light yellow, while brown denotes the intersected cells that are not part of the SBM computation.

    The Taylor expansion plays a fundamental role in the SBM for accurately shifting boundary conditions from the true boundary to the surrogate boundary. This process involves performing an m-th order Taylor expansion of the variable of interest, u, at the surrogate boundary. Assuming u is sufficiently smooth within the strip between Γ and the surrogate​, the expansion incorporates directional derivatives along the vector d, capturing the influence of the distance between the true and surrogate boundaries.

    The m-th order Taylor expansion includes a remainder term that becomes negligible as the distance ∥d∥ approaches zero. Dirichlet conditions defined on Γ are extended​ using the boundary shift operator, enabling the surrogate boundary to approximate the true boundary conditions. 

    Neglecting the remainder term simplifies the boundary condition approximation to an m-th order shifted expression, which is subsequently enforced weakly in the SBM framework.

    This approach ensures both computational efficiency and numerical accuracy while avoiding the complexities associated with direct boundary alignment [31].

    Moreover, we have taken advantage of two ideas recently published:

    1. Penalty-Free formulation: A penalty-free weak formulation [29] eliminates the need for penalty calibration, which is traditionally required for Nitsche’s method. The SBM variational form ensures numerical stability and accuracy without the overhead of parameter tuning.
    2. Optimal surrogate boundary selection: the surrogate boundary minimizes the distance function between the surrogate and true boundaries [30]. A parameter λ determines the inclusion of cut elements in the surrogate domain, with λ=0.5 found to be optimal for minimizing error.

    Figure 2: Selection of the surrogate boundary. Two-dimensional complex embedded holes treated using the SBM. The red solid line denotes the surrogate boundary corresponding to 𝜆 = 0 and the green one to 𝜆 = 0.5 (optimal surrogate boundary). 

    In Figure 3 we report the convergence studies of the two cases in Figure 2 when applying only Dirichlet boundary conditions.

    Figure 3: Selection of the surrogate boundary. Convergence study of two cases curvilinear and diamond of Figure 2  respectively, considering 𝜆 = 0 (red solid line) and 𝜆 = 0.5 (blue solid line). The order of the basis functions is indicated by the following symbols: triangles 𝑝 = 1, squares 𝑝 = 2, and stars 𝑝 = 3.

    As previously mentioned, B-Spline basis functions are utilized to represent both the geometry and the solution fields. The implementation of the SBM avoids the use of trimmed knot spans, which are characteristic of the IBRA approach (refer to the left panel of Figure 4). The right panel of Figure 4 presents a comparison between the SBM and IBRA methodologies using linear, quadratic, and cubic B-Spline basis functions. The results demonstrate that both methods achieve very similar error norms and exhibit optimal convergence rates across all orders of the basis functions.

    Figure 4: SBM vs IBRA. Left: the diamond shape applying the IBRA technique; the green dots denote the integration points, while the solid blue line highlights the true boundary. Right: the convergence analysis of the same case, using 𝑝 = 1 (triangular symbols), 𝑝 = 2 (square symbols), and 𝑝 = 3 (star symbols) and comparing the IBRA approach (red solid line) and the SBM (blue solid line) with the optimal surrogate boundary.

    For more details and examples we refer to the paper “The Shifted Boundary Method in Isogeometric Analysis” [31], where the authors also analyze the case of Neumann boundary conditions.

    SBM-IGA for Stokes fluids

      This section explores the application of the SBM within the IGA framework to solve Stokes flow problems. The Stokes equations describe the behavior of viscous, incompressible fluids under low Reynolds number conditions, where inertial forces are negligible compared to viscous forces. These equations involve coupling between the velocity field and the pressure field, which introduces additional numerical complexity compared to the Poisson problem discussed in Section 1.

      The Stokes problem consists of solving the following equations:

      together with the boundary conditions. The weak form after integration by parts:

      where v is the test function for the velocity vector field, and q is the test function for the pressure scalar field.

      The Stokes problem is solved within a Variational Multi-Scale (VMS) [32-33] framework, which provides a stabilized formulation to handle the coupling of the velocity and pressure fields. In the VMS approach:

      • The solution space is decomposed into resolved (coarse) and unresolved (fine) scales.
      • Stabilization terms are introduced to account for the effects of unresolved scales on the resolved solution, enhancing numerical robustness and accuracy.

      where

      The VMS method, combined with B-Spline basis functions in the IGA framework, allows for arbitrary order discretizations of the velocity and pressure fields. Unlike the simpler Poisson problem, the Stokes equations involve higher-order terms that do not simplify directly, making the implementation of VMS more intricate. The weak formulation now is

      Therefore we have three additional terms with respect to a formulation without the sub-scale:

      The implementation of the SBM for Stokes flow in IGA follows a similar strategy as in Section 1:

      • B-Spline basis functions: The velocity and pressure fields are approximated using B-Splines of arbitrary order, ensuring smooth and accurate solutions.
      • Shifted boundary conditions: Dirichlet and Neumann boundary conditions are imposed weakly on the surrogate boundary using Taylor expansions.
      • Stabilization terms: VMS stabilization terms are added to the variational formulation to address numerical instabilities arising from the coupling between velocity and pressure fields.

      The combination of SBM and VMS within the IGA framework enables the efficient and accurate solution of Stokes flow problems.

      In addition to Stokes fluids, we introduced the use of  non-Newtonian fluids, specifically the Bingham fluid model. The Bingham constitutive law describes a material that behaves as a solid when the stress magnitude is below a certain yield stress and flows as a viscous fluid when the stress exceeds the yield stress. The stabilized formulation of the Bingham fluid incorporates an exponential factor that smooths the transition.

      For small values of the regularization parameter, the behavior transitions smoothly between solid-like and fluid-like phases. Figure 5 shows the Stress vs Strain diagram.

      Figure 5: Stress vs. Strain Rate for the Bingham Model with Various Regularization Parameters (m). The graph illustrates the stress response of a Bingham fluid as a function of strain rate (gamma_dot​) for different values of the regularization parameter mmm.

      We can show some preliminary results obtained with a manufactured solution for a Stokes fluid having a Bingham constitutive law. The manufactured solution is a smooth non-polynomial function in both the components of the velocity and in the pressure field. The convergence studies in the L2 norm of the error for velocity and pressure is shown in Figure 6 and Figure 7 (m = 10 and m = 300 respectively) for linear, quadratic, and cubic B-Splines basis functions.

      Figure 6: Convergence Study for Velocity and Pressure Errors with Yield Stress τ_0 = 100 and regularization parameter m=10. The figure shows the error convergence for the velocity (left) and pressure (right) fields as the mesh size h decreases. Results are presented for B-Spline basis functions of degree p=1 (blue), p=2 (red), and p=3 (cyan).

      Figure 7: Convergence Study for Velocity and Pressure Errors with Yield Stress τ_0 = 100 and regularization parameter m=300. The figure shows the error convergence for the velocity (left) and pressure (right) fields as the mesh size h decreases. Results are presented for B-Spline basis functions of degree p=1 (blue), p=2 (red), and p=3 (cyan).

      We observed that the exact solution must be smooth, including its derivatives, to achieve optimal convergence when using B-Splines. Unlike classical FEM, where the solution is typically continuous between elements, the IGA framework with B-Splines ensures continuity throughout the domain. If the solution exhibits discontinuities, particularly under nonlinear constitutive laws such as the Bingham model, the method fails to achieve the expected optimal convergence rates. Therefore, when applying IGA to problems with nonlinear laws, it is essential to ensure the smoothness of the solution to fully leverage the advantages of B-Splines.

      We have also validated our code by testing it on the well-known lid-driven cavity problem, as illustrated in Figure 8, adapted from [34]. This benchmark problem is widely recognized in computational fluid dynamics as a standard test for assessing the accuracy and robustness of numerical methods applied to fluid flow simulations. The lid-driven cavity problem involves a square domain where the top boundary (lid) moves at a constant velocity, inducing circulation within the fluid. The remaining boundaries are stationary and impose no-slip conditions. This configuration creates a primary vortex at the center of the domain.

      Figure 8: Cavity problem. The flow in a 2D cavity filled with a Bingham fluid. The problem is set up following Mitsoulis and Zisis [34]. Defining a square domain Ω = (0, H) × (0, H), we impose a horizontal velocity u_x = 1 m/s on the y = H.

      Since no analytical solution is available for the lid-driven cavity problem, we have compared our results with the benchmark studies presented in [34] and [35]. This comparison focuses on two critical aspects of the flow behavior for different values of the yield stress: the extent of the yielded and unyielded regions within the domain, and the position of the eye of the main vortex that forms in the central-upper zone.
      The yielded region, defined by areas where the fluid behaves as a viscous material, and the unyielded region, where the fluid behaves like a solid due to insufficient stress to induce flow, are key features of viscoplastic flows. Accurately capturing these regions provides insight into the robustness and precision of the numerical method (see Figure 9 and Figure 10). 

      Figure 9: Yielded region. Comparison of the yielded region in the cavity flow problem. The Bingham number is equal to 20.0 and m = 300. Left: results with B-Splines with p=2 and 50×50 knot spans. Right: results in FEM from [35] using local refinement.

      Figure 10: Yielded region. Comparison of the yielded region in the cavity flow problem. The Bingham number is equal to 20.0 and m = 300. Left: results with B-Splines with p=2 and 100×100 knot spans. Right: results in FEM from [35] using local refinement.

      Additionally, the location of the vortex eye serves as an important indicator of the method’s ability to predict the flow dynamics accurately. The vortex position is influenced by both the yield stress and the imposed lid motion, making it a sensitive measure for evaluating the performance of our implementation (see Figure 11).

      Figure 11: Eye vertex position. Comparison of the eye vertex y-coordinate in the cavity flow problem varying the Bingham number, m = 300. Left: results with B-Splines with p=1 & p=2. Right: results in FEM from [34] with 40×40 divisions.

      Our results exhibit strong agreement with those in [34] and [35], confirming the capability of our Shifted Boundary Method (SBM) within the IsoGeometric Analysis (IGA) framework to handle complex boundary conditions and accurately simulate viscoplastic flow behavior in this classical benchmark problem.

      REFERENCES

        1. Thomas J.R. Hughes, John A. Cottrell, Yuri Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (39) (2005) 4135–4195.
        2. John A. Cottrell, Alessandro Reali, Yuri Bazilevs, Thomas J.R. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Engrg. 195 (41) (2006) 5257–5296.
        3. Yuri Bazilevs, Lourenco Beirao da Veiga, John A. Cottrell, Thomas J.R. Hughes, Giancarlo Sangalli, Isogeometric analysis: Approximation, stability and error estimates for h-refined meshes, Math. Models Methods Appl. Sci. 16 (07) (2006) 1031–1090.
        4. Yuri Bazilevs, Victor M. Calo, Yongjie Zhang, Thomas J.R. Hughes, Isogeometric fluid–structure interaction analysis with applications to arterial blood flow, Comput. Mech. 38 (4) (2006) 310–322.
        5. Yongjie Zhang, Yuri Bazilevs, Samrat Goswami, Chandrajit L. Bajaj, Thomas J.R. Hughes, Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow, Comput. Methods Appl. Mech. Engrg. 196 (29) (2007) 2943–2959.
        6. John A. Cottrell, Thomas J.R. Hughes, Alessandro Reali, Studies of refinement and continuity in isogeometric structural analysis, Comput. Methods Appl. Mech. Engrg. 196 (41) (2007) 4160–4183.
        7. Josef Kiendl, Kai-Uwe Bletzinger, Johannes Linhard, Roland Wüchner, Isogeometric shell analysis with Kirchhoff–Love elements, Comput. Methods Appl. Mech. Engrg. 198 (49) (2009) 3902–3914.
        8. David J. Benson, Yuri Bazilevs, Ming-Chen Hsu, Thomas J.R. Hughes, Isogeometric shell analysis: The Reissner–Mindlin shell, Comput. Methods Appl. Mech. Engrg. 199 (5) (2010) 276–289, Computational Geometry and Analysis.
        9. Kenji Takizawa, Yuri Bazilevs, Tayfun E. Tezduyar, Ming-Chen Hsu, Takuya Terahara, Computational cardiovascular medicine with isogeometric analysis, J. Adv. Eng. Comput. 6 (3) (2022) 167–199.
        10. Massimo Carraturo, Carlotta Giannelli, Alessandro Reali, Rafael Vázquez, Suitably graded THB-spline refinement and coarsening: Towards an adaptive isogeometric analysis of additive manufacturing processes, Comput. Methods Appl. Mech. Engrg. 348 (2019) 660–679.
        11.  Carlotta Giannelli, Bert Jüttler, Hendrik Speleers, THB-splines: The truncated basis for hierarchical splines, Comput. Aided Geom. Design 29 (7) (2012) 485–498, Geometric Modeling and Processing 2012.
        12. Dominik Schillinger, Luca Dedè, Michael A. Scott, John A. Evans, Michael J. Borden, Ernst Rank, Thomas J.R. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 116–150, Higher Order Finite Element and Isogeometric Methods.
        13. Martin Ruess, Dominik Schillinger, Yuri Bazilevs, Vasco Varduhn, Ernst Rank, Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method, Internat. J. Numer. Methods Engrg. 95 (10) (2013) 811–846.
        14. Jamshid Parvizian, Alexander Düster, Ernst Rank, Finite cell method: h-and p-extension for embedded domain problems in solid mechanics, Comput. Mech. 41 (1) (2007) 121–133.
        15. Ernst Rank, Martin Ruess, Stefan Kollmannsberger, Dominik Schillinger, Alexander Düster, Geometric modeling, isogeometric analysis and the finite cell method, Comput. Methods Appl. Mech. Engrg. 249–252 (2012) 104–115.
        16. Michael Breitenberger, Andreas Apostolatos, Philipp Bucher, Roland Wüchner, Kai-Uwe Bletzinger, Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures, Comput. Methods Appl. Mech. Engrg. 284 (2015) 401–457, Isogeometric Analysis Special Issue.
        17. Tobias Teschemacher, Anna M. Bauer, Thomas Oberbichler, Micheal Breitenberger, Riccardo Rossi, Roland Wüchner, Kai-Uwe Bletzinger, Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis, Adv. Model. Simul. Eng. Sci. 5 (2018) 1–54.
        18. Tobias Teschemacher, Anna M. Bauer, Ricky Aristio, Manuel Meßmer, Roland Wüchner, Kai-Uwe Bletzinger, Concepts of data collection for the CAD-integrated isogeometric analysis, Eng. Comput. 38 (6) (2022) 5675–5693.
        19. Manuel Meßmer, Tobias Teschemacher, Lukas F. Leidinger, Roland Wüchner, Kai-Uwe Bletzinger, Efficient CAD-integrated isogeometric analysis of trimmed solids, Comput. Methods Appl. Mech. Engrg. 400 (2022) 115584.
        20. Manuel Meßmer, Stefan Kollmannsberger, Roland Wüchner, Kai-Uwe Bletzinger, Robust numerical integration of embedded solids described in boundary representation, Comput. Methods Appl. Mech. Engrg. 419 (2024) 116670.
        21. Erik Burman, Ghost penalty, C. R. Math. 348 (21–22) (2010) 1217–1220.
        22. Santiago Badia, Eric Neiva, Francesc Verdugo, Linking ghost penalty and aggregated unfitted methods, Comput. Methods Appl. Mech. Engrg. 388 (2022) 114232.
        23. Alex Main, Guglielmo Scovazzi, The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems, J. Comput. Phys. 372 (2018) 972–995.
        24. Alex Main, Guglielmo Scovazzi, The shifted boundary method for embedded domain computations. Part II: Linear advection–diffusion and incompressible Navier–Stokes equations, J. Comput. Phys. 372 (2018) 996–1026.
        25. Nabil M. Atallah, Guglielmo Scovazzi, Nonlinear elasticity with the shifted boundary method, Comput. Methods Appl. Mech. Engrg. 426 (2024) 116988.
        1. Efthymios N. Karatzas, Giovanni Stabile, Leo Nouveau, Guglielmo Scovazzi, Gianluigi Rozza, A reduced-order shifted boundary method for parametrized incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 370 (2020) 113273.
        2. Nabil M. Atallah, Claudio Canuto, Guglielmo Scovazzi, The second-generation shifted boundary method and its numerical analysis, Comput. Methods Appl. Mech. Engrg. 372 (2020) 113341.
        3. Nabil M. Atallah, Claudio Canuto, Guglielmo Scovazzi, Analysis of the shifted boundary method for the Poisson problem in domains with corners, Math. Comp. 90 (331) (2021) 2041–2069.
        4. Haydel Collins, Alexei Lozinski, Guglielmo Scovazzi, A penalty-free shifted boundary method of arbitrary order, Comput. Methods Appl. Mech. Engrg. 417 (2023) 116301, A Special Issue in Honor of the Lifetime Achievements of T. J. R. Hughes.
        5. Cheng-Hau Yang, Kumar Saurabh, Guglielmo Scovazzi, Claudio Canuto, Adarsh Krishnamurthy, Baskar Ganapathysubramanian, Optimal surrogate boundary selection and scalability studies for the shifted boundary method on octree meshes, Comput. Methods Appl. Mech. Engrg. 419 (2024) 116686.
        6. Antonelli, N., Aristio, R., Gorgi, A., Zorrilla, R., Rossi, R., Scovazzi, G., Wüchner, R. (2024). The Shifted Boundary Method in Isogeometric Analysis. Computer Methods in Applied Mechanics and Engineering, 430, 117228.
        7. Codina, R., Badia, S., Baiges, J., & Principe, J. (2018). Variational multiscale methods in computational fluid dynamics. Encyclopedia of computational mechanics, 1-28.
        8. Codina, R. (2000). Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Computer methods in applied mechanics and engineering, 190(13-14), 1579-1599.
        9. Mitsoulis E, Zisis T. Flow of Bingham plastics in a lid-driven square cavity. Journal of Non-Newtonian Fluid Mechanics 2001; 101(1–3):173 – 180.
        10. Cotela Dalmau, Jordi, Riccardo Rossi, and Antonia Larese De Tetto. “Simulation of two-and three-dimensional viscoplastic flows using adaptive mesh refinement.” International journal for numerical methods in engineering (2017).

        DC10 Implementation of IGA in the design and analysis workflow of machine elements and systems:

        Although the basic principles of IGA are known for some decades, its application in industry is less spread than one would expect considering its fundamental advantage: the representation of real surfaces in a well-defined mathematical way instead of a discrete non-smooth finite element mesh.

        Progress was made in FEM and IGA simulation studies using advanced workflows, including ANSA-Epilysis-META, ANSA-LSDYNA-META, and ANSA-Kratos-Rhino, for tensile test specimen analysis. The results validates against theoretical and experimental data.

        1 Test different workflows

        According to [Tuo21], the test specimens, shown in Fig 1(b), are fabricated from a DP900 steel sheet with a thickness of 1.5mm, a gauge length 55mm, a Young’s modulus E=207,600MPa, and a Poisson’s ratio ν=0.3. The uniaxial tensile test was conducted on an MTS universal tensile machine with a maximum load capacity of 100kN, Fig 1(a). Boundary conditions were applied such that one end of the specimen was fully fixed, while the opposite end was subjected to axial displacement at a loading rate of 2.0 mm/min. The tensile specimen exhibits complex mechanical behavior under tensile loading, including elastic deformation, yielding, plastic deformation, and fracture. The force–displacement curves is shown in Fig 1(c).

        Our study focuses on the elastic deformation stage. As shown in Fig 1(c), for displacements ≤0.04 mm and uniaxial forces ≤3kN, the response demonstrates a clear linear elastic behavior.

        (a) Experimental equipment
        (b) specimen (mm)
        (c) Force–displacement curve of uniaxial tensile test

        Figure 1: Test specimen model [1]

        Experimentally, Fig 1(c) shows when F = 3KN, the displacement is 0.0398mm. While Fig 1(c) is obtained by point-pick from the reference paper, which contains acceptable physical error.

        Theoretically, the elongation Δx and stress σ are calculated by formula:

        Δx = FL/EA

        σ = F/A

        Where F is the tensile force, L is the effective length, E is the Young’s modulus, A is the initial cross-sectional area.

        According to Fig 1(a), only the elongation of the gauge length 55 mm is considered. The theoretical elongation calculation is organized by three parts, as shown in table 1.

        Table 1: Theoretical elongation calculation

        Cross-sectional area(mm^2)Effective length (mm)Elongation Δx (mm)
        14mm*1.5mm(55-4)/2 = 25.50.0175
        5mm *1.5mm40.0077
        14mm*1.5mm(55-4)/2 = 25.50.0175
        sum/0.0427

        The nominal stress at the notch cross-section is 400 Ν/mm^2. The stress concentration factor 1.62[Pe08] should be considered, so the maximum stress at the notch is 648 Ν/mm^2.

        In the subsequent analysis, the elastic deformation of the notched tensile specimen will be simulated using different workflows.

        1.1 ANSA-Epilysis-META Workflow

        In this section, a comprehensive FEM process ANSA–Epilysis (Nastran SOL101)–Meta workflow is used to simulate the linear static analysis of the specimen. The left end of the specimen is fully fixed, and a 3 kN force is applied along the axial direction. The volume elements simulation obtains the elongation of the gauge 0.04779 mm and the maximum stress in the notch region 692.878 Ν/mm^2. The shell elements simulation obtains the elongation of the gauge 0.04795 mm and the maximum stress in the notch region 675.179 Ν/mm^2. Results shows shell elements simulation will fulfill the analysis requirement.

        1.2 ANSA-Lsdyna-Meta workflow

        In this section, trimmed-IGA-based ANSA-Lsdyna-Meta workflow is used. Lower polynomial order and refined mesh can be a good strategy to avoid the cross-talk effect in the notched part. In this case, biquadratic polynomial order and 1mm span are used. The simulation obtains the elongation of the gauge 0.04785 mm and the maximum stress in the notch region 635.151 Ν/mm^2.

        1.3 ANSA- Kratos -Rhino workflow

        Additionally, IGA-based ANSA–Kratos–Rhino workflow is used. The IGA trimmed model and the boundary conditions are consistent with Section 1.2. The ANSA–Kratos plugin generates the input data (5 JSON files) required by the Kratos solver. Then Kratos performs the simulation and produces the result files which will be visualized by Rhino–Crocodile plugin. The simulation obtains the elongation of the gauge 0.0476 mm and the maximum stress in the notch region 655.35 Ν/mm^2.

        Figure 2: ANSA-Kratos-Rhino workflow

        Table 2: Results from experiment, theory and different workflows

        MethodMin mesh size (mm)Elongation (mm)Max stress (N/mm^2)Deviation from Experimental elongation (%)
        Theoretical/0.04276487.29
        Experiment/0.0398/0
        FEM Volume elements0.150.04779692.87820.08
        FEM Shell elements0.150.04784675.17920.2
        IGA (LSDYNA)10.04785635.15120.23
        IGA (Kratos)10.0476655.3519.6

        Table 2 shows details of simulation results. The findings demonstrate that IGA provides acceptable results, closely matching the experimental data.

        2 REFERENCES

        • [Tuo21] Tuo, Z., Yue, Z., Zhuang, X., Min, X., Badreddine, H., Qiu, L., Gao, J., 2021. Comparison of two uncoupled ductile damage initiation models applied to DP900 steel sheet under various loading paths. International Journal of Damage Mechanics 30, 25–45.
        • [Pe08] Peterson, R. E., & Pilkey, W. D., 2008. Peterson’s stress concentration factors (3rd ed., pp. 57–134). Wiley.

        DC9 Model Order Reduction of coupled vibro-acoustic systems:

        The vibro-acoustic field is a prime example where geometry-enhanced methods can have a strong influence on current applications. At mid- and high-frequency problems, where the acoustic wavelength is typically smaller or of the same order than the geometrical details within the object considered, the accuracy of the geometric representation is the driving factor for the analysis. To suppress the high computational burden, DC9 will work towards efficient model order reduction (MOR) schemes for one-way coupled and fully coupled vibro-acoustic analysis where the acoustic domain is described using an isogeometric boundary discretisation

        Model order reduction of isogeometric BEM systems

        In the first part of the PhD, DC9 worked on model order reduction strategies for the acoustic boundary element method within the IGA-framework (IGABEM) to speed-up the computational time and alleviate the memory costs. Specifically, an automatic MOR scheme based on Krylov subspaces recycling [1] is applied in conjunction with IGABEM to solve acoustic problems. The automatic MOR method automates the selection of Krylov subspaces to be recycled and creates a projection basis which sufficiently approximates the solution of the full-order model (FOM). The projection basis is used in combination with a Chebyshev polynomial approximation of the IGABEM system to create a reduced-order model (ROM), thus alleviating the computational burden.

        In the following the presented method is investigated and verified by an analytical solution, the sphere. A plane wave with an unit amplitude is propagating through an unbounded air volume and is impinging from the left on a rigid sphere (particle velocity = 0 m/s for the entire boundary surface). The acoustic pressure is studied in a frequency range from 1 to 1000 Hz with a step size of 1 Hz. The sphere consists of 6 identical conforming patches and is modeled with NURBS polynomial of degree four. The sphere consists of 726 degrees of freedom (DoFs).

        The following Figure shows that the sound pressure level (SPL) of the ROM matches that of the FOM, implying that the set threshold for the MOR method is sufficient. By inspecting the normalized residual for the defined frequency range, the residual is below the predefined tolerance. To generate the projection basis, the AKR algorithm only requires 3 full assemblies and 2 partial solutions, producing a basis spanning a subspace of 62 dimensions.

        Figure 1: (a) Sound pressure level at an evaluated point; (b) Normalized residual error of the double layer potential of ROM; (c) Configuration for AKR with x(ω), ω ∈ Ω

        References

        [1] . Panagiotopoulos, W. Desmet, and E. Deckers, “An automatic krylov subspaces recycling technique for the construction of a global solution basis of non-affine parametric linear systems,” Computer Methods in Applied Mechanics and Engineering, vol. 373, p. 113510, 2021.

        DC8 Efficient unbounded acoustic analysis starting from CAD:

        DC8 is focusing on implementing the fundamentals of IGA workflows for the acoustic simulations in time domain for unbounded (external) domains with the help of model order reduction (MOR) schemes.

        Update 1

        To achieve unbounded time domain analysis, time domain boundary element method (TD-BEM) is chosen such that Sommerfeld radiation condition is satisfied exactly and time simulations are made with applying convolution as matrix multiplication. Since there were no readily available or open source TD-BEM application, a 2D external/internal code was build on top of the open source framework OpenBEM1 in Matlab. An academic example of a 2D rectangular geometry with a sudden application of pressure example is examined and the results are given below and compared with analytical solutions.

        1. http://www.openbem.dk/ ↩︎

        DC7 Complex Constitutive modelling for immersed shell discretizations:

        Immersed Isogeometric Analysis for structural dynamics applications.

        Introduction and Motivation

        Isogeometric analysis (IGA) is able to provide numerous advantages compared to conventional Finite Element Analysis (FEA). IGA introduced originally as an idea of integrating Computer Aided Design (CAD) and FEA, providing efficiency by eliminating the meshing and re-meshing processes. Another strong motivation in its initial conception was the reduction or elimination of errors due to geometric approximation. The effect of higher-order inter-element continuity due to the smoothness of its spline based basis functions, leads to improved spectral accuracy over classical finite elements analysis. For this reasons IGA has been established as a strong high-order competitor in the field of structural dynamics [1].

        The maximum discrete eigenfrequency is a value dictating the critical time-step ensuring stability in hyperbolic problems with explicit time integration. In standard finite elements having only the possibility of C 0 continuity, using higher-order discretizations will result in high maximum discrete eigenfrequencies (both in the conforming and immersed setting). This will lead to tiny time steps and infeasible simulation times. Immersed IGA allows for the analysis of complex geometries without the need for conforming meshes, simplifying the meshing process. Employing immersed IGA and a Lumped Mass Matrix demonstrated promising results restricting the largest discrete eigenfrequency, as presented originally in [2]. Thus, with this configuration the time step size is practically independent of the presence of small cut elements, making immersed methods an attractive alternative for explicit dynamics simulations.

        Inspired by this observation, the studies conducted in [3] show that although this statement is valid, allowing the use of larger time steps, the lumping of the mass presents significant disadvantages in terms accuracy, especially in the immersed setting.

        These studies conducted on simple example problems, in particular for the case of a bar with free-free boundary conditions. The physical domain of a bar is immersed in an extended domain where the parameter ζ is dictating the size of the fictitious domain on each side, and for ζ = 0 the boundary fitted scenario applies.

        Figure 1: Maximum eigenvalue (left) and error in the time domain (right) for computed without stabilization (α min = 0), consistent mass matrices and ∆t = 0.9 ∆t critical . Dashed lines indicate ∆t.

        In Figures 1 and 2 we can observe the compromise between efficiency and accuracy for IGA as stated in [3]. On the one hand, consistent mass matrix provides better accuracy but smaller time steps and, on the other hand, lumped mass matrix, being able to bound the max eigenfrequency independently of the cut element size (parameter z), is providing larger time steps (dashed line), but with an accuracy being at least an order of magnitude less.

        Figure 2: Maximum eigenvalue (left) and error in the time domain (right) for computed without stabilization (α min = 0), lumped mass matrices (row summing for B-splines, HRZ lumping for Lagrange polynomials) and ∆t = 0.9 ∆t critical . Dashed lines indicate ∆t.

        In conclusion, mass lumping in combination with immersed IGA has some interesting advantages, although seems like this topic is not fully understood at the moment, with room for improving the established workflows. Additionally, some interesting research directions would be towards an improved mass lumping technique, as well as introducing adaptive refinement, resulting in a better alternative for the efficiency-accuracy dilemma.

        Research progress

        The initial scope of this work is to extend the studies performed in [3] to more complicated problem scenarios and soon including also a Kirchhoff-Love shell, where IGA has demonstrated its success. With that in mind, the first step was the extension to Bernoulli-Euler beams, since they are the first step towards shells, and it was the simplest way to see the possible behavior changes moving from 2 nd to 4 th order PDEs (see Fig. 3 and Fig. 4).

        In brief, we can comment that the extension of these results to the Bernoulli-Euler beam problem demonstrates a similar behavior compared to the bar with some minor differences.

        The next steps were the extension to two-dimensional problems, namely the equivalent of the bar, which is the membrane, and of the Bernoulli-Euler beam, which is the Kirchhoff plate.

        As we can see, the 2D results (Fig. 5) as expected, follow the same trend as the equivalent 1D results. Most of the obtained results are not herein presented for the sake of brevity but, as a general comment, we can state that they all mimic the equivalent 1D results. In addition, the transition from consistent to lumped mass matrices results in the bounding of the maximum eigenfrequency, in a similar way as in 1D.

        Figure 3: Largest eigenvalue ω max over fictitious domain size ζ in the Bernoulli-Euler beam problem, for computations without stabilization (left) and with stabilization (right) and a consistent mass matrix.
        Figure 4: Largest eigenvalue ω max over fictitious domain size ζ in the Bernoulli-Euler beam problem, for computations without stabilization (left) and with stabilization (right) and a lumped mass matrix.
        Figure 5: Largest eigenvalue ω max over fictitious domain size ζ in the membrane problem, for computations without stabilization (left) and with stabilization (right) and a consistent mass matrix.

        References

        [1] J. Cottrell, A. Reali, Y. Bazilevs, and T. Hughes, “Isogeometric analysis of structural vibrations,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 41, pp. 5257–5296, 2006.

        [2] L. Leidinger, “Explicit isogeometric b-rep analysis for nonlinear dynamic crash simulations: Integrating design and analysis by means of trimmed multi-patch shell structures,” Ph.D. dissertation, Technische Universitat Munchen, 2020.

        [3] L. Radtke, M. Torre, T. J. Hughes, A. Duster, G. Sangalli, and A. Reali, “An analysis of high order fem and iga for explicit dynamics: Mass lumping and immersed boundaries,” International Journal for Numerical Methods in Engineering, vol. 125, no. 16, p. e7499, 2024.

        DC6 Mathematical tools for immersed IGA:

        Development of mathematical tools for immersed IGA, related, in particular, to higher order phase-field modeling.

        Hello World!

        #include <iostream>
        
        int main() {
            std::cout << "Hello World!";
            return 0;
        }

        My name is Lucas and here I will write some highlights of my research within the Gecko Doctoral Network!

        During the past months, I have been studying higher order phase-field formulations, notably the Cahn-Hilliard equation, used to simulate phase separation phenomena. The smooth basis functions of Isogeometric Analysis (IGA) simplify the resolution of fourth-order spatial derivatives, enabling enhanced accuracy and computational efficiency compared to second-order models, even with coarser meshes. Local adaptive refinement schemes using truncated hierarchical B-splines (THB-splines) further optimize computational costs by increasing the resolution only in regions where it is needed; in phase-field problems, these regions are the interfaces.

        Recently, we developed an IGA-based phase-field implementation of the Cahn-Hilliard equation within G+Smo, an open-source C++ library that brings together mathematical tools for geometric design and numerical simulation. This implementation includes an adaptive scheme for solving the Cahn-Hilliard equation using truncated hierarchical B-spline (THB-spline) basis functions, which are particularly well-suited for adaptive refinement in the context of IGA, as they preserve the properties of hierarchical B-splines, such as linear independence and non-negativity, and also form a partition of unity.

        This implementation provides the foundation for introducing other higher-order phase-field models in the library and serves as a prototype for solving the Cahn-Hilliard equation with a locally adaptive scheme on any single-patch geometry parametrization that exists therein.

        Stay tuned for more news coming up soon!

        DC5 Large deformation structural elements (beams and shells) modeled with IBRA, including trimming and multiple coupled patches:

        Hello there!


        I’m Maram Alkhlaifat, one of the PhD candidates in the GECKO network, and I’m thrilled to share a glimpse into my research journey so far.

        My work revolves around structural elements for IBRA. What makes this work exciting? Well, imagine trying to design a component of Harry Potter’s flying car. Traditionally, you’d need to bounce back and forth between your design software and analysis tools, losing precision and time in the process. Our implementation aims to bridge this gap, allowing seamless integration between design and analysis.

        Started with Reissner-Mindlin shell element formulation in Kratos Multiphysics and tested to validate the coupled membrane-bending behavior, paving the way for even more complex challenges. Up next? Think nonlinear materials, large deformation theories, locking and yes—trimming and multipatch coupling to handle even the most intricate geometries. Here’s a sneak peek:

        P.S. Whether you’re interested in shell or beam elements, curious about CAD-CAE integration, or just want to geek out about computational mechanics, let’s connect—I’d love to chat!

        DC4 Co-simulation strategies involving IBRA for solution of multi-field problems:

        Complex technical systems often require a partitioned approach to enable disciplinary modelling and simulation with bestsuited solution approaches and discretization techniques in each domain.

        Developing Flexible and Accurate Strategies for Multiphysics Partitioned Simulations

        As part of the GECKO Project, my PhD research focuses on developing computational methods for partitioned multi-disciplinary simulations, with particular attention to Fluid-Structure Interaction (FSI) problems. The aim is to address the challenges of coupling different discretization techniques while ensuring robust, accurate, and scalable simulations that can be applied to a variety of engineering problems.

        One of the key objectives of this work is to enhance the efficiency and flexibility of partitioned simulation approaches. These methods allow different physical domains, such as fluids and structures, to be modeled using specialized solvers, which are then coupled and exchange data across their interfaces. This requires the development of advanced data transfer operators that ensure accurate mapping between domains, even when the underlying discretizations differ significantly.

        During the first year, significant progress was made in tackling these challenges. Among the main achievements is the implementation of robust coupling schemes that integrate Isogeometric B-Rep Analysis (IBRA) with other solvers, leveraging IBRA’s ability to provide highly accurate geometric representations. Additionally, efforts were focused on improving data transfer methods, particularly through the implementation of the Mortar Mapper, which is capable of handling both body-fitted and unfitted discretizations.

        Another important outcome of this research is a new approach to addressing singularity issues that arise when imposing strong boundary conditions on unfitted meshes. This development expands the applicability of partitioned simulations, enabling accurate and reliable modeling in scenarios involving unfitted discretizations. The methods developed during this first phase have been implemented within the Kratos-Multiphysics open-source framework, ensuring accessibility and further development opportunities.

        The results achieved so far provide a strong foundation for addressing more complex engineering problems in fields such as aerospace, biomechanics, and civil engineering. Future work will involve publishing these methods, comparing the performance of different data transfer operators, and applying the techniques to benchmark cases to validate their efficiency and accuracy further.

        The research represents a step forward in providing tools that enhance the reliability and flexibility of simulations for multi-physics problems, contributing to the design of safer and more efficient engineering solutions.

        DC3 Application of IBRA-type discretizations in implicit contact mechanics:

        Use of smooth CAD discretizations in contact mechanics is known to be beneficial.

        A brief comparison: IGA and FEM

          A key step in advancing the use of IGA was the implementation of a general body-fitted problem, where a flexible, nonlinear mapping between the parameter and physical spaces was developed. This mapping leverages IGA’s inherent capability to describe CAD geometries perfectly, ensuring precise simulation of even the most complex domains.

          Advantages of IGA over FEM in Body-Fitted Scenarios

          Comparative studies between IGA and FEM have revealed several advantages of IGA in body-fitted scenarios:

          1. Reduced Degrees of Freedom (DOFs):
            IGA requires fewer DOFs to achieve the same error level compared to FEM. This is due to the higher-order continuity of NURBS, which reduces the number of elements and control points needed for accurate approximations.
          2. Exact Geometry Representation:
            Unlike FEM, where mesh generation can introduce geometric inaccuracies, IGA ensures that the computational domain is an exact replica of the CAD model (Figure 1). This exactness is particularly beneficial for problems where geometric fidelity is critical, such as those involving contact mechanics.
          3. Simplified Refinement:
            Refinement in IGA can be achieved without altering the underlying geometry, making it more straightforward to adapt the computational model to different levels of precision.
          4. Higher Convergence Velocity:
            Both FEM and IGA exhibit similar convergence orders when using the same polynomial degree. However, IGA’s ability to easily increase the order of basis functions without remeshing provides a significant edge. Higher-order basis functions lead to enhanced convergence rates, enabling IGA to achieve desired accuracy more efficiently and lower numbers of degrees of freedom (Figure 2).

          Figure 1: Comparison of IGA and FEM discretizations.

          Figure 2: Comparison of IGA and FEM convergence and DOFs utilization.



          The Shifted Boundary Method

          The Shifted Boundary Method (SBM) offers a flexible framework to address challenges in numerical integration over complex domains, including contact mechanics. As outlined in, SBM shifts the imposition of boundary conditions from the true boundary Γ to a surrogate boundary h ​, composed of edges of a computational grid.

          Boundary conditions are modified using Taylor expansions, ensuring optimal convergence rates. This approach avoids challenges associated with small cut cells and simplifies numerical integration, making it particularly suited for embedded methods and large deformation problems.

          In this context, the SBM complements Isogeometric Analysis (IGA) by leveraging exact geometry descriptions and avoiding trimmed knot spans. This synergy improves the representation of physical geometries and enhances computational efficiency.

          Figure 3:  SBM main characteristics.

          The SBM for IGA has been implemented inside the Kratos Multiphysics framework (KratosMultiphycsGithub) to handle 2D fluid and structural mechanics problems with complex geometries.

          The SBM implementation preserves the optimal convergence of body-fitted cases under Dirichlet boundary conditions, though it experiences a one-order reduction in convergence for Neumann (load) conditions. This functionality has also been extended to 3D problems, broadening its applicability (Figure ).

          Figure 4: IGA+SBM for 2D/3D problems.

          Finally, Figure 5 shows a convergence comparison between standard IGA and IGA with SBM for the case of circular geometry. No big loss is evident for the use of the SBM against the exact geometry of the body-fitted scenario.

          Figure 5. Comparison of IGA body-fiitted and IGA + SBM convergence for a circle.

          Towards a “penalty-free” Contact mechanics

            Contact mechanics investigates the interaction of surfaces under load, with emphasis on stress distribution, deformation, and phenomena such as friction and adhesion. This field is essential in various engineering applications, including mechanical systems, aerospace structures, and biomechanics, where precision is paramount for safety and performance.

            1. Governing equations and variational formulation

            Contact mechanics deals with the study of stresses and deformations arising at the interface between contacting bodies. In computational mechanics, the weak form of the governing equations is often employed to facilitate numerical implementation.

            The formulation of contact mechanics begins with the total elastic potential energy of the system, expressed as:

            where Ψ is the strain energy density, C​ denotes the contact interface, λ is the Lagrange multiplier associated with the contact stress, and g​ is the gap function, which measures the normal separation between contact surfaces.

            By introducing the variation of the potential energy, the weak form of the contact problem becomes:

            where σ is the Cauchy stress tensor, and ϵ is the strain tensor.

            1. Approaches to Enforcing Contact Constraints

            Several methods have been developed to enforce contact constraints, each with its own advantages and limitations. Two of the most widely used are the Lagrange multiplier method and the penalty method, both of which are foundational to modern contact mechanics algorithms.

            Lagrange Multiplier Method

            The Lagrange multiplier method introduces an additional field, λ\lambdaλ, to explicitly enforce the contact constraints. The augmented variational form becomes:

            where λ acts as a contact force. This approach ensures exact constraint satisfaction but increases the number of unknowns in the system. The method is robust but computationally expensive due to the saddle-point nature of the resulting system. Furthermore, poor conditioning of the system matrix can complicate numerical solution strategies.

            Penalty Method

            The penalty method simplifies the implementation by replacing the constraint with a penalty term in the variational form:

            where ϵ > 0 is the penalty parameter. Larger values of ϵ enforce the constraints more strictly but can lead to numerical instability, while smaller values introduce constraint violations. Striking the right balance requires careful calibration, which can be problem-specific and nontrivial.

            Nitsche’s “penalty free” Method

            Nitsche’s method combines the strengths of both approaches by weakly enforcing contact constraints without introducing a Lagrange multiplier or requiring explicit penalty parameters. This formulation introduces stabilization terms that ensure numerical consistency and stability. The main idea is to substitute the lagrangian multiplier, , which physically represents the normal stress at the contact boundary, exactly with the stress at the contact, , that it’s directly computed from the displacements of the bodies involved in the contact and does not require additional degrees of freedom in the system. We can choose  

            to set the contact stress closer to the master or slave contact stresses. In this way the perturbation to the potential results as 

            1. Results

            To validate the implementation, two classical benchmark problems and a case comparison with the classical Contact Lagrangian FEM were analyzed:

            1. Patch Test

            The patch test evaluates the ability of the algorithm to maintain stress and displacement continuity across a contact interface between two squares with non-coincident meshes. The Nitsche penalty-free formulation demonstrated excellent accuracy, achieving stress and displacement continuity within numerical tolerances. This result confirms the robustness of the nearest projection search and the contact activation/deactivation strategy.

            Figure 6:  Patch test, Mesh (above), vertical-displacent (left ) vertical stress (right)

            2. Hertz Contact Problem

            The Hertz contact problem involves the interaction of a cylinder (or circle in 2D) with a rigid wall, serving as a well-known analytical reference for contact mechanics. The implementation accurately reproduces the contact pressure distribution, confirming the capability of the algorithm to capture normal stresses and displacements. In this case:

            • The pressure distribution follows the classical parabolic shape, consistent with the analytical solution.
            • The vertical displacement along the contact interface matches the derived theoretical solution, further validating the accuracy of the contact formulation.

            These benchmarks highlight the robustness and accuracy of the contact algorithm, demonstrating its applicability to a wide range of contact problems, from simple linear cases to more complex non-linear scenarios. The results also underscore the advantages of the penalty-free Nitsche formulation in maintaining stability and avoiding the pitfalls associated with traditional penalty-based methods.

            Figure 7: Hertz Circle-Wall contact, stress comparison with the true solution on the contact boundary. 

            3. Comparison with FEM

            Finally a punch test has been checked against a well tested FEM contact solver with ALM. The penalty-free/IGA body fitted algorithm proves to converge fast the FEM accurate solution. 

            Figure 8:Punch test horizontal displacements, FEM (above) vs. IGA(below)

            DC2 IBRA-type discretizations in computational solid mechanics:

            COMPARISON OF IGA & FEM IN COMPUTATIONAL SOLID MECHANICS

            In the field of computational solid mechanics, the prediction of mechanical behavior and damage evolution in materials is a critical challenge. Among the most widely used numerical techniques for solving solid mechanics problems, Isogeometric Analysis (IGA) and Finite Element Method (FEM) stand out due to their ability to model complex material behaviors, including damage, cracks, and deformations. This chapter provides an in-depth comparison of IGA and FEM in the context of computational solid mechanics, focusing on damage mechanics and the ability to simulate material failure and crack propagation.

            1 Introduction to Isogeometric Analysis and Finite Element Method

            Isogeometric Analysis (IGA), introduced by Hughes et al. (2005), integrates Computer-Aided Design (CAD) representations directly into the analysis process, enabling a seamless representation of geometry and the solution fields. The basis of IGA is the use of Non-Uniform Rational B-Splines (NURBS) or B-splines for both the geometry and the solution field (e.g., displacements, stress, etc.). This exact representation of geometry provides significant advantages over traditional finite element methods.

            The Finite Element Method (FEM), on the other hand, has been the standard numerical method for solving solid mechanics problems for decades. FEM discretizes the domain into small elements, where the problem is solved numerically using polynomial approximations. While FEM has proven to be highly versatile and effective in a wide range of problems, it typically involves approximations of the geometry, leading to potential errors, especially in complex structures.

            Both methods are widely used in computational solid mechanics, including in the study of damage mechanics, where the goal is to predict the evolution of damage in materials under various loading conditions. However, their differences in geometry representation, element formulation, and crack modeling make them suitable for different types of damage problems.

            2 Geometry Representation and Continuity

            • Isogeometric Analysis (IGA): A key feature of IGA is its exact representation of geometry. Using NURBS or B-splines, IGA can represent curved and complex geometries exactly without the approximation errors seen in FEM. This is particularly beneficial when modeling intricate geometries, such as those encountered in damage mechanics simulations involving cracks or material interfaces. IGA also provides higher-order continuity (C1, C2, etc.) in its shape functions, making it ideal for simulating smooth deformation and material behavior.
            • Finite Element Method (FEM): FEM approximates geometry using a discretized mesh, which can lead to errors when modeling complex geometries. The accuracy of the geometry representation in FEM depends on the mesh resolution and element type used. FEM typically employs lower-order continuity (C^0), meaning that while displacement continuity is maintained, higher-order derivatives such as stress and strain may exhibit discontinuities. This lower continuity can be a limitation when dealing with smooth material behavior but is suitable for handling discontinuous damage such as cracks.
            • Comparison: While IGA provides a more accurate and continuous representation of the geometry, FEM is more flexible in handling different types of element formulations, which makes it more adaptable to problems involving localized damage and cracks.

            3 Element Design and Mesh Sensitivity

            • Isogeometric Analysis (IGA): IGA uses higher-order elements that are capable of capturing complex material behavior with fewer degrees of freedom compared to FEM. The higher-order basis functions in IGA provide better accuracy for problems involving smooth deformations and damage evolution. In many cases, IGA requires fewer elements to achieve the same level of accuracy, especially for problems involving complex geometries.
            • Finite Element Method (FEM): FEM typically uses lower-order elements that require finer meshing and local refinements to achieve the desired level of accuracy. In problems involving damage localization, such as crack propagation, fine mesh refinement near the crack tip is often necessary. Mesh sensitivity is a common issue in FEM, especially when dealing with stress concentrators or sharp damage features.
            • Comparison: While IGA generally requires fewer elements for high accuracy in smooth problems, FEM’s flexible meshing and adaptive refinement capabilities make it better suited for problems involving sharp damage features, such as cracks, where localized refinement is needed.

            4 Modeling Cracks and Damage Propagation

            • Isogeometric Analysis (IGA) and Crack Propagation: A significant challenge in damage mechanics is simulating the propagation of cracks. IGA’s smooth basis functions are not naturally suited for capturing the discontinuities introduced by cracks, as they lead to sharp changes in displacement and stress fields. The inherent higher continuity of IGA can create difficulties in modeling crack initiation and propagation in a realistic manner, especially in mode I fracture where stress singularities are present.
              • Approaches to Handle Cracks in IGA: To address this challenge, several strategies have been developed:
                • Phase-Field Models: These models introduce a continuous representation of cracks, simulating crack growth as a gradual transition from intact to damaged material. Phase-field models are particularly well-suited to IGA, as they allow for the smooth evolution of damage without the need for sharp discontinuities (Miehe et al., 2010).
                • Extended Finite Element Method (XFEM): XFEM can be combined with IGA to enrich the basis functions near the crack region, allowing for the representation of cracks without the need for remeshing (Belytschko et al., 2009).
            • Finite Element Method (FEM) and Crack Propagation: FEM, due to its lower continuity nature, is better suited for simulating cracks and discontinuous damage. FEM can easily accommodate cracks using various methods:
              • Extended Finite Element Method (XFEM): XFEM extends the FEM by enriching the displacement field near crack regions to capture discontinuities in displacement and stress without the need for remeshing (Moës et al., 1999).
            • Comparison: While IGA faces challenges in directly simulating cracks due to its smoothness and higher continuity, techniques like phase-field models, and XFEM allow IGA to simulate cracks effectively. On the other hand, FEM is naturally suited to handle discontinuous damage and crack propagation due to its ability to introduce sharp discontinuities and perform adaptive meshing.

            5 Computational Efficiency

            • Isogeometric Analysis (IGA): IGA’s exact representation of geometry and higher-order elements make it computationally efficient for problems with smooth damage evolution and complex geometries. The use of fewer elements can lead to significant computational savings. However, when applied to highly nonlinear damage models, the integration over higher-order elements can be computationally expensive (Hughes et al., 2005).
            • Finite Element Method (FEM): FEM is highly adaptable and can be computationally efficient, especially when adaptive meshing or parallel processing techniques are used. However, problems involving cracks or localized damage often require fine mesh refinements, increasing computational cost. Remeshing during crack propagation can add further overhead (Zienkiewicz et al., 2005).

            Comparison: IGA is more efficient for problems with smooth damage evolution and less complex crack behavior, while FEM is better suited for large-scale problems involving crack propagation and local damage.

            ISOTROPIC DAMAGE MODELS

            1 Isotropic Damage Models in One-Dimensional Systems

            A critical area of study in damage mechanics is the development of isotropic damage models, which describe the degradation of material properties due to loading and unloading cycles. A commonly used formulation for such models is that introduced by Oliver et al. (Oli+90), which characterizes the stress-strain relationship of materials undergoing damage.

            In the 1D isotropic damage model, the material response is governed by the evolution of a damage variable, which reduces the effective stiffness of the material as damage accumulates. This model accounts for both loading and unloading behavior, and its implementation involves solving the system of equations that govern the damage evolution process. For instance, by considering a material with properties such as Young’s modulus (E) of 1E+08 Pa, a yield stress of 2 MPa, and a fracture energy of 5E+04 J/m², the material’s degradation can be studied by solving the stress-strain relations for different loading paths, as illustrated by the stress-strain curve (Figure 3.1).

            Figure 1. Stress-strain curve during both loading and unloading phases for the given material properties, with a Young’s modulus (E) of 1.00E+08 Pa, yield stress of 2.00E+06 Pa, and fracture energy of 5.00E+04 J/m2

            2 Extension to Two-Dimensional Damage Models

            Building on the 1D isotropic damage model, a natural progression is to extend the formulation to two-dimensional systems. In these models, the material exhibits strain-softening behavior, which requires careful attention to regularization techniques to prevent mesh dependence in simulations. The characteristic length (lc) is a critical parameter in this context, as it determines the scale at which damage is regularized, ensuring that damage and fracture processes are not overly influenced by the numerical mesh.

            In traditional FEM formulations, the characteristic length is computed based on the maximum distance from the centroid of an element to its nodes. For a triangular element, the characteristic length lc is given by:

            lc=max⁡(∥𝐱c−𝐱i∥),i=1,2,3

            where xc\mathbf{x}_c is the centroid and xi\mathbf{x}_i are the coordinates of the element’s vertices. This approach is effective in FEM, where the mesh is discretized with nodes. However, this definition is not directly applicable to Isogeometric Analysis (IGA), where elements are defined in terms of NURBS (Non-Uniform Rational B-Splines) basis functions rather than nodal points.

            3 Characteristic Length in Isogeometric Analysis (IGA)

            In IGA, the elements are not tied to a discrete mesh of nodes, but rather defined by the continuous nature of the NURBS basis functions. To maintain the role of the characteristic length in regularizing damage mechanics models within IGA, a new formulation is needed.

            The characteristic length lc in IGA is defined based on the parametric spans of the NURBS basis functions, rather than relying on the geometry’s nodal discretization. For instance, for a quadrilateral element, the characteristic length can be expressed as:

            lc=⁡(॥𝐱(ui+1,vj)−𝐱(ui,vj)॥,॥𝐱(ui,vj+1)−𝐱(ui,vj)॥),

            where 𝐱(u,v) is the mapping from the parametric space to the physical space, and uu and vv are the parameters defining the NURBS surface. This definition ensures that the characteristic length in IGA reflects the smooth, continuous nature of the geometry and maintains the regularization properties needed for accurate damage modeling.

            4 Role of Regularization in Damage Mechanics

            In damage mechanics, the regularization of strain-softening behavior is critical to obtaining stable and physically meaningful solutions. Without regularization, damage models can exhibit mesh dependence, where the results vary significantly with the mesh size, leading to unphysical crack patterns and non-converging solutions. The characteristic length lc, whether defined in FEM or IGA, serves as a scale factor that regularizes the damage evolution, ensuring that the results are independent of the discretization used.

            The regularization process is particularly important in the simulation of fracture and material degradation, where local damage must be captured without introducing numerical artifacts. By defining lc in a way that is consistent with the underlying geometry, as in the case of the IGA-based definition, the approach mitigates mesh dependency and allows for more accurate simulations of complex material behavior.