1 Introduction

In the design process for commercial and industrial products, cost-effectiveness and the accuracy of the numerical simulation will be crucial. Although finite element analysis (FEA) is an effective numerical method for structural analysis, further improved efficiency and accuracy are desired. Finite element (FE) discretization based on a simple polygon often requires numerous meshes to express complicated curved objects. In particular, the surface expressed by the polygons may cause considerable discrepancy between computer-aided design (CAD) and FE discretization. Therefore, the representation of the geometry by FE mesh will be cumbersome. In addition, a stress evaluation procedure based on the conventional FEA will induce an inaccuracy, in which a Lagrange polynomial shape function of the \(C^{0}\) inter-element continuity will be used. Therefore, a computation of the nodal stress from the displacement requires an additional recovery technique such as an extrapolation above such geometrical and numerical discrepancies, and it will be greater for further complicated curved geometries such as an industrial component used in various engineering fields.

Isogeometric analysis (IGA) is one of the methods that alleviate such limitations. IGA has been proposed as a means to combine traditional CAD with FEA [1]. This method uses the same basis function as that employed by CAD, non-uniform rational B-spline (NURBS), to describe both the geometry and physical quantities in the field. The detailed CAD geometry reduces the approximation discrepancy induced from the polygon-based FE discretization. Because the NURBS basis function has a higher inter-element continuity, a more accurate numerical solution will then be obtained. IGA has been successfully applied in many engineering fields, such as structural vibration [2, 3], fluids [4], electromagnetics [5], and fluid–structure interaction analysis [6]. IGA results were more accurate and exhibited good convergence with fewer discretization refinements, compared to traditional FEA results. In the structural analysis, a precise and continuous stress field will be obtained because of its higher inter-element continuity from the NURBS basis function.

However, the application of NURBS-based IGA to arbitrary engineering geometries will not be straightforward because of the pre-processing difficulties. Because the NURBS basis expresses objects as a tensor product, many patches will be required for the complicated geometry. Moreover, an additional procedure will be accompanied to obtain higher-order continuity at the intersection among the patches. An analysis on a volumetric object encounters other difficulties. In most CAD software, a three-dimensional solid object is represented as a boundary-representation (B-rep). This is inappropriate for the numerical analysis of a volumetric solid owing to the absence of the inner information for the discretization. Numerous attempts have been made to overcome such limitation. Aigner et al. [7] proposed a swept volume parameterization for IGA based on the given boundary condition and guiding curves. Hsu et al. [8] developed a parametric modeling platform by implementing plug-ins into CAD software, which generated a CAD suitable for IGA. Bazilevs et al. [9] attempted IGA on the fluid-structural vibration of a gas-turbine blade of a helicopter using the NURBS mesh-generation method [10]. In addition, the solid IGA via the other splines, e.g., T-splines [11] and U-splines [12], was proposed.

As an alternative, FE discretization was considered. Although the traditional mesh generation was required, it generated the volumetric parameterization effectively regardless of the geometry. Engvall and Evans [13, 14] proposed surface-to-volume parameterization using various Bernstein–Bézier FEs. Ainsworth et al. [15] presented algorithms for an optimal assemblage of matrices for FEA framework with various Bernstein–Bézier simplicial elements. Although it was not an exact IGA, Kadapa [16] implemented the quadratic Bézier mesh into an FEA framework and proposed improved accuracy over traditional FEA. Those methods used lower-order FE and were relatively efficient to implement. However, the analysis was conducted using the Bézier basis function, which contains \(C^{0}\) inter-element continuity. Jaxon and Qian [17] and Xia and Qian [18] used the triangular and tetrahedral Bézier elements for outer NURBS surface mapping and an inner volumetric parameterization, respectively. They used a rational Bézier element to represent the exact geometry and construct a \(C^{1}\)-continuous basis function to approximate the geometry and physical quantities. However, such a method required the fifth-order element, which was expensive and rarely supported by the traditional FE pre-processor.

Inspired by IGA framework based on FE discretization, this paper will attempt to enhance structural FEA framework utilizing the advantageous features of IGA. The quadratic tetrahedron from the conventional FEA software will be converted into Bézier tetrahedron, which will be applicable to the complicated geometry by the unified pre-processing. The resulting surface reconstruction will significantly decrease the geometrical discrepancy of the traditional FE polygon. Although such concept was proposed by Kadapa [16] and Ainsworth [15], the present method includes the enhancement of the inter-element continuity. Additionally, the present approach is distinguished from Jaxon [17] and Xia [18] by using general tetrahedral FE.

An approximate \(C^{1}\) Bézier basis function will be constructed through the linear combination of \(C^{0}\) Bézier basis function and continuity coefficients. The continuity coefficients will be obtained by Hermite interpolation to intersections among tetrahedrons. The macro element split technique will be employed, and a minimal determining set (MDS) will be defined to represent the entire domain by continuity condition. The accuracy of the present method will be examined by several curved geometries for which an exact solution exists. Moreover, a large-scale analysis on the curve-dominant industrial geometry will be performed. This paper is organized as follows. Section 2 introduces Bézier mesh generation for enhancing geometric representation with a brief explanation of the Bernstein–Bézier FE. Section 3 describes constructing an approximate \(C^1\) basis function with general tetrahedral FE. In Sect. 4, the outline of the present approach is described. In Sect. 5, the present approach is validated by applying it to various curved geometries. Conclusions and future works are addressed in Sect. 5.

2 Numerical discretization

In this section, a pre-processing procedure for the numerical discretization will be introduced. Direct use of the general solid CAD for the numerical analysis will be a challenge due to the absence of the inner volumetric expression. Therefore, the present technique will utilize FE discretization from an existing mesh generator to obtain the inner volumetric information. To reduce the geometric discrepancy between FE discretization and CAD at the exterior curved surface, polygonal FE on the outer boundary will be replaced by Bézier element. First, the exterior elements will be identified by using FE mesh connectivity. Then, the identified elements will be replaced by Bézier elements through the relationship between the physical and parametric spaces. Regarding the conversion, the facial nodes will be converted into Bézier control points, which will be used in the remaining analysis.

2.1 Bernstein–Bézier finite element

Before presenting the conversion procedure, a few preliminaries about Bézier element will be introduced. The basis functions will be presented in Bernstein–Bézier form (B-form) with several appropriate properties for the numerical analysis [19]. Then, the quadratic triangular and tetrahedral Bézier elements will be presented with the corresponding control net.

2.1.1 Bernstein polynomials

Bernstein polynomials are used to construct Bézier curve and surface. The relevant univariate, bivariate, and trivariate forms are expressed in Eqs. 1a to 1c, respectively.

$$\begin{aligned}&B_{ij}^d\left( \lambda _1\right) =\frac{d!}{i!j!} \lambda _1^i\left( 1-\lambda _1\right) ^j, \quad i+j=d \end{aligned}$$
(1a)
$$\begin{aligned}&B_{ijk}^d\left( \lambda _1,\lambda _2\right) =\frac{d!}{i!j!k!} \lambda _1^i\lambda _2^j\left( 1-\lambda _1-\lambda _2\right) ^k,\nonumber \\&\qquad \qquad \qquad \qquad i+j+k=d \end{aligned}$$
(1b)
$$\begin{aligned}&B_{ijkl}^d\left( \lambda _1,\lambda _2,\lambda _3\right) =\frac{d!}{i!j!k!l!} \lambda _1^i\lambda _2^j\lambda _3^k \left( 1-\lambda _1-\lambda _2-\lambda _3\right) ^l,\nonumber \\&\quad \qquad \qquad \qquad \qquad \quad i+j+k+l=d \end{aligned}$$
(1c)

Herein, B is Bernstein polynomials, and d is the degree of Bernstein polynomials. \(\lambda _1\), \(\lambda _2\), \(\lambda _3\) are the parametric coordinates.

B-spline and NURBS, which are generally utilized in CAD, are the generalized version of Bernstein polynomials [20]. The basis functions are represented using Cox–deBoor recursion as in Eq. 2.

$$\begin{aligned}&N_{i}^0\left( r\right) ={\left\{ \begin{array}{ll} 1 & \text {if } r_i\le r < r_{i+1} \\ 0 & \text {otherwise} \end{array}\right. } \end{aligned}$$
(2a)
$$\begin{aligned}&\begin{aligned}&N_{i}^d\left( r\right) =\frac{r-r_i}{r_{i+d}-r_{i}}N_{i}^{d-1}\left( r\right) \\&\quad +\frac{r_{i+d+1}-r}{r_{i+d+1}-r_{i+1}}N_{i+1}^{d-1}\left( r\right) \end{aligned} \end{aligned}$$
(2b)
$$\begin{aligned}&R_{i}^d\left( r\right) =\frac{N_{i}^{d}\left( r\right) w_i}{\sum _{i=1}^{n+d+1}N_{i}^{d}\left( r\right) w_i} \end{aligned}$$
(2c)

Herein, r is the parametric coordinates. w is the weight for each control point. N and R are the basis functions of B-spline and NURBS, respectively.

The curves and surfaces expressed via the above basis functions have a relationship with each other. As expressed in Eq. 2, B-spline is the non-rational form of NURBS. In other words, the B-spline is equivalent to NURBS with a weight of unity for all the control points. Also, Bézier curve and surface, which are constructed via Bernstein polynomials, correspond to each knot span of the B-spline curve and surface. This suggests that an entire domain constructed with Bézier surfaces will be identical to NURBS surfaces with a weight of unity. Despite employing the non-rational form of NURBS, the representation of a surface using Bernstein polynomials will be more accurate than that using a simple polygonal element. Besides, the Bézier element contains a strong convex hull property, and its basis function is formulated on the parametric space while satisfying the partition of unity. Furthermore, because Bernstein polynomials are element-wisely defined, they will be easily adopted into the conventional FEA approach.

2.1.2 Bézier elements

Through the conversion from FE to Bézier elements, the exterior triangular surfaces will be replaced by Bézier triangles. All the elements included in the interiors are assumed to be Bézier tetrahedron. In this work, a quadratic element, which has been frequently used in the conventional FEA, is used. Bézier elements in the physical space will be represented by a linear combination of Bernstein basis and control point. Equations 3 and 4 define the physical point \(x^{phy}\) of Bézier triangle and tetrahedron, respectively.

$$\begin{aligned}&x^{phy}=\sum _{i+j+k=d}B_{ijk}\left( {\lambda _1,\lambda _2}\right) p_{ijk} \end{aligned}$$
(3)
$$\begin{aligned}&{x}^{phy}=\sum _{i+j+k+l=d}B_{ijkl}\left( {\lambda _1,\lambda _2,\lambda _3}\right) p_{ijkl} \end{aligned}$$
(4)

Herein, p is the control point with i, j, k, l \(\in \ [0,d]\).

The quadratic Bézier elements in the parametric space with the numbered control points are depicted in Fig. 1.

Fig. 1
figure 1

Quadratic Bézier a triangular and b tetrahedral elements with corresponding control points in parametric space

Full size image

2.1.3 Derivatives in a B-form

The derivatives of the Bernstein polynomial will be required for Hermite interpolation. The directional derivative at point \(\varvec{\lambda } = \left( \lambda _1,\lambda _2,\lambda _3\right) ^\top \) with respect to direction \(\textbf{u}\) is defined in Eq. 5.

$$\begin{aligned} D_\textbf{u}f\left( \varvec{\lambda }\right) = d\sum _{i+j+k+l=1}c_{ijkl}^{\left( 1\right) }\left( \textbf{a}\right) B_{ijkl}^{d-1}\left( \varvec{\lambda }\right) \end{aligned}$$
(5)

Herein, f is a polynomial function. Direction vector \(\textbf{u}\) is defined in terms of the barycentric coordinates \(\textbf{a}=\left( a_1,a_2,a_3,a_4\right) ^\top \). The quantity \(c_{ijkl}^{\left( 1\right) }\left( \textbf{a}\right) \) will be obtained from the first step of the de Casteljau algorithm as expressed in Eq. 6.

$$\begin{aligned} \begin{aligned} c_{ijkl}^{\left( 1\right) }\left( \textbf{a}\right)&=c_{i+1,j,k,l} a_1+c_{i,j+1,k,l} a_2\\&\quad +c_{i,j,k+1,l}a_3+c_{i,j,k,l+1}{a}_4,\\&\quad i+j+k+l=d-1 \end{aligned} \end{aligned}$$
(6)

Herein, \(c_{ijkl}\) is the control point.

2.2 Bézier mesh generation

As mentioned in the previous section, FE discretization will be used to obtain the inner volumetric information. Because it uses a polygon to represent an object, geometric discrepancy will occur at the exterior curved surfaces. Thus, the conversion from FE to the Bézier element will be required to represent a curved surface with a reduced geometric discrepancy. First, an FE located on the exterior surface will be identified. Then, the FE and the relevant nodes will be converted into a Bézier element and corresponding control points. The element conversion at the boundary parametric space is based on the relationship between the physical point of the FE mesh and that of the Bézier surface. The physical points of the FE mesh will be the same as the linear combination of Bézier control points and Bézier basis at the boundary parametric coordinate. For a single tetrahedron, such relationship is expressed by Eq. 7. Then, Bézier control points can be obtained using Eq. 8 for each element. The matrix form is expressed as Eq. 9.

$$\begin{aligned}&{x}_j^{phy}=\sum _{i=1}^{n_f}{B}_{i}^d\left( {\varvec{\lambda }_j}\right) {p}_{i} ,\quad j = 1, ..., n_s \end{aligned}$$
(7)
$$\begin{aligned}&\textbf{P}=\mathbf {{B}^{-1}}\left( {\varvec{\lambda }}\right) {\textbf{P}^{phy}} \end{aligned}$$
(8)
$$\begin{aligned} \textbf{B}({\varvec{\lambda }})=\begin{bmatrix}{B}_{1}^{d}({\varvec{\lambda }}_{1})& \dots & {B}_{n_s}^{d}({\varvec{\lambda }}_{1})\\ \vdots & \ddots & \vdots \\ {B}_{1}^{d}({\varvec{\lambda }}_{n_s})& \dots & {B}_{n_f}^{d}({\varvec{\lambda }}_{n_s})\end{bmatrix} \end{aligned}$$
(9a)
$$\begin{aligned} {\textbf{P}^{phy}}=\begin{Bmatrix}{x}_1^{phy}\\\vdots \\{x}_{n_s}^{phy}\end{Bmatrix} ,\quad \textbf{P}=\begin{Bmatrix}{p}_1\\\vdots \\{p}_{n_s}\end{Bmatrix} \end{aligned}$$
(9b)

Herein, \({n_f}\) and \({n_s}\) are the number of boundary elements and nodes per face, respectively. \(\textbf{B}\) is the matrix of Bernstein polynomials at the parametric coordinate \(\varvec{\lambda }\), \(\textbf{P}^{phy}\) is a physical point, and \(\textbf{P}\) is the Bézier control point. Figure 2 shows a two-dimensional curved configuration discretized by Bézier elements.

Fig. 2
figure 2

a Two-dimensional curved CAD configuration and b discretization by the Bézier element. The blue dot indicates the location of control points

Full size image

As other methods for the surface element replacement, Bézier extraction [21] and Bézier projection [22] may be used. However, those are the local operators and have a restriction that each face of Bézier elements should be in a single NURBS knot span. Therefore, those methodologies will not be compatible with an existing mesh generator. Those will further require a cumbersome process to be adopted in a complicated trimmed geometry. Oppositely, the present method will be easily applicable to complicated geometries with the compatibility to an existing mesh generator.

Even when applying the element conversion, the exterior surface expressed by the non-rational Bézier elements may not be identical to an exact CAD geometry. However, as the number of non-rational elements comprised in a single surface increases, the geometric discrepancy against the rational surface will drastically decrease. The situation is visualized in Fig. 3. As the number of the elements increases, a curve discretized by Bézier elements will show higher improvement in terms of the geometric discrepancy than that discretized by the polygonal elements.

Fig. 3
figure 3

Approximation of the surface of a unit quarter circle. The number of elements is a 1, b 2, and c 4

Full size image

3 Continuity enhancement

Although the Bézier elements represent the object with improved geometric fidelity, the basis function for numerical approximation will still be \(C^0\) inter-element continuous. The higher-order continuity of the basis function contributes to predicting physical phenomena more accurately [23]. This section introduces the enhancement of the continuity condition of the domain composed of \(C^0\) continuous Bézier elements. In this study, a higher-order continuous Bézier basis function will be obtained through Hermite interpolation on the parts where Bézier tetrahedrons meet. For improved compatibility and computational efficiency, a macro element technique with relaxed constraints will be used. Then, the approximate \(C^1\) basis function will be constructed on the minimal determining set (MDS).

3.1 Macro element split

Constructing the \(C^1\) basis function requires one greater than 9th-order tetrahedral element for the three-dimensional analysis. Meanwhile, by splitting a single original tetrahedron into multiple tetrahedrons, the basis function of the higher continuity will be obtained with the lower-order tetrahedral element. Various macro element techniques have been utilized to efficiently construct the basis function of higher inter-element continuity. Powell-Sabin split and Clough-Tougher split were utilized for two-dimensional analysis with the 2nd-order and 3rd-order triangles, respectively [17, 24]. For three-dimensional analysis, Alfeld split was used with the 5th-order tetrahedron [18]. Conventionally, Worsey–Piper, Worsey-Farin, and Alfeld split techniques have been employed to construct the \(C^1\)-continuous basis function of the 2nd, 3rd, and 5th-order tetrahedrons [25].

In this study, the construction of a higher-order continuous basis function will be attempted using the 2nd-order tetrahedral element, which has been frequently used in traditional FEA. Worsey–Piper technique in Fig. 4 will be employed, which divides a single macro tetrahedron into 24 micro tetrahedrons for the construction of \(C^1\) basis function [26]. The macro vertex is a node of an original macro tetrahedron. Edge, facial, and tetrahedral split point are the node of micro tetrahedrons, which divide macro edge, face, and tetrahedron, respectively. For the exact split, the split point should satisfy the two requirements. First, the facial split point should lie on the line connecting tetrahedral split points. Second, tetrahedral split points around the edge should be co-planar, and the edge split point should lie on the plane. To satisfy those conditions, entire facial angles of each tetrahedron should be acute. However, such highly restricting conditions will suppress the applicability for various geometry. Therefore, an approximated Worsey–Piper split will be developed to alleviate the split condition related to the edge split point.

Fig. 4
figure 4

Nomenclature of split points of Worsey–Piper split technique

Full size image

3.2 Approximate \(C^1\) basis function

To construct a higher-order continuous basis function, the inter-element continuity condition will be applied to each micro tetrahedron. It is quoted as Hermite interpolation, for which the formulation is described in Eqs. 10a to 10d.

$$\begin{aligned}&s(x_i,y_i,z_i)=k_i,\quad i = 1,\cdots ,~n \end{aligned}$$
(10a)
$$\begin{aligned}&D_xs(x_i,y_i,z_i)=k_i^x,\quad i = 1,\cdots ,~n \end{aligned}$$
(10b)
$$\begin{aligned}&D_ys(x_i,y_i,z_i)=k_i^y,\quad i = 1,\cdots ,~n \end{aligned}$$
(10c)
$$\begin{aligned}&D_zs(x_i,y_i,z_i)=k_i^z,\quad i = 1,\cdots ,~n \end{aligned}$$
(10d)

Herein, \(s(x_i,y_i,z_i)\) is a \(C^0\) Bézier spline at the intersection point. And, \(D_xs\), \(D_ys\), and \(D_zs\) are its x, y, and z-directional derivatives. n is the number of points to be interpolated.

Such directional derivatives will be used at the interfacial components, among the tetrahedrons such as a node, edge, and face. Consequently, the control properties of one tetrahedron will be determined from the others. An example of a directional derivative related to two tetrahedrons sharing the same node is presented in Fig. 5. Then, Bézier ordinates in one tetrahedron will be defined as the linear combination of those in another tetrahedron through their directional derivatives based on the continuity condition. A similar relationship will be established for the interfacial components of the other types. Using the Bernstein–Bézier technique, the derivatives of the basis function will be obtained in terms of the barycentric coordinates as in Eqs. 5 and 6.

Fig. 5
figure 5

Two tetrahedrons sharing the same node and directional derivatives for the continuity condition. The blue arrow indicates the directional derivative vector

Full size image

Through enforcing the continuity, some set \(\Gamma \subseteq \textbf{S}_d^0\) of the total domain \(\textbf{S}_d^0\) that determine all other set \(\Gamma ^c=\textbf{S}_d^0-\Gamma \) will be defined. Then, \(\Gamma \) will be the determining set for \(\textbf{S}_d^0\), and the smallest set of \(\Gamma \) will be the MDS for \(\textbf{S}_d^0\). Therefore, the entire control points \(\textbf{c}\subseteq \textbf{S}_d^0\) will be expressed as a linear combination of MDS with the continuity coefficients as described in Eq. 11.

$$\begin{aligned}&\textbf{c}=\textbf{A}\tilde{\textbf{c}} \end{aligned}$$
(11a)
$$\begin{aligned}&\textbf{A}=\begin{bmatrix}\textbf{C}\\ \textbf{I}\end{bmatrix} \end{aligned}$$
(11b)

Herein, \(\textbf{c}\subseteq {\textbf{S}}_d^0\) and \(\tilde{\textbf{c}}\subseteq \mathcal {M}\subseteq {\textbf{S}}_d^0\) are the entire and MDS control points, in which \(\mathcal {M}\) denotes MDS. \(\textbf{C}\) is a continuity coefficient matrix.

In other words, the entire control points will be determined by employing MDS only, and it will be the unknown degree of freedom (DOF) to be obtained for the present framework. A detailed explanation for the selection of MDS and its validation are found in [25,26,27]. One example of a proper MDS for a two-dimensional example geometry that conducts the Powell-Sabin split, two-dimensional form of the Worsey–Piper technique, is depicted in Fig. 6. In the Powell-Sabin split, entire domain points will be determined using the macro vertex and its surrounding two points [25].

Fig. 6
figure 6

Two-dimensional schematics of the macro element split. Discretization based on a macro and b micro elements

Full size image

Then, a higher-order continuous basis function will be constructed using the continuity coefficient matrix. In the present approximated Worsey–Piper split, the continuity at the edge split point of a macro tetrahedron will not be the exact \(C^1\) continuity. Despite such a weakness, a moderate accuracy enhancement against the conventional FEA will be shown in the later section.

The approximate \(C^1\) basis function for each Bézier tetrahedral element will be defined on MDS control points included in and adjacent to the tetrahedron. This is obtained from the linear combination of the \(C^0\) Bézier basis function at the parametric point with the continuity coefficients, as described in Eq. 12.

$$\begin{aligned} \psi _{\xi \mid {T}}\left( \hat{\textbf{x}}\right) =\sum _{\textbf{c}_\eta \in \textbf{S}_{d,T}^0}a_{\eta \xi }B_\eta \left( \hat{\textbf{x}}\right) \end{aligned}$$
(12)

Herein, \(\psi _{\xi \mid {T}}\) is the approximate \(C^1\) basis of MDS point \(\xi \) associated with the tetrahedron T, and \(B_\eta \) is the \(C^0\) Bézier basis at the parametric point \(\textbf{c}_\eta \). \(\textbf{S}_{d,T}^0\subseteq \textbf{S}_{d}^0\) is the set of entire points included in the tetrahedron T. \(\hat{\textbf{x}}\) is the arbitrary parametric location in tetrahedron T. \(a_{\eta \xi }\) is the coefficient of the matrix \(\textbf{A}\) in Eq. 11. \(\xi \) and \(\eta \) are the indices for the points in MDS, and those in the set of the entire domain, respectively.

The constructed approximate \(C^1\) basis function is unity for each \(\tilde{\textbf{c}}_{\xi }\in \mathcal {M}\) as expressed in Eq. 13.

$$\begin{aligned} \gamma _{\eta }\psi _{\xi }=\delta _{\eta \xi },\quad \textbf{c}_\eta \in \mathcal {M} \end{aligned}$$
(13)

Herein, \(\gamma _\eta \) is the linear functional that selects the coefficient associated with the point \(\textbf{c}_\eta \in {\textbf{S}}_d^0\). \(\delta \) is a Dirac delta function.

The continuity coefficient for the MDS point is unity only on the corresponding points, and it is zero for all the other components. Consequently, the basis function satisfies the property of partition of unity, as well as that of the local support. The constructed basis function will be used in traditional FE formulation easily because it satisfies the condition for shape function.

4 Analysis procedure with an approximate \(C^1\) Bernstein–Bézier finite element

In this section, the formulation for structural analysis will be introduced. Sequentially, the comprehensive framework of the presented methodology will be explained in detail, in the context of a linear elastic structure. In addition, it will be compared against the analysis procedures of conventional FEA.

4.1 Structural analysis formulation

In this section, we aim to provide a concise formulation of the subsequent linear structural analysis that will be applied. First, the governing equation and boundary condition for linear elasticity are expressed in Eq. 14.

$$\begin{aligned}&\nabla \cdot \varvec{\sigma }+\textbf{b}_{f} = 0,\quad \text {on} \,\quad \Omega \end{aligned}$$
(14a)
$$\begin{aligned}&\varvec{\sigma }\cdot \textbf{n}=\textbf{t}_{f},\quad \text {on} \,\quad \partial \Omega _{t} \end{aligned}$$
(14b)
$$\begin{aligned}&\textbf{u}=\bar{\textbf{u}},\quad \text {on} \,\quad \partial \Omega _{u} \end{aligned}$$
(14c)

Herein, \(\Omega \) is the entire domain. \(\partial \Omega _{t}\) and \(\partial \Omega _{u}\) are the surfaces imposed by the Neumann and Dirichlet boundary condition, respectively. \(\textbf{b}_f\) is a body force and \(\textbf{t}_f\) is surface traction with traction vector \(\textbf{n}\).

Using the approximate \(C^1\) basis function, which is defined on the MDS points, arbitrary physical position \(\textbf{x}\) will be approximated by Eq. 15.

$$\begin{aligned} \textbf{x}=\sum _{\tilde{\textbf{c}}_{\xi }\in \mathcal {M}}\tilde{\textbf{c}}_{\xi }\psi _{\xi } = \sum _{\tilde{\textbf{c}}_{\xi }\in \mathcal {M}_{b}}\tilde{\textbf{c}}_{\xi }\psi _{\xi } + \sum _{\tilde{\textbf{c}}_{\xi }\in \mathcal {M}_0}\tilde{\textbf{c}}_{\xi }\phi _{\xi } \end{aligned}$$
(15)

Herein, \(\mathcal {M}_{0}\) is MDS points corresponding to \(\partial \Omega _{u}\). The rest of the MDS set will be defined as \(\mathcal {M}_{b} = \mathcal {M} - \mathcal {M}_{0}\). \(\phi _\xi \) is the \(C^0\) basis function.

The material and structural properties will be defined in Eq. 16.

$$\begin{aligned}&\tilde{\textbf{K}}_e = \textbf{A}_{e}^{\top }\left( \sum _{i=1}^{24}\textbf{L}_{ie}^{\top }\int _{\Omega _{ie}}\textbf{B}_{mat}^{\top }\textbf{D}\textbf{B}_{mat} d\Omega \textbf{L}_{ie}\right) \textbf{A}_{e} \end{aligned}$$
(16a)
$$\begin{aligned}&\tilde{\textbf{M}}_e = \textbf{A}_{e}^{\top }\left( \sum _{i=1}^{24}\textbf{L}_{ie}^{\top }\int _{\Omega _{ie}}\rho \textbf{N}_{mat}^{\top }\textbf{N}_{mat} d\Omega \textbf{L}_{ie}\right) \textbf{A}_{e} \end{aligned}$$
(16b)

Herein, \(\tilde{\textbf{K}}_e\) and \(\tilde{\textbf{M}}_e\) are the elemental stiffness and mass matrices of the e-th macro element. \(\tilde{\square }\) denotes that the quantity \(\square \) is defined in terms of MDS points. \(\Omega _{ie}\) is the volume of i-th micro element in the e-th macro element. In addition, \(\textbf{B}_{mat}\) and \(\textbf{N}_{mat}\) are the matrix of the strain–displacement relationship and shape function constructed for \(C^0\) continuity, respectively. \(\textbf{D}\) is the linear elastic constitutive matrix. \(\textbf{L}_{ie}\) is the assembling Boolean matrix from the i-th micro element to the e-th macro element. \(\textbf{A}_{e}\) is the elemental continuity matrix of the e-th macro element.

Fig. 7
figure 7

Graphical representation of the present analysis procedures

Full size image

Sequentially, the numerical examples will be examined in Sect. 5, in which the rotating effect is included in Sect. 5.3. The relevant formulation is written in advance in Eqs. 17 and 18 for the steady-state and frequency response analyses, respectively.

$$\begin{aligned} \left( \tilde{\textbf{K}}-\tilde{\textbf{K}}^{sp}\right) \tilde{\textbf{u}} = \tilde{\textbf{f}} \end{aligned}$$
(17)
$$\begin{aligned} \left( -\omega ^2\tilde{\textbf{M}} + \tilde{\textbf{K}} - \tilde{\textbf{K}}^{sp}\right) \tilde{\textbf{u}}_{hr} = \tilde{\textbf{f}}_{hr} \end{aligned}$$
(18)

Herein, \(\tilde{\textbf{K}}\) is the global stiffness matrix. \(\tilde{\textbf{K}}^{sp}\) is the global spin-softening matrix due to the rotation [28]. The spin-softening relevance matrix will be included for the case of a rotating component. \(\tilde{\textbf{u}}\) and \(\tilde{\textbf{f}}\) are the displacement and external force vector, respectively. For the frequency response analysis in Eq. 18, damping effect is ignored. \(\tilde{\textbf{M}}\) is the global mass matrix. \(\tilde{\textbf{u}}_{hr}\) and \(\tilde{\textbf{f}}_{hr}\) are the amplitude vectors of displacement and external excitation load, respectively.

4.2 Procedure of the present computation

Figure 7 represents the analysis procedure of the present methodology for a two-dimensional linear elastic example. The main features of the present method include the reduction of the geometric discrepancy and the utilization of higher continuity for shape function.

An existing FE mesh generator will discretize the CAD geometry by polygonal FEs. Such discretization will bring geometric discrepancy, particularly on the curved surface. To reduce it, the polygonal elements will be converted into \(C^0\) Bézier elements. The resulting Bézier elements will exhibit significantly diminished geometric discrepancy compared against the original FE discretization when representing a curved geometry.

Subsequently, an inter-element continuity of the basis function will be enhanced through the element split technique on the parametric space. The FE discretized shape will serve as parametric space. This paper will utilize the approximated Worsey–Piper split technique to alleviate the restraints of the original Worsey–Piper split. Through Hermite interpolation for all the micro elements, the basis function with enhanced continuity will be defined on MDS, and those will be employed to represent the entire domain.

After the pre-processing, the present method will be conducted as illustrated in Fig. 8. The formulation of the present method is similar to the conventional FEA, except that the created solution is primarily associated on the MDS. In pursuance of the continuity enhancement, the relationships between the Bézier ordinates of each micro element and those of a specific MDS will be considered. From the local support property of the basis function, such MDS will be adjacent to the corresponding micro element. Based on the relationship, displacements for the entire domain will be approximated by using the values which defined only on the MDS. The material and structural properties will be obtained through the Gaussian quadrature and assemblage of the elemental quantities which defined on the corresponding MDS. By using the obtained displacements on the MDS and continuity coefficients, nodal stress values will be acquired directly without any additional techniques such as the extrapolation. However, a spurious wriggling will appear because of the approximate \(C^1\) continuity. Following the numerical processing, the obtained solutions will be approximated for the entire physical domain based on the continuity coefficients as well.

Fig. 8
figure 8

Overview of the present procedure

Full size image
Fig. 9
figure 9

Comparison of the algorithms of a the conventional FEA, and b the present method

Full size image

4.3 Comparison of the analysis procedures

To describe the capabilities of the present method clearly, the algorithms for the analysis are to be compared with conventional FEA as depicted in Fig. 9. Herein, \(n_e\) and \(n_e^m\) is the numbers of elements and macro elements. \(\textbf{L}_e\) and \(\tilde{\textbf{L}}_e\) are the global assemblage matrix from e-th element. \(\Omega _e\) is the volume of e-th element. In this paper, a sparse direct solver in MATLAB will be utilized for the linear algebraic solver. The conventional FEA method employs displacement formulation with Lagrange polynomial interpolation [29] and 10-node tetrahedral element.

In the conventional FEA, CAD geometry will be discretized using 10-node tetrahedral FEs. From the Gaussian quadrature and assemblage of each element, the entire stiffness matrix will be constructed. Then, the nodal displacement values will be estimated by using a linear algebraic solver. After that, stress recovery techniques which are required to obtain the values from the nodal displacement will be conducted such as extrapolation and nodal averaging.

From the same input as the conventional FEA, the present method will use Bézier mesh generation to reduce a geometric discrepancy induced from the polygonal representation. Furthermore, the approximate \(C^1\) basis function will be constructed through the element split and the MDS with the continuity coefficients. From the Gaussian quadrature each elemental stiffness matrix will be constructed with corresponding MDS. Beyond an assemblage of the matrices, displacements on the control points will be predicted by using a linear algebraic solver as used in the FEA. From the higher-order continuity of the basis function, nodal stress will be obtained directly.

5 Numerical applications

In this section, the present approach will be applied to the analysis of the curved geometry, thereby validating its advantages against the conventional FEA. Initially, in Sects. 5.1 and 5.2, a linear elasticity analysis will be implemented on curved geometries with exact solutions. A comprehensive comparison between the numerical results by the present method and conventional FEA will be attempted regarding the varying mesh size relative to the analytical solution. The convergence of the results in terms of mesh refinement will be discussed by employing \(H_0\)-norm error of the stress [30] as a metric. In this case, the accuracy of the present method and conventional FEA will be compared on an identical mesh obtained by the micro elements. The \(H_0\)-norm error is defined as

$$\begin{aligned} \epsilon =\sqrt{\sum _{e=1}^{n_e}\epsilon _e^2} \end{aligned}$$
(19)

where the contribution of the e-th element, \(\epsilon _e\), is given as follows:

$$\begin{aligned} \epsilon _e=\sqrt{\frac{\int _{\Omega _e}\left( \tilde{\varvec{\sigma }}-\tilde{\varvec{\sigma }}_{ref}\right) ^\top \left( \tilde{\varvec{\sigma }}-\tilde{\varvec{\sigma }}_{ref}\right) d\Omega }{\int _{\Omega _e}\tilde{\varvec{\sigma }}_{ref}^\top \tilde{\varvec{\sigma }}_{ref} d\Omega }}. \end{aligned}$$
(20)

Herein, \(\tilde{\varvec{\sigma }}\) and \(\tilde{\varvec{\sigma }}_{ref}\) are the stress by the numerical prediction and analytical solution, respectively, along with Voigt notation.

Then, the computational cost will be evaluated. In this case, different sets of discretization will be employed for the present method and conventional FEA, while producing similar accuracy on each other. The computational time will be measured under MATLAB R2023b environment with a single processor of Intel i9-13900K 3.00 GHz CPU.

Furthermore, in Sect. 5.3, the current method and the conventional FEA will be applied to an intricate turbine blade geometry under the computational routine of Fig. 9, demonstrating the geometric applicability of the present method. Moreover, to demonstrate the adaptability of the present method for cumbersome simulation, steady-state and frequency response analysis including the rotational consideration will be attempted.

5.1 Inner pressurized thick-walled cylinder

As the first example to demonstrate the accuracy of the present approach, a linear elastic cylindrical object is selected. From the symmetry of the geometry, its one-quarter division will only be considered as in Fig. 10. Elastic modulus and Poisson’s ratio are 1Pa and 0, respectively. Uniform pressure is applied on the internal surface. Plane strain conditions are applied by fixing the axial direction. An analytic solution for this situation is presented in [1].

Fig. 10
figure 10

Geometry of inner pressurized thick-walled cylinder: a 1/4 configuration and b von-Mises stress distribution of analytic solution

Full size image
Fig. 11
figure 11

Inner pressurized thick-walled cylinder: numerical discretization of Case 1 with a macro element and b micro element

Full size image
Table 1 Inner pressurized thick-walled cylinder: mesh information for the analysis
Full size table

In Fig. 11, a graphical representation of discretization with macro and micro elements is presented. In this example, both conventional FEA and the present approach are attempted on an identical discretization based on micro elements. Detailed mesh information is summarized in Table 1. The convergence of the computation in comparison with the analytic solution, under varying levels of mesh refinement, is plotted in Fig. 12. \(H_0\)-norm error distribution of stress for each discretization is plotted in Figs. 13 and 14. It is found that the present approach demonstrates enhanced accuracy and convergence rate.

Fig. 12
figure 12

Inner pressurized thick-walled cylinder: convergence comparison between the present method and conventional FEA

Full size image
Fig. 13
figure 13

Inner pressurized thick-walled cylinder: stress \(H_0\)-norm error distribution of the present analysis for each discretization

Full size image
Fig. 14
figure 14

Inner pressurized thick-walled cylinder: stress \(H_0\)-norm error distribution of the conventional FEA for each discretization

Full size image

Table 2 summarizes the computational cost by the present method when it provides a similar accuracy to that of the conventional FEA. As shown, the present method is 15 times faster than the conventional FEA despite additional computational procedures such as the Bézier mesh generation. Such result demonstrates that the present method will alleviate the requirement of FE discretization while producing a relatively accurate solution.

Table 2 Inner pressurized thick-walled cylinder: computational time comparison
Full size table

5.2 Inner pressurized hollow sphere

A linear elastic hollow spherical example is selected for the second geometry to validate the accuracy of the present method against the conventional FEA. To alleviate the expense of the analysis, its one-eighth division will only be considered as illustrated in Fig. 15. Elastic modulus and Poisson’s ratio are 1Pa and 0.3, respectively. Uniform pressure is applied on the internal surface. The analytic solution for this situation is presented in [31].

Fig. 15
figure 15

Geometry of inner pressurized hollow sphere: a 1/8 configuration and b von-Mises stress distribution of analytic solution

Full size image

In Fig. 16 and Table 3, graphical representation and detailed mesh information are provided. Both conventional FEA and the present approach are attempted on an identical discretization which will be the same as the previous example. Stress predictions for both analyses against the analytic solution are plotted in Fig. 17. \(H_0\)-norm error distribution of stress for each discretization is plotted in Figs. 18 and 19.

In Figs. 13 and 18, the present approach shows partially scattered \(H_0\)-norm error distribution. Its main reasons will be degraded mesh quality due to the element split and the approximations in continuity. However, the convergence and accuracy of the predictions by the present approach are found to be superior to those of the conventional FEA on an overall domain.

Fig. 16
figure 16

Inner pressurized hollow sphere: numerical discretization of Case 1 with a macro element and b micro element

Full size image
Table 3 Inner pressurized hollow sphere: mesh information for the analysis
Full size table
Fig. 17
figure 17

Inner pressurized hollow sphere: convergence comparison between the present method and conventional FEA

Full size image
Fig. 18
figure 18

Inner pressurized hollow sphere: stress \(H_0\)-norm error distribution of the present analysis for each discretization

Full size image
Fig. 19
figure 19

Inner pressurized hollow sphere: stress \(H_0\)-norm error distribution of the conventional FEA for each discretization

Full size image

As in the previous example, the comparison of the computational cost is summarized in Table 4. The speed-up capability of the present method is found to be 26, which is even larger than that in Sect. 5.1. Therefore, it is inferred that the present method becomes more efficient as the number of unknowns is increased.

Table 4 Inner pressurized hollow sphere: computational time comparison
Full size table

5.3 GT11N turbine blade

Finally, a rotating turbine blade will be selected. The hardware shown in Fig. 20 is the first stage turbine blade of a 75MW GT11N gas turbine. The complexity of the configuration, including multiple inner holes and curved surfaces due to internal coolant channels, enforces discretization for NURBS-based IGA to become cumbersome. Also, prediction of its mechanical behavior using the conventional FEA with high-fidelity may require a large number of DOFs due to the fine discretization, as shown in Sects. 5.1 and 5.2. The primary emphasis of the application is on demonstrating the geometrical applicability and adaptability to more realistic handling capability of the present method.

Fig. 20
figure 20

Representation of the \(1^{st}\) stage turbine blade of a 75MW GT11N by a CAD with opaque surfaces, b CAD with transparent surfaces, and c FE discretization

Full size image

First, steady-state response of the blade under the centrifugal force will be compared. The rotating speed is 3,600 RPM with a fixed boundary condition as shown in Fig. 21, and the material properties of Alloy In-738LC are employed [32]. The blade will be discretized into 170,236 10-node tetrahedral elements and 278,310 FE nodes.

von-Mises stress distribution of the turbine blade is depicted in Fig. 22 and Fig. 23. The red box locations of the figures represent the maximum differences between the results. Although a spurious stress wriggling appears owing to the approximate \(C^1\) continuity, a large stress gradient on the dovetail and slot structures is predicted more clearly using the present method. On the contrary, the result by the conventional FEA shows a scattered stress gradient due to Gaussian extrapolation and nodal averaging.

Fig. 21
figure 21

GT11N turbine blade: analysis condition for steady-state response prediction

Full size image
Fig. 22
figure 22

GT11N turbine blade: isometric view of von-Mises stress distribution by a conventional FEA and b the present method for steady-state response prediction

Full size image
Fig. 23
figure 23

GT11N turbine blade: magnified view on the trailing edge slot of von-Mises stress distribution by a conventional FEA and b the present method for steady-state response prediction

Full size image

Sequentially, the present method will be used to execute the frequency response analysis, and the relevant results will be compared against those by the FEA as in Fig. 24. Except the transverse harmonic excitation on the tip of \(1\cos {\left( 2\pi f_{hr} t\right) }\) kN with frequency \(f_{hr}\) of 1 kHz, other conditions will be the same for the steady-state analysis. Herein, t denotes time. In Figs. 25 and 26, the predicted responses are depicted. Similar to the previous results, spurious stress wriggling was observed in the present analysis. However, large stress gradient near the blade root was predicted much more significantly than the conventional FEA.

Fig. 24
figure 24

GT11N turbine blade: analysis condition for frequency response prediction

Full size image
Fig. 25
figure 25

GT11N turbine blade: isometric view of von-Mises stress distribution at \(t=0\) of the entire geometry by a conventional FEA and b the present method for frequency response prediction

Full size image
Fig. 26
figure 26

GT11N turbine blade: rear view of von-Mises stress distribution at \(t=0\) by a conventional FEA and b the present method for frequency response prediction

Full size image

In the comparison of stress distributions for the turbine blade geometry, the difference will be investigated regarding the stress wriggling and stress gradient using the present method. In the process for continuity enhancement, the present method uses the element split technique by relieving the requirement for original Worsey–Piper split. Therefore, the continuity of the points related with the edge split point as in Fig. 4 will not become an exact \(C^1\). Consequently, the obtained values in that point exhibit discrepancies compared with the others, manifesting as wriggling in the stress distribution.

The numerical approximation of the present method retains a higher level of continuity than the conventional FEA. The utilization of basis functions, which maintain high inter-element continuity for the approximation of displacement is advantageous in obtaining a more accurate stress prediction. Furthermore, while FEA is a well-established technique, it may require a large number of DOFs due to the fine discretization to predict stress distribution for the complicated curved geometry. Consequently, when predicting by the identical discretization, the present method will show superiority under a large stress gradient.

6 Conclusion

In this study, an improved three-dimensional structural analysis framework, emphasizing in accuracy and applicability, was presented. The proposed approach attained the advantageous features of FE-based volumetric IGA through the pre-processing procedures. Notably, the method is seamlessly compatible with the conventional FEA input, which is based on a 10-node tetrahedron, using a \(2^{\textrm{nd}}\)-order Bernstein–Bézier basis function. Through Bézier mesh generation, one of the pre-processing procedures, geometric accuracy was enhanced by approximated geometrical mapping without CAD information. Additionally, an approximate \(C^1\) continuity of the basis function was constructed via an approximated Worsey–Piper split, which alleviates the restriction of the original Worsey–Piper split. Although the accuracy of the presented method is not as precise as that of IGA, it surpasses that of the conventional FEA.

To demonstrate the features, the presented method was applied to various numerical applications. After validating the method by comparing it with the FEA against the analytic solution, the method was implemented to complicated configuration. In comparing the results with FEA, the present method exhibits spurious wriggling due to the approximate \(C^1\) continuity. However, the method not only excels conventional FEA in numerical accuracy on an overall domain but predicts the large stress gradient that the FEA does not.

In the future, study on the optimal convergence rate of approximate \(C^1\) continuity will be conducted. Additionally, the implementation will be extended to more intricate simulations such as elasto-plastic thermal analyses.