Isogeometric Element Research Articles

Method of virtual sources using on-surface radiation conditions for the Helmholtz equation

We develop a novel method of virtual sources to formulate boundary integral equations for exterior wave propagation problems. However, by contrast to classical boundary integral formulations, we displace the singularity of the Green’s function by a small distance h>0. As a result, the discretization can be performed on the actual physical boundary with continuous kernels so that any naive quadrature scheme can be used to approximate integral operators. Using on-surface radiation conditions, we combine single- and double-layer potential representations of the solution to arrive at a well-conditioned system upon discretization. The virtual displacement parameter h controls the conditioning of the discrete system. We provide mathematical guidance to choose h, in terms of the wavelength and mesh refinements, in order to strike a balance between accuracy and stability. Proof-of-concept implementations are presented, including piecewise linear and isogeometric element formulations in two- and three-dimensional settings. We observe exceptionally well-behaved spectra, and solve the corresponding systems using matrix-free GMRES iterations. The method is implemented for problems with known analytical solutions for comparison. We conclude that the proposed method leads to accurate solutions and good stability for a wide range of wavelengths and mesh refinements.

On the use of mixed basis function degrees within a convective isogeometric element formulation

AbstractSolving linear elasticity problems using standard finite element methods, different locking phenomena can occur. In order to counteract these effects, mixed methods or formulations using higher approximation orders can be employed, for instance. As a consequence, this leads to an increased computational effort. Hence, a selective elevation of orders in decisive directions within a purely displacement‐based element formulation is proposed in this contribution. Within isogeometric analysis (IGA), the geometry is discretized using non‐uniform rational B‐splines (NURBS), which simultaneously represent the basis functions for the analysis. Due to the fact that the geometry can be preserved exactly during analysis, this can increase the accuracy of results. In this contribution, convective basis systems that are aligned with the local geometry are employed combined with selective order elevation. The required convective basis systems are interpolated from those determined in each control point.

Three-dimensional bending and free vibration analyses of laminated cylindrical panel with/without elastic foundation using two-dimensional discrete method

In this paper, a novel isogeometric scaled boundary finite element method (IGSBFEM) is developed to provide a comprehensive investigation into bending and free vibration of laminated cylindrical panels resting on Pasternak foundation. The significant novelty of this method lies in that the proposed formulation is derived based on the three-dimensional (3D) theory of elasticity, while only two-dimensional (2D) reference surface associated with the curvilinear coordinate system needs to be discretized by using isogeometric elements constructed with the non-uniform rational B-splines (NURBS) basis functions, which allows to reduce the computation cost and maintain the exact geometry of problem with desired precision flexibly. The SBFEM governing equation of each single lamina is reduced to a matrix form of first-order ordinary differential equations by introducing a vector of dual variables made of nodal displacements and forces, and the solution can be expressed analytically in terms of the radial coordinate. Therefore, shear locking phenomenon can be overcome naturally without any shear correction factors or reduced integration. The global solving equations of laminated cylindrical panels resting on Pasternak foundation are assembled with a layer-wise method. In addition, the high-order continuous features of isogeometric elements are insensitive to distorted mesh. Therefore, it is compatible with the research of complex geometrical structures with holes. The superior performances of the current approach are testified in the analyses of cylindrical panels with or without elastic foundation, solutions obtained by many other researchers are used for comparison purpose, and effects of several parameters including length and thickness, fiber orientation, stacking sequence and foundation stiffnesses coefficient on the bending or free vibration characteristics are also discussed.

Three-dimensional vibration analysis of rotating pre-twisted variable thickness blades composed of spanwise graded functional materials

In this study, the isogeometric method is used to obtain a three-dimensional (3-D) solution for the vibration of rotating pre-twisted blades composed of spanwise functionally graded materials (SFGM). The material properties are assumed to change continuously in the spanwise direction according to a mixture rule. Nonuniform rational B–splines (NURBS) are used to describe the unknown displacement and geometry of the rotating blades, whose thickness linearly tapers along the spanwise direction. A typical isogeometric element for a rotating pre-twisted variable thickness SFGM blade is introduced to prevent transformations among multiple coordinate systems. This model considers centrifugal stiffening, centrifugal softening, and Coriolis force effects. Several cases are studied, including rotating rectangular plates, partial cylindrical shells and conical shells with different material properties and geometric parameters. The predicted vibration behaviors are verified by comparisons with results from the literature and commercial finite element software. The comparisons show that the current model has good accuracy. Finally, the influence of material properties, rotating speeds, installation angles, pre-twist angles, and variable thicknesses on the vibrational characteristics of rotating pre-twisted SFGM blades are examined.

Vibration and buckling analysis of agglomerated CNT composite plates via isogeometric analysis using non-polynomial shear deformation theory

The present study employs the non-polynomial shear deformation theory within the framework of isogeometric analysis to study the vibration and buckling behavior of carbon nanotubes reinforced composite plates. Owing to the higher aspect ratio of carbon nanotubes, they are prone to agglomeration in matrix. Therefore, to analyze the effect of carbon nanotubes agglomeration in composite plates, a two-parameter micromechanics model is used to determine the effective mechanical properties of the composite plate. Current study opts isogeometric analysis as numerical method, as it can represents the complex shape geometries, and provides higher-order element continuity. Higher-order element continuity is desirable for implementation of non-polynomial shear deformation theory. By the virtue of non-uniform rational B-splines, isogeometric analysis turns out to be better analysis numerical method for complex and shell structures. Governing isogeometric element equations are derived using Hamilton's principle. A detailed parametric study is also performed to investigate the effect of degree of agglomeration, thickness to length/radius ratio, different boundary conditions, and complex shape cutouts in a composite plate. The presented formulation is validated across analytical and finite element method solutions with available literature.

A boundary-penalized isogeometric analysis for second-order hyperbolic equations

Explicit time-marching schemes are popular for solving time-dependent partial differential equations; one of the biggest challenges these methods suffer is increasing the critical time-marching step size that guarantees numerical stability. In general, there are two ways to increase the critical step size. One is to reduce the stiffness of the spatially discretized system, while the other is to design time-marching schemes with larger stability regions. In this paper, we focus on the recently proposed explicit generalized-α method for second-order hyperbolic equations and increase the critical step size by reducing the stiffness of the isogeometric-discretized system. In particular, we apply boundary penalization to lessen the system’s stiffness. For pth order Cp−1 isogeometric elements, we show numerically that the critical step size increases by a factor of p2−3p+64, which indicates the advantages of using the proposed method, especially for high-order elements. Various examples in one and two dimensions validate the performance of the proposed technique.

A novel hybrid isogeometric element based on two-field Hellinger–Reissner principle to alleviate different types of locking

In the present work, novel hybrid elements are proposed to alleviate the locking anomaly in non-uniform rational B-spline-based isogeometric analysis (IGA) using a two-field Hellinger–Reissner variational principle. The proposed hybrid elements are derived by adopting the independent interpolation schemes for displacement and stress fields. The key highlight of the present study is the choice and evaluation of higher-order terms for the stress interpolation function to provide a locking-free solution. Furthermore, the present study demonstrates the efficacy of the proposed elements with the treatment of several two-dimensional linear-elastic benchmark problems alongside the conventional single-field IGA, Lagrangian-based finite element analysis (FEA), and hybrid FEA formulation. It is shown that the proposed class of hybrid elements performs effectively for analyzing the nearly incompressible problem domains that are severely affected by volumetric locking along with the thin plate and shell problems where the shear locking is dominant. A better coarse mesh accuracy of the proposed method in comparison with the conventional formulation is demonstrated through various numerical examples. Moreover, the formulation is not restricted to the locking-dominated problem domains but can also be implemented to solve the problems of general form without any special treatment. Thus, the proposed method is robust, most efficient, and highly effective against both shear and volumetric locking.

Discontinuous Galerkin isogeometric analysis with peridynamic model for crack simulation of shell structure

A discontinuous Galerkin (DG) isogeometric analysis (IGA) formulation with peridynamic (PD) model is proposed for the simulation of cracks in shell structure. A micro-beam PD bond which considers the planar, bending and shearing deformations of shell is applied. Then, based on the micro-beam bond, the peridynamic integral equation is implemented via DG formulation of a conventional isogeometric element. The proposed DGIGA implementation of PD model (DGIGA PD) assures the analysis of crack propagation with spline based smooth geometric model, which is important in the simulation of shell structure. The DGIGA PD calculation improves the computational precision through multiple points based Gaussian integration and provides derivative fields such as strain and stress. To further improve the computational efficiency, the local/nonlocal coupling algorithm for the DGIGA PD model and IGA model is proposed. The computational domain is described fully by spline mesh of IGA, within which the PD domain is assigned where the cracks are propagating. The force and moment balance principal is applied, and an element based stiffness replacing technique is proposed for the coupling process. The performance of the proposed formulation is verified by benchmark examples.

Unconditional stability in large deformation dynamic analysis of elastic structures with arbitrary nonlinear strain measure and multi-body coupling

Stability is a fundamental feature of a time integration method in long simulations of elastic bodies in large deformations. It is well known that energy conservation is the key to achieve it. The implicit one-step conserving method of Simo and Tarnow is simple and effective for quadratic strain measures. However, the issue of unconditional stability does not yet have a definitive solution for general nonlinear strain measures and multi-body couplings. This article shows how to reach this goal by a simple modification of the Simo and Tarnow method. In practice, the mean internal force vector of the time step is evaluated using the average value of the stress at the end-points and an integral mean of the strain–displacement tangent operator over the step computed by time integration points. Compared to other proposals, the approach conserves energy for any structural model regardless of its spatial finite element discretization and does not require Lagrange multipliers, discrete derivatives, approximations of the strain–displacement law or special iterative solutions. Long term stability is proved using dynamic analyses of Reissner beams and Kirchhoff–Love shells, discretized spatially with Lagrangian and isogeometric elements respectively. Moreover, the proposed energy-conserving scheme is extended to multi-body analyses with nonlinear couplings by means of a penalty formulation, with focus on the rotational continuity for generic multi-patch Kirchhoff–Love shells.

Isogeometric-analysis-based stiffness spreading method for truss layout optimization

A novel isogeometric analysis-based stiffness spreading method (IGA-based SSM) is proposed for the truss layout optimization design within a continuum, where the topology, geometry, and cross-sectional areas of the truss elements can be simultaneously optimized. The IGA-based SSM based on energy conservation is employed to obtain the spreading stiffness matrices of the truss elements that are embedded in weak IGA background grids. The IGA background grids are constructed using high-order continuous isogeometric elements that can overcome the limitation of discontinuous sensitivity derived by traditional C0 elements. Because the sensitivity required for the optimization can be obtained analytically, gradient-based optimization algorithms can be easily implemented. To verify the effectiveness of the IGA-based SSM for truss layout optimization design, typical illustrative examples are used, considering the compliance minimization optimization problem. The proposed method can provide a continuous and smooth sensitivity field and overcome the numerical inaccuracy caused by interpolation methods, such as the radial basis function. The different background mesh sizes, initial layouts, and orders p of nonuniform rational B-splines basis functions are investigated in truss layout optimization. The results indicate that the selection of the appropriate background mesh size can avoid the problem of discontinuous truss connections, C1-continuous isogeometric elements can provide better optimization designs with a lower computational cost, and the proposed method can provide a clear layout for different initial layouts. Moreover, the spatial truss layout optimization problem and the min–max displacement optimization problem are also studied to further verify the effectiveness of the proposed method. Finally, the proposed method is extended to solve the concentrated force diffusion problem that widely exists in the interstage adapters of launch vehicles by adding the variance constraint of reaction forces. The results indicate that the proposed method can obtain required concentrated force diffusion results with reasonable stiffness distributions.

A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes

We introduce a boundary penalization technique to improve the spectral approximation of isogeometric analysis (IGA). The method removes the outliers appearing in the high-frequency region of the approximate spectrum when using the Cp−1,p-th (p≥3) order isogeometric elements. We focus on the classical Laplacian (Dirichlet) eigenvalue problem in 1D to illustrate the idea and then use the tensor-product structure to generate the stiffness and mass matrices for multidimensional problems. To remove the outliers, we penalize the higher-order derivatives from both the solution and test spaces at the domain boundary. Intuitively, we construct a better approximation by weakly imposing features of the exact solution. Effectively, we add terms to the variational formulation at the boundaries with minimal extra-computational cost. We then generalize the idea to remove the outliers for the isogeometric analysis of the Neumann eigenvalue problem (for p≥2). The boundary penalization does not change the test and solution spaces. In the limiting case, when the penalty goes to infinity, we perform the dispersion analysis of C2 cubic elements for the Dirichlet eigenvalue problem and C1 quadratic elements for the Neumann eigenvalue problem. We obtain the analytical eigenpairs for the resulting matrix eigenvalue problems. Numerical experiments show optimal convergence rates for the eigenvalues and eigenfunctions of the discrete operator.

Efficient computation of nonlinear isogeometric elements using the adjoint method and algorithmic differentiation

We present a consistent and efficient approach to the formulation of geometric nonlinear finite elements for isogeometric analysis (IGA) and isogeometric B-Rep analysis (IBRA) based on the adjoint method. IGA elements are computationally expensive, especially for high polynomial degrees. Using the method presented here enables us to reduce this disadvantage and develop a methodical framework for the efficient implementation of IGA elements. The elements are consistently derived from energy functionals. The load vector and stiffness matrix are obtained from the first and second order derivatives of the energy. Starting from the functional, we apply the concept of algorithmic or automatic differentiation to compute the precise derivatives. Here, we compare the direct (forward) and adjoint (reversed) methods. Analysis of the computational graph allows us to optimize the computation and identify recurring modules. It turns out that using the adjoint method leads to a core-congruential formulation, which enables a clean separation between the mechanical behavior and the geometric description. This is particularly useful in CAD-integrated analysis, where mechanical properties are applied to different geometry types. The adjoint method produces the same results but requires significantly fewer operations and fewer intermediate results. Moreover, the number of intermediate results is no longer dependent on the polynomial degree of the NURBS. This is important for implementation efficiency and computation speed.The procedure can be applied to arbitrary element formulations and coupling conditions based on energy functionals. For demonstration purposes, we present the proposed approach specifically for use with geometrically nonlinear trusses, beams, membranes, shells, and coupling conditions based on the penalty method.

Split generalized-[formula omitted] method: A linear-cost solver for multi-dimensional second-order hyperbolic systems

We propose a variational splitting technique for the generalized-α method to solve hyperbolic partial differential equations. We use tensor-product meshes to develop the splitting method, which has a computational cost that grows linearly with the total number of degrees of freedom for multi-dimensional problems. We use the generalized-α method for the temporal discretization while standard C0 finite elements as well as isogeometric elements for spatial discretization. . We perform spectral analysis on the amplification matrix to establish the unconditional stability of the method and to show how splitting affects the overall behavior of the time marching scheme for finite time step sizes, including the standard stability analysis limits 0 and ∞ as particular cases. We use various examples to demonstrate the performance of the method and the optimal approximation accuracy. In these examples, we compute the L2 and H1 norms of the error to show the optimal convergence of the discrete method in space and second-order accuracy in time. Lastly, we also use these tests to demonstrate the linear cost of the solver as the number of degrees of freedom grows.

A New Lighting on Analytical Discrete Sensitivities in the Context of IsoGeometric Shape Optimization

Isogeometric shape optimization has been now studied for over a decade. This contribution aims at compiling the key ingredients within this promising framework, with a particular attention to sensitivity analysis. Based on all the researches related to isogeometric shape optimization, we present a global overview of the process which has emerged. The principal feature is the use of two refinement levels of the same geometry: a coarse level where the shape updates are imposed and a fine level where the analysis is performed. We explain how these two models interact during the optimization, and especially during the sensitivity analysis. We present new theoretical developments, algorithms, and quantitative results regarding the analytical calculation of discrete adjoint-based sensitivities. In order to highlight the versatility of this sensitivity analysis method, we perform eight benchmark optimization examples with different types of objective functions (compliance, displacement field, stress field, and natural frequencies), different types of isogeometric element (2D and 3D standard solids, and a Kirchhoff-Love shell), and different types of structural analysis (static and vibration). The numerical performances of the analytical sensitivities are compared with approximate sensitivities. The results in terms of accuracy and numerical cost make us believe that the presented method is a viable strategy to build a robust framework for shape optimization.

A mixed isogeometric analysis approach for the transient swelling of hydrogel

In this paper, a mixed isogeometric analysis approach is developed to realize accurate and efficient simulations of the transient swelling of hydrogel strongly coupled with large deformation and fluid diffusion. In this method, the mixed isogeometric element is presented for the coupled multi-physics problems accompanied with the volumetric locking caused by the incompressibility of polymeric chains and solvent molecules. Such type of element is based on the subdivision properties of the Non-Uniform Rational B-Spline (NURBS) basis functions, in which the fine and coarse basis functions with the same order are used to discretize the displacement and chemical potential fields, respectively. Meanwhile, an auxiliary configuration for the multiplicative decomposition of deformation gradient is introduced to eliminate the difficulty due to certain specific cases, such as a negative infinite chemical potential for the hydrogel in a dry state. Then, the virtual work equations are derived based on the strain energy function with the ability of describing the constitutive relation of the hydrogel for large deformation from any swollen state and the mass conservation. The discrete formulations of the linearized virtual work equations for two fields on the base of the mixed isogeometric element are obtained. The efficiency, accuracy and validity of the proposed mixed isogeometric analysis approach are illustrated by solving several representative problems. The applicability of the proposed method is further demonstrated through two practical applications, which indicates that it has great potential to serve as an effective numerical tool for mechanical analysis and design of many hydrogel-based bionic devices.

Isogeometric analysis of the new integral formula for elastic energy change of heterogeneous materials

In recent years, more and more attention has been paid to isogeometric methods, in which shape functions are used to accurately describe CAD models and approximate unknown fields. The isogeometric boundary element method (IGABEM) realizes the integration on the exact boundary of the region, which means there is no geometric discretization error. In this paper, a more general interface integral formula for elastic energy increment of heterogeneous materials is extended from previous work (Dong, 2018), in which the only unknown variable is the displacement located on the interface between matrix and inclusion. This feature makes it more compatible with boundary element method (BEM) because of none of volume parametrization. However, the geometry discontinuity on the boundary called ‘corner point problem’ increases the difficulty and decreases the accuracy when solving complex numerical examples when heterogeneous structures are considered. The discontinuous element method combined with IGABEM is presented and applied to deal with ‘corner point problem’. In the numerical examples, the interface between matrix and inclusion is discretized by quadratic isogeometric elements. Compared with the analytical solution, the numerical results show that the method has higher accuracy and efficiency.

New 25-point stencils with optimal accuracy for 2-D heat transfer problems. Comparison with the quadratic isogeometric elements

A new approach for the increase in the order of accuracy of high order elements used for the time dependent heat equation and for the time independent Poisson equation has been suggested on uniform square and rectangular meshes. It is based on the optimization of the coefficients of the corresponding discrete stencil equation with respect to the local truncation error. By a simple modification of the coefficients of 25-point stencils, the new approach exceeds the accuracy of the quadratic isogeometric elements by four orders for the heat equation and by twelve orders for the Poisson equation. Despite the significant increase in accuracy, the computational costs of the new technique are the same as those for the conventional quadratic isogeometric elements on a given mesh. The numerical examples are in a good agreement with the theoretical results for the new approach and also show that the new approach is much more accurate than the conventional isogeometric elements at the same number of degrees of freedom. Hybrid methods that combine the new stencils with the conventional isogeometric and finite elements and can be applied to irregular domains are also presented.

On the condition number of high order finite element methods: Influence of p-refinement and mesh distortion

In this article, the condition number of the stiffness matrix κ(K) is compared for three high order finite element methods (FEMs), i.e., the p-version of the FEM, the spectral element method (SEM), and the NURBS-based isogeometric analysis (IGA). Note that only problems in linear elasticity are considered in the analysis. It is well-known that the condition number is one factor strongly influencing the number of significant digits for direct solvers or the required iteration count for iterative solution schemes. Therefore, it is important to investigate the effect of the choice of the shape functions and the element distortion on κ(K). Based on numerous one- and two-, and three-dimensional examples, these influences are comprehensively studied, and the numerical results are compared with condition number estimates extracted from the literature. Overall, a good agreement is observed for p-FEM, SEM and two-dimensional IGA, while discrepancies are noted for three-dimensional isogeometric elements. These are important findings as theoretical results may be only available for very restricted scenarios, where one-element geometries, constant Jacobi matrices of the element maps, etc. are considered.

Accurate numerical solutions of 2-D elastodynamics problems using compact high-order stencils

A new approach for the increase in the order of accuracy of high-order numerical techniques used for 2-D time-dependent elasticity (structural dynamics and wave propagation) has been suggested on uniform square and rectangular meshes. It is based on the optimization of the coefficients of the corresponding semi-discrete stencil equations with respect to the local truncation error. It is shown that the second order of accuracy of the linear finite elements with 9-point stencils is optimal and cannot be improved. However, we have calculated new 9-point stencils with the second order of accuracy that yield more accurate results than those obtained by the linear finite elements. We have also developed new 25-point stencils (similar to those for the quadratic finite and isogeometric elements) with the optimal sixth order of accuracy for the non-diagonal mass matrix and with the optimal fourth order of accuracy for the diagonal mass matrix. The numerical results are in good agreement with the theoretical findings as well as they show a big increase in accuracy of the new stencils compared with those for the high-order finite elements. At the same number of degrees of freedom, the new approach yields significantly more accurate results than those obtained by the high-order (up to the eighth order) finite elements with the non-diagonal mass matrix. The numerical experiments also show that the new approach with 25-point stencils yields very accurate results for nearly incompressible materials with Poisson’s ratio 0.4999.

Weak coupling of nonlinear isogeometric spatial Bernoulli beams

This paper proposes a coupling scheme for isogeometric elements, which is valid for large rotations and displacements. Two elements can be coupled not only at interpolated control points but also inside the parametric domain of a NURBS patch. Furthermore, each cross section may be arbitrarily oriented in space. The formulation is applicable for all structural elements, where an orthogonal system of base vectors can be derived. An Euler–Bernoulli beam is used as an example in order to illustrate the demands on the coupling conditions in the proposed approach. Moreover, the conditions are chosen such that different types of joints, e.g. scissor joints, can be modeled and such that they are easy and fast to implement. The coupling is incorporated by an additional term in the Principle of Virtual Work which is here computed using a penalty approach. The Euler–Bernoulli beam is summarized in order to introduce the notations in this contribution. This is followed by the proposed coupling conditions in the weak form. Subsequently, alternatives for the coupling conditions are briefly introduced. Eventually, benchmark and demonstrator examples validate and illustrate the proposed coupling methodology.