Numerical Aspects of Isogeometric Boundary Element Methods: (Nearly) Singular Quadrature, Trimmed NURBS and Surface Crack Modeling

Xuan Peng; Haojie Lian

doi:10.32604/cmes.2022.017410

1Department of Mechanical Engineering, Suzhou University of Science and Technology, Suzhou, 215009, China
2Key Laboratory of In-situ Property-improving Mining of Ministry of Education, Taiyuan University of Technology, Taiyuan, 030024, China
*Corresponding Authors: Xuan Peng. Email: pengxuan89@sina.com; Haojie Lian. Email: lianhaojie@gmail.com
Received: 08 May 2021; Accepted: 29 June 2021

Abstract: This work presents some numerical aspects of isogeometric boundary element methods (IGABEM). The behavior of hyper-singular and nearly-singular integration is first explored on the distorted NURBS surface. Several numerical treatments are proposed to enhance the quadrature in the framework of isogeometric analysis. Then a numerical implementation of IGABEM on the trimmed NURBS is detailed. Based on this idea, the surface crack problem is modeled incorporation with the phantom element method. The proposed method allows the crack to intersect with the boundary of the body while preserving the original parametrization of the NURBS-based CAD geometry.

Keywords: Isogeometric analysis; trimmed NURBS; singular integration; boundary element method; surface crack

The isogeometric analysis (IGA) uses the spline-basis functions, which are adopted to describe CAD geometry, to approximate the physical fields in analysis. The corporation of the IGA and boundary element methods (BEM) (abbreviated as IGABEM) provides a direct link between CAD and analysis since both methodologies are related to the boundary representation of the body. The IGABEM has gained more developments recently. The contributions regarding the numerical methods include the T-spline or subdivision surface based IGABEM [1,2], the extended IGABEM (XIGABEM) [3,4], the trimmed NURBS [5,6], the fast solution [7,8] and the Galerkin form [9,10], etc.

The study of IGA for trimmed NURBS geometry has been reported by [11–13]. Moreover, within the BEM community, Beer et al. [5] proposed a double mapping method to perform the integration on NURBS. This method relies on the establishment of the mapping area and is not able to be applied directly for a closed trimming curve, but later they proposed a stable extended B-splines scheme for IGABEM with trimmed geometry [14]. Wang et al. [6] developed an approach based on multi-patch non-singular BEM, and the Lagrange-multiplier was used to enforce the displacement continuity between the patches.

The numerical integration for trimmed elements is one of the key ingredients in IGA for trimmed CAD geometry. In the work by Kim et al. [11,12], the method which was applied in NURBS-Enhanced FEM [15] was adopted. This method requires the triangulation of the parametric domain of the trimmed element. Schmidt et al. [13] reconstructed the trimmed element by using new patch via interpolation or least square approximation. Beer et al. [5] created a mapping from an area bordered by two trimming curves and straight lines which connecting the ends of the trimming curves to the trimmed surfaces. The implementation is simple but excludes the case of holes cutout where closed trimming curves exist.

The IGABEM modeling of trimmed geometry provides an alternative approach to perform the surface crack simulation. In the boundary element method, cracks can be modeled by pairs of coinciding surfaces as external boundaries of the body. Then dual boundary integral equations are applied to form the linear system [16]. Fig. 1 illustrates a surface crack (or breaking crack) model. From (a) to (b), the two coinciding surfaces (crack) are inserted into a corner of a cube, thus breaking the boundary surfaces of the cube. If we take the front surface separately as in (c), the intersection curve of the crack surfaces and the front surface OA leads to a displacement discontinuity on the front surface. Hence the surface crack problem has two manifolds in numerical implementation: one is the coinciding crack surfaces inside the body domain, the other is the surface discontinuity along the geometry boundary. The latter can be considered as a problem of cracks in 2D plane or 3D shell conditions in the finite element method (FEM). In Lagrange-based elements (triangle or quadrilateral), one way to initiate and propagate the surface cracks is remeshing [17]. The extend FEM (XFEM) proposed by Belytschko and his team [18] allows the discontinuities modeled without changing the mesh discretization by introducing the enrichment functions such as the Heaviside function into the original basis functions.

Figure 1: Surface crack model by boundary element method. (a) The crack is modeled by two coincided surfaces as external boundary of the body geometry; (b) The crack is inserted to a corner of the cube, breaking the boundary surfaces of the cube; (c) If only the front surface is concerned, the intersection curve OA of the crack and front surface will create a discontinuity on the surface

When the problem comes into the Isogeometric analysis, analogous treatments can be done. The work contributing to the enriched IGA to model discontinuities can be referred to in, for example [19–21]. The original parametrization of the geometry is preserved by introducing the enrichment functions to describe the crack. Verhoosel et al. [22] proposed a reparametrization scheme via T-Splines to explicitly represent the crack. It should be noted that in their method, the original parametrization is lost due to the reparametrization. In order to form an analysis-suitable parametrization, the elements are distorted when the crack needs to take a turn. This will deteriorate the system condition number when the crack has a sharp turn.

In this paper, an implementation of IGABEM on trimmed NURBS surfaces is first outlined. The method presented in our work circumvents restrictions in [5] regarding the closed trimming curve and make improvements in terms of accuracy compared to previous work for the trimmed NURBS. Then a surface crack modeling technique is realized thanks to the development in trimmed NURBS.

Section 2 outlines the description and integration regarding trimmed NURBS, with details on collocation and prescribed boundary exertion. Section 3 introduces the way to model surface crack problem. The methods of singular and nearly singular integration are briefed in Section 4. Then follows the numerical examples on the study of (nearly) singular integration, the trimmed NURBS and the surface crack in Section 5. Section 6 concludes the work.

The geometry is generally created via the trimming operations on the NURBS surfaces. And the data of the NURBS surface and the relative trimming information is stored in .IGES file in CAD systems. Fig. 2a illustrates a trimmed surface

where Ri,j, Rk and Pi,j, Qk are the NURBS basis functions and control points of the surface and trimming curve respectively. Fig. 2b shows the corresponding parametric space of the trimmed surface and Cp(u) is the parametric trimming curve.

Although the trimming curve provided in .IGES file is in both physical and parametric form; there is no analytical relation between the parameter set (ξ, η) and u. This gives rise to challenges to numerical integration in analysis. For quadrature, the trimmed elements need to be determined for the first step, which requires a searching algorithm to find the parametric coordinates of the intersection points. The bisection searching algorithm by Schmidt et al. [13] is adopted in this work. After finishing searching the intersection points, the trimmed elements can be picked out and categorized as three types: triangle (‘3’), quadrangle (‘4’) and pentagon (‘5’) as in Fig. 2b. More unknown types of cutting can be transformed into the named types by knot insertion (mesh refinement). Note that the trimming orientation (Fig. 2b) will determine which side of the surface is maintained. In the present work, the left-hand side will remain.

Then the untrimmed elements will be judged to be cropped (‘−1’) or remained (‘1’) via the following procedure:

(1) Linearize the parametric trimming curve and find the shortest distance vector d from the central point of the element of interest to the line-segment represented trimming curve (Fig. 2b);

(2) Take the cross product of the distance vector d and trimming orientation vector t:

Figure 2: An example of trimmed surface, (a) in physical space; (b) in parametric space. The arrow in (b) denotes the direction of the trimming curve. The trimmed elements can be classified into three types: ‘3’ denotes a triange, ‘4’ denotes a quadrangle and ‘5’ denotes a pentagon. ‘1’ represents the untrimmed element and ‘−1’ is the cropped element

In this work, an approximated relation is built up between the (ξ, η) and u which will simplify the integration procedure. Inspired by the work of Beer et al. [5], the mapping method is applied locally for the trimmed elements. The parametric trimming curve is reconstructed from the same basis functions of the physical trimming curve. Then, the segmentation is done at every intersection point by knot insertion until C0 continuity is obtained. As presented in Fig. 3, the i-th segment of the parametric trimming curve Cip(u) is obtained by reducing to C0 at the intersection points ui0 and ui1. The next step would be to recover the parametric control points of this segment. The interpolation technique is adopted for this as following:

(1) Using Greville Abscissae to generate the sample points on the i-th physical trimming curve, i.e.,

where p is the order of the trimming curve and n is the number of the basis functions of the trimming curve;

(3) Using the point inversion algorithm [23] to retrieve the parametric coordinates Pks=ξks(ξks,ηks) of each xks;

The approximated relation between (ξ, η) and u is achieved by the reconstructed parametric trimming curve as ∀P(ξ,η)∈Cip(u),

where ξj and ηj are the components of the parametric control points Quj. The approximation property is illustrated in Fig. 4. In Fig. 4a, since the mapping from parametric space to physical space is linear, for curve order p = 2, only 3 sample points can exactly capture the trimmed geometry. In Fig. 4b, since the mapping from parametric space to physical space is nonlinear, the trimmed surface is bad approximated using 3 sample points, however, using more sample points by refining the trimming curve will improve the approximation as in Fig. 4c.

Now we are at the stage to establish the mapping from parent space to parametric space. For the case in Fig. 3, the Cip(u) would be η¯ = 1 (curve II), and the opposite edge η¯ = 0 (curve I) is a straight line. Taking u=ξ¯, using Eq. (7), we have

ξI=∑j=1nIRjI(ξ¯)ξjI,ηI=∑j=1nIRjI(ξ¯)ηjI,ξII=∑j=1nIIRjII(ξ¯)ξjII,ηII=∑j=1nIIRjII(ξ¯)ηjII. (8)

Figure 3: The mapping from parent space to parametric space. η¯ = 0 refers to curve I and η¯ = 1 is curve II

Figure 4: Approximation of the trimming surface, the green line is the physical trimming curve, the red dots are sample points, (a) the mapping from parametric to physical space is linear, the trimmed surface is exactly represented; (b) the mapping from parametric to physical space is nonlinear with 3 sample points on the trimming curve; (c) nonlinear with 6 sample points on the trimming curve

And the Jacobian transformation matrix from parent space to parametric space would be

J|ξ¯→ξ=[∂ξ∂ξ¯∂ξ∂η¯∂η∂ξ¯∂η∂η¯]=[(1−η¯)∂ξI∂ξ¯+η¯∂ξII∂ξ¯−ξI+ξII(1−η¯)∂ηI∂ξ¯+η¯∂ηII∂ξ¯−ηI+ηII], (10)

Note that for the pentagon-type trimmed elements, the parametric domain will be subdivided into two sub-quadrilaterals. Then the mapping scheme is set up for each sub-quadrilateral.

Compared to the work by Beer et al. [5], the proposed scheme for the trimmed NURBS can handle the closed trimming without further subdivision of the original patch. Integration for the trimmed NURBS is simplified and no triangulation on the parametric domain is needed as was done by Kim et al. [11]. This would facilitate the implementation of singular integration for the trimmed elements.

For a closed domain composed by trimless and compatible NURBS patches, the Greville abscissae (GA) is proved to be a nice way to generate the collocation points [24]. In IGABEM implementation, for those collocation points lie in sharp edges or corners, or when discontinuous basis functions are needed, these collocation points will be offset from the original place by (Fig. 6a)

ξs,i=ξs,i+α(ξs,i+1−ξs,i),orξs,i=ξs,i−α(ξs,i−ξs,i−1),0<α<1. (11)

The GA collocation can be used for trimmed NURBS as well. One only needs to be aware that some collocation points associated with the basis functions on the cropped part of a patch should be inactivated due to the trimming operation, providing a one to one correspondence is specified between the basis functions and collocation points (see [6] for details). Considering that in that case, some collocation points still locate outside the original patch, which will lead to an uncertain parametric position of the collocation point (it may locate on other patches with a different parametrization). We thus tried a mixed collocation scheme tailored for the trimmed NURBS patch based on the above modified GA as stated in the following step:

(1) Remove the collocation points located in the cropped elements and trimmed elements as denoted by blue dots in Fig. 5a;

(2) For each trimmed elements, (p + 1) (q + 1) collocation points are placed equally inside the elements (the yellow dots in Fig. 5b).

We note that in this way the number of collocation points is more than the number of basis functions. This leads to an over-determined system equations. Hence the boundary integral equations from collocation points in each trimmed element will be merged into the equations numbered by the global index of the basis functions of each trimmed element. For the pentagon-type trimmed elements, the (p + 1) (q + 1) collocation points are just simply placed in one of the sub-quadrilaterals (more details are outlined in Section 5.3.1) in this work. However, a more dedicated scheme can be realized to place the collocation points by comparing the area of the sub-quadrilaterals or uniformly distribute them into both sub-quadrilaterals. How to obtain a more efficient collocation scheme would be further explored in future work. Fig. 6 compares the GA collocation and mixed collocation methods for a square trimmed by a circle. Note that the points locating in the trimmed elements are moved to the center of the parent space to ensure a robust singular and nearly-singular integration.

Figure 5: Mixed collocation for the trimmed NURBS surface of order p = q = 2; (a) the collocation points generated by modified Greville abscissae, the blue ones are located in the cropped and trimmed elements and will be removed; (b) the final collocation points for trimmed NURBS patch by adding (p + 1) (q + 1) collocation points in each trimmed element

Figure 6: Two collocation methods for trimmed NURBS of p = q = 2. (a) The Greville Abscissae (GA) collocation (the points located in the trimmed elements are moved to the central of the parent space) (b) The mixed collocation

In this work, the prescribed Dirichlet or Neumann boundary conditions for each single trimless or trimmed NURBS patch is applied by the L2 projection. Suppose u¯ is the prescribed displacement field on a boundary patch S. The finite element approximation uh∈ Uh can be found by minimizing the L2 norm of the error between the prescribed and approximated displacement fields, i.e.,

We use (f, g)S to represent the inner product ∫s f· gdS. The minimization can be realized by letting the residual u¯−uh be orthogonal to any arbitrary vh ∈ Uh,

Substituting Eqs. (14) into (13), the discretized linear system can be obtained as

where the component of K is Kij = (Ri, Rj)S, and Fi=(Ri,u¯)S. Then the control coefficients d can be found by solving the linear system equations.

In this work, we outlined a simple approach named phantom element method (PEM, or in literature also called phantom node method [25]) to model the surface crack by IGABEM. In PEM, a crack will cut the element of interest into two parts. The degrees of freedom (DOFs) associated with this element will be duplicated; then each part has its independent DOFs to describe the primary physical fields (Fig. 7). The PEM has been investigated in finite element-based methods to represent strong and weak discontinuities. The PEM is also a way to model discontinuities without changing the mesh by integrating the split parts independently while no additional enrichment function is introduced. more details and applications can be referred in for example [25–28].

Fig. 7a illustrates the PEM in IGABEM to model the surface crack problem. The intersection curve (red curve in the figure) of the crack surfaces and the boundary surfaces can be considered as the trimming curve. By doubling the DOFs of the cut element, the surface discontinuity is described by performing the quadrature for each part in the same way as was done for trimmed elements as detailed in previous sections. For the cracked element, the displacement field is approximated as

u+=∑jNeRjdj,x∈Se+,u−=∑kNeRkdk,x∈Se−,Ne is the number of basis functions for the crack element. (16)

For the element containing the crack tip (the endpoint of the crack front), a local knot insertion (green line in the figure) is done to cut the element into two new elements. The element I is completely cut by the crack and the PEM is applied. However, we note that the NURBS does not allow local knot insertion due to the tensor-product nature of its basis functions. In order to do so, we first perform the knot insertion on each knot of the NURBS patch until its multiplicity is equal to p + 1 (p is the degree of each direction). Thus the continuity between the element reduces to C−1 such that discontinuous NURBS basis functions (or discontinuous rational Bezier basis) are obtained.

Figure 7: The phantom element method to model surface discontinuity. The red curve denotes the crack, which is regarded as a trimming curve to split the surface. The yellow rectangle represents the degrees of freedom (DOFs) associated with the element of interest. (a) For the completely cut element, the DOFs will be doubled directly. One group is associated with the upper part of the element; the additional group is with the lower part. The way of quadrature for trimmed NURBS can be applied for each part. (b) For the element containing the crack tip (an endpoint of the crack front), a knot (green line) is inserted to reduce the continuity such that the crack tip halted inside the element can be represented. In this way, the old element will be split into two new elements, and the obtained element I can then be applied with phantom element method

Although the continuity between the untrimmed elements is lost in order to halt the crack inside an element. The original parametrization is preserved and the extension to cracks with sharp turns or multi-cracks are straightforward. The future work will focus on performing the local knot insertion to reduce the continuity only for the crack tip element by using T-Splines in which the continuity in uncracked area can remain.

To overcome the degeneration in system matrices of BEM due to the overlapped crack surface, dual boundary integral equations are used. One can refer to the work in [29] for details.

Numerical integration for for IGABEM encounters is challengeable because the NURBS bases are rational, and the mesh can occasionally be highly distorted. This also complicates the singular and nearly-singular behaviors in the integrand of BEM. Efforts have been made to tackle this problem [30,31]. In this work, we investigate the conformal transformation technique [29] to the improve the (hyper-)singular integration in detail, then an L − 1/5 transformation is implemented to handle nearly singular integration.

The singularity subtraction technique is used for both untrimmed and trimmed elements. For the hyper-singular integral of the form

where H (s, x(ξ¯)) is the hyper-singular kernel, R(ξ¯) is the NURBS basis function and J¯(ξ¯) is the Jacobi transformation from parent space to physical space. By subtracting the Lorentz term and adding back semi-analytically, the integral becomes

I=∫02π∫0ρ^(θ)[F(ρ,θ)−F−2(θ)ρ2−F−1(θ)ρ]dρdθ+∫02πI−1(θ)ln⁡ρ^(θ)β(θ)dθ−∫02πI−2(θ)[γ(θ)β2(θ)+1ρ^(θ)]dθ=∫02π∫0ρ^(θ)F(ρ,θ)dρdθ−∫02π∫0ρ^(θ)F−2(θ)ρ2dρdθ−∫02π∫0ρ^(θ)F−1(θ)ρdρdθ+∫02π{I−1(θ)ln⁡ρ^(θ)β(θ)−I−2(θ)[γ(θ)β2(θ)+1ρ^(θ)]}dθ=I0−I2−I1+Iline (18)

The curvi-linear basis vectors at the source points mis=mi|ξ¯=ξ¯s,(i=1,2) and are calculated as

m1=[∂x∂ξ¯,∂y∂ξ¯,∂z∂ξ¯],m2=[∂x∂η¯,∂y∂η¯,∂z∂η¯]. (19)

such that the conformal transformation can be established. It can be concluded that λ reflects the local aspect ratio of the element at source point and cosψ indicates distortion of the element. Influence on singular integration by these two factors will be investigated in detail in the section of examples.

For nearly singular integration, two methods are realized. The first one is the recursive subdivision in the parametric domain and the other is a variable transformation technique. Note that many variable transformation techniques exist for nearly singular integration in Lagrange-based elements in literature (readers can refer to [32] for more details). However, as far as we know, these techniques are merely verified for NURBS elements. We adopt an L − 1/5 transformation, which was proposed by Hayami et al. [33]. The general procedure in Hayami’s work is as following:

(1) Find the closest point s′(ξ¯,η¯) to the source point s in the element of interest and the distance between them is d = ∥s − s′∥;

(2) Create a projection plane composed by sub-triangles Δ˜i=xis′xi+1. These triangles can be obtained by projecting the each vertex xi into the tangent plane containing s′ [33] or simply connecting the vertex and the closest point s˜′ which can be found by s′(ξ¯,η¯) [34];

(3) Map the parent space (ξ¯,η¯) to the sub-triangles Δ˜i, which is in Cartesian coordinate. The Jacobi transformation would be J˜;

(4) Introduce the polar coordinate transformation in Δ˜i (Fig. 8). And for the nearly singular integral of the type

In this work, we simply do the polar coordinate transformation in parent space and avoid the construction of Δ˜i. This would cause two problems as stated in [34]:

(1) In R(ρ)=(ρ+d)−15, d is in Cartesian space. Then the meaning of R is vague. This is overcome by regularizing d with the characteristic element size h, i.e., d0 = d/h;

(2) The parent space is insensitive to the element distortion. This is improved by using the conformal and Sigmoidal transformation which are used for singular integration in this work.

is checked for various geometries in this section. The related reference solutions are obtained by Mathematica©. For singular integration, three cases are involved: conformal transformation (‘con’), conformal transformation and Sigmoidal transformation (‘con + sig’), and the original SST (‘ori’). For nearly-singular integration, two cases, the adaptive subdivision and the L − 1/5 transformation are compared. Gauß-Legendre rule is used for quadrature. For the integration performed in sub-triangles (the singular integration and the nearly-singular integration by L − 1/5 transformation), ‘ngp_s’ denotes the number of Gauß points in the angular direction in each sub-triangle and ‘ngp_t’ represents the number of Gauß points in the radial direction.

We first verify the hyper-singular integration over a quarter of a disc given by Coons parametrization as in Fig. 9. For the source point ξs moving from (0.5, 0.5) to (1, 1), ψ varies from 90° to 180° while λ keeps to be 1. Thus the integral undergoes an increasing near-singularity. We take ξs as (0.7, 0.7), (0.9, 0.9) and (0.99, 0.99), and ψ is equal to 123.2°, 157.3° and 177.6° respectively. Since one element is used, only the singular integration is involved here. It can be seen from Fig. 10 that all the three cases will finally converge to a stable precision with increasing ngp_s. For ξs (0.7, 0.7), ‘con’ and ‘con + sig’ show a very close convergence rate and reaches O(10−11) at ngp_s = 14, while the original SST uses 20 Gauß points. For ξs (0.9, 0.9), the discrepancy of convergence rate among the three cases becomes significant. And the ‘con + sig’ performs the best with ngp_s = 12 to reach O(10−8), which costs 28 Gauß points in the angular direction for original SST. The use of only conformal transformation presents an intermedium convergence rate in this case. However, the convergence error increases with the source points ξs approaches to (1, 1) for a fixed ngp_t = 10, especially for ξs (0.99, 0.99) although the improved methods show much faster convergence rates than the original method, the error however is only at O(10−2). Fig. 11 illustrates the study with increasing the Gauß points in radial direction ngp_t, while fixing the number of Gauß point in the angular direction ngp_s. We fixed ngp_s = 18 for the improved method and ngp_s = 36 for original method since we note that the original method needs more Gauß points in angular direction to get converged. It can be observed that the improved method shows the same convergence trend as the original method with regard to the number of Gauß Points in the radial direction. And both methods achieve O(10−10) for ξs at (0.7, 0.7) and (0.9, 0.9) when ngp_t is 14. However, further increasing the number of Gauß points in the radial direction will accumulate the integration error. For ξs at (0.99, 0.99), the original method converges to O(10−2) with increasing the ngp_t while the improved method reaches O(10−7). This is due to the fact that even 36 Gauß points in the angular direction is insufficient to circumvent the near-singularity in the integrand for the original method.

Figure 9: A quarter of a disc by Coons parametrization. (a) The parameter line; (b) the single element for singular integration with source points P1: ξs (0.7, 0.7), P2: ξs (0.9, 0.9) and P3: ξs (0.99, 0.99); (c) the parametric space

Figure 10: Convergence check with respect to ngp_s for hyper-singular integral over a quarter of disc by Coons parametrization (a) ξs(0.7, 0.7) (b) ξs(0.9, 0.9) (c) ξs(0.99, 0.99)

Figure 11: Convergence study with respect to ngp_t for hyper-singular integral over a quarter of disc by Coons parametrization

We still perform the hyper-singular integral over a quarter of disc, but with the parametrization degenerated at the pole as in Fig. 12 such that the distortion angle ψ is always 90° regardless of ξs. For ξs moving from (0.5, 0.5) to (0.5, 0), the local aspect ratio λ changes from 1 to + ∞ (assume λ ≥ 1). We take ξs at (0.5, 0.1), (0.5, 0.01), (0.5, 0.001) and (0.001, 0.001), and λ is 6.0, 60.3, 603.6 and 706.7, respectively. The convergence results with respect to the number of Gauß points in the angular direction ngp_s while keeping ngp_t = 10 and the errors of the integral are compared in Fig. 13. We can see that for all cases, the ‘con + sig’ scheme outperforms with the fastest convergence rate and the smallest error (for the former three positions of ξs, the precision at O(10−12) with about 20 Gauß points in the angular direction). The conformation transformation only shows a better convergence than the original method, but both methods fail to get converged at a stable precision with increasing λ. For ξs (0.001, 0.001), the integral by the original method results in arbitrary value thus its error curve is not plotted in Fig. 13d. The ‘con + sig’ scheme converges at O(10−6) with 14 Gauß points in the angular direction and the conformal transformation only presents a slow convergence and higher error than ‘con + sig’ finally although the error is much lower than ‘con + sig’ when ngp_s < 12.

Figure 12: A quarter of a disc by parametrization degenerated at the pole. (a) The parameter line; (b) the single element for singular integration with source points P1: ξs (0.5, 0.1), P2: ξs (0.5, 0.01), P3: ξs (0.5, 0.001) and P4: ξs (0.001, 0.001); (c) the parametric space

Figure 13: Convergence study with respect to ngp_s for hyper-singular integral over a quarter of disc by parametrization degenerated at the pole (a) ξs(0.5, 0.1) (b) ξs(0.5, 0.01) (c) ξs(0.5, 0.001) (d) ξs(0.001, 0.001)

We also studied the error convergence trend concerning the number of Gauß points in the radial direction ngp_t with a fixed ngp_s = 18 for an improved method and ngp_t = 36 for the original method. We can conclude from Fig. 14 that the integral shows a stable behavior for these cases, while with too many Gauß points in the radial direction, the error starts to accumulate slowly. The accuracy is poor by the original method for ξs at (0.5, 0.01) and (0.5, 0.001), due to the fact that ngp_s = 36 is far from sufficient to overcome the near-singularity.

First of all, we give some explanations on ‘complex distortion’: (1) the distortion angle ψ gets close to π and the local aspect ratio λ deviated from 1, or the λ cannot reflect the distortion of the geometry; (2) The parameter line shows an obvious change in its direction for an element. For example, a quarter of an ellipse with parametrization degenerated at the pole has a close-to-unit λ at its parametric centre point regardless of the variation of the ratio of semi-major and semi-minor axes (a/b) as illustrated in Fig. 16. We use a single element to discretize the quarter of an ellipse and place ξs at (0.5, 0.5) (Fig. 15), then perform the hyper-singular integral over the domain with a/b = 2, 10 and 20. The results are plotted in Fig. 17. It can be observed that both the improved method and original method present close convergence trends with respect to ngp_s. Moreover, they can achieve low error at O(10−12) for small a/b, but the integral becomes diverged when a/b is larger. Tables 1 and 2 list the value of the components of the integral. It can been observed that for a/b = 2, all terms achieves an error below O(10−7) within 22 Gauß points. While for a/b = 20, the terms I−1, I−2 and Iline easily reaches the O(10−7) within 22 Gauß points, nevertheless I0 stays in poor precision. It can be referred that for complex distortion of the geometry, it is the bad quadrature in I0 that leads to the poor accuracy in the final integration.

Figure 14: Convergence study with respect to ngp_t for the hyper-singular integral over a quarter of a disc by parametrization degenerated at the pole

Figure 15: A quarter of an ellipse by parametrization degenerated at the pole. (a) The parameter line; (b) the single element for singular integration with source points P: ξs(0.5, 0.5); (c) the parametric space

Figure 16: A quarter of an ellipse by parametrization degenerated at the pole with source point P: ξs (0.5, 0.5). (a) a/b = 2, λ = 1.2 and ψ = 53.1° for P; (b) a/b = 10, λ = 1.2 and ψ = 11.4° for P; (c) a/b = 20, λ = 1.2 and ψ = 5.7° for P

Figure 17: Convergence check with respect to ngp_s for hyper-singular integral over a quarter of ellipse with varied a/b

Figure 18: Convergence check with respect to ngp_s for hyper-singular integral over a quarter of ellipse with a/b = 20

A remedy for this deterioration would be to per form certain mesh refinement as in Fig. 19c with ω = 0.01 (see Fig. 19c for the definition), i.e., knots insertion in NURBS. Then multiple elements are used to discretize the domain, and near-singular integration is introduced. We use the adaptive subdivision scheme for nearly-singular integration. The convergence is checked for the case of a/b = 20 with respect to the ngp_s in the singular element, which is shown in Fig. 18. It is seen that with certain mesh refinement, the improved method achieves a fast convergence as previous case studies and reaches O(10−8) with about 20 Gauß point in the angular direction, in contrast, the original method stays in poor accuracy and slow convergence rate. This example indicates a quadrature scheme for the geometry with complex distortion, which we name as the ‘transformation+subdivision’ scheme. We note that the final result of singular integration by this method is determined by both singular and nearly-singular integration, and it will be studied further in the next section.

Figure 19: The hyper-singular integral over a quarter of an ellipse with different mesh design. (a) Mesh A with knot vectors ξ: [0,0,0,0.5−ω,0.5+ω1,1,1],η: [0,0,0,1,1,1]; (b) mesh B with knot vectors ξ: [0,0,0,1,1,1],η: [0,0,0,0.5−ω,0.5+ω1,1,1]; (c) mesh C with knot vectors ξ: [0,0,0,0.5−ω,0.5+ω1,1,1],η: [0,0,0,0.5−ω,0.5+ω1,1,1]

In this section, the nearly-singular integration by the proposed L − 1/5 transformation is studied in detail. We take a quarter of an ellipse with a/b = 2 as an example. Three mesh refinement configurations are set up as in Fig. 19 such that the near-singularity arises in the neighbor of the singular element by introducing a relative distance factor ω. Note that different spatial position of the source point and nearly-singular element can be checked simultaneously. For the singular element, ngp_s = 18 and ngp_t = 10 will be fixed. And for the nearly-singular element, we increase the ngp_s and ngp_t from 4 to 80 respectively. The integration error for all the cases are plotted in Fig. 20. For all the mesh configurations, the L − 1/5 transformation can achieve O(10−6) with about 20 Gauß points in both the radial and the angular directions for each sub-triangle. Moreover, the convergence trend is similar in both directions. The accuracy can be maintained with ω = 0.001. While further reduce the value of ω, the error will increase. How to improve the precision for nearly-singular integration will be further explored in future work.

Table 3 presents the reduction in total number of Gauß points if compared with the adaptive subdivision scheme. It can referred that the total number of Gauß points is reduced by two orders of magnitude when ω ≤ 0.01 for all the mesh configurations.

Figure 20: Convergence study for nearly-singular integration by L − 1/5 transformation (a) mesh A with ω=0.1 (b) mesh B with ω=0.1 (c) mesh C with ω=0.1 (d) mesh A with ω=0.01 (e) mesh B with ω=0.01 (f) mesh C with ω=0.01 (g) mesh A with ω=0.001 (h) mesh B with ω=0.001 (i) mesh C with ω=0.001

Two case studies are performed to evaluate the proposed methods for trimmed NURBS. The relative error in displacement or traction L2 norm over the boundary of the domain is used to measure the results. They are given as

∥eu∥L2=∫S(u−uext)(u−uext)TdS∫SuuextTdS,∥et∥L2=∫S(t−text)(t−text)TdS∫SttextTdS. (25)

The default order of the basis functions is 2. In the degree elevation process, the higher-order will be 3.

A cube with edge L = 1 cut a cylindrical surface with radius r at its central is studied in this section as in Fig. 21. The faces x = 0, y = 0 and z = 0 are applied with normal displacement constraints and the top face z = 1 is applied by a uniform traction in z direction. Remaining degrees of freedom are traction-free. The material constants E = 1000 and ν = 0.3 such that the analytical displacement field over the domain would be

Since the stress field is constant in this example and the geometry is exactly represented, the source of numerical error comes from the integration.

Figure 21: The central of a cube trimmed by a cylindrical surface (red). The unspecified degrees of freedom are zero tractions

Figure 22: Collocation scheme for pentagon-type trimmed element in top face of the cube. The green line is the untrimmed element boundary; the blue line denotes the parameter line and red dots are collocation points (a) Coarse mesh with Scheme A (b) Coarse mesh with Scheme B (c) Fine mesh with Scheme A (d) Fine mesh with Scheme B

The aforementioned collocation schemes are first investigated on a quarter of the trimmed cube with r = 0.15 of cylindrical surface to check the singular integration for the pentagon-type elements in the trimmed NURBS geometry. Fig. 22 illustrates collocation Scheme A (the collocation points are in the sub-quadrilateral with all straight edges) and Scheme B (the collocation points in the sub-quadrilateral with the curved edge) for a coarse mesh (4 elements) and a fine mesh (9 elements). When the collocation points are placed in one sub-quadrilateral, the integration over the other sub-quadrilateral will be performed as a regular one (the nearly singular quadrature is used). It can be seen that for the collocation points in Scheme A, good parametrization is obtained for singular integration. While in Scheme B, the collocation points are located in the sub-quadrilateral with a highly distorted parametrization, which will increase the difficulty for singular integration. Two singular integration methods are compared. One is performing the singular integration with the conformal and Sigmoidal transformations directly (denoted by ‘trans’); The other will subdivide the parametric space as was done in Fig. 19, and then the relative transformations will be adopted for singular integration, and the remaining sub-domains will be treated as nearly singular integration (i.e., the ‘transformations + subdivision’ or shorted as ‘trans + subdi’). The coefficient ω = 0.01 is taken for the subdivision. Table 4 presents the results for the test. It can be observed that Scheme A achieves orders of magnitude higher precision than Scheme B in singular integration scheme ‘trans’ due to the complex distortion of the parametrization. Analogous to the example in Section 5.1.3, the ‘trans + subdi’ method improves the results when the collocation points are in the distorted sub-quadrilateral. It is seen that for mesh (b), the ‘trans + subdi’ method gains the same order (O(10−7)) for both collocation Schemes A and B, while for mesh (c), the ‘trans + subdi’ only reaches O(10−4) although two orders of magnitude higher than ‘trans’. Reducing ω is supposed to further improve the precision. However, due to the restriction in nearly singular integration, the final result would be difficult to improve. This requires a further local refinement closed to the trimming curve in the untrimmed mesh to reduce the distortion in trimmed elements.

Figure 23: Meshes of the trimmed cube by uniform refinement. (a) 2 × 2 × 2 elements; (b) 4 × 4 × 4 elements; (c) 6 × 6 × 6 elements; (d) 8 × 8 × 8 elements

Then the whole model is used, and four mesh configurations are obtained by uniform knot insertion for the untrimmed cube (Fig. 23). The collocation Scheme A for the pentagon-type trimmed elements in the mixed collocation method and ‘trans + subdi’ singular quadrature method for the trimmed elements are used for all four meshes. Table 5 presents the relative error in the displacement L2 norm for both mixed and GA collocation. It can be observed that the smallest error can reach O(10−8) in both collocation methods. The GA collocation maintains the error below O(10−5). The mixed collocation method obtains better accuracy however, for mesh 3 in Fig. 23, the error is of O(10−4). This is because the collocation points in the small triangle-type trimmed elements will lead to nearly-singular integration. We note that both cases give a smaller error than the result in [6] where the error stays at O(10−3) although it was given in a kind of discretized L2 norm of some selected points.

In this section, the convergence behavior of the IGABEM for the trimmed NURBS is studied by applying the Kelvin fundamental solution to the trimmed geometry. The cube cut by the cylindrical surface in the patch test is used. And a unit point force in z direction is applied at sP (0, 0, 1.5), then corresponding displacement and traction fields in the closed domain Ω¯ invoked by the point force are

ui(x)=116πμ(1−ν)r[(3−4ν)δi3+r,ir,3],ti(x)=−18π(1−ν)r2{∂r∂n[(1−2ν)δi3+3r,ir,3]−(1−2ν)(r,in3−r,3ni)},∀x∈Ω¯, (27)

where r = |sP − x| and n is the unit out normal. The material constant μ=E2(1+ν). We take μ = 1 and ν = 0.3 for the case study.

Fig. 24 compares the two collocation methods for the pure Dirichlet boundary condition (BC). Since the displacement fields for all the faces are prescribed by the L2 projection at the beginning of the analysis, the error can be considered as the approximation error, and it is the same for the two collocation methods. It can be observed that both methods can obtain a convergence result in the solution of traction fields. Moreover, the result of the GA collocation is more accurate than the mixed collocation. Based on this, we select the GA collocation for further studies.

Fig. 25 shows the convergence curves for a degree elevation of the basis functions from p = 2 to 3 for the pure Dirichlet BC. The optimal order of convergence rates (oc) of the prescribed displacement fields for both p = 2 (oc = 3.04) and p = 3 (oc = 4.26) are obtained. For the traction fields, the order elevation will reduce the error. However, we note that for both cases, the convergence results are sub-optimal. And the oc for p = 3 (1.95) is worse than p = 2 (2.46). One possible reason could be that the integration error is not small enough compared to the approximation error to give the optimal convergence. For this example, the integration error is known at O(10−5) from the patch test. And the approximation error of the problem can be known from the prescribed displacement fields. It can be seen that after degree elevation, the approximation error approaches closely to O(10−5), which is almost at the same level as the integration error. This may explain the deterioration in oc for p = 3 compared to p = 2. More theoretical and numerical problem on the integration error, approximation error and their influences on the convergence rate is still an open topic. It will be further investigated in the future work.

Figure 24: Comparison of the two collocation methods for a pure Dirichlet problem in the convergence study

Figure 25: Relative error in displacement and traction fields for pure dirichlet boundary condition

The mixed BCs are checked where the top z = 1 and bottom z = 0 faces are applied with Dirichlet BC, and the remaining faces are applied with Neumann BC. Fig. 26 compares the convergence results for the displacement and traction fields in pure Dirichlet BC and mixed BC. It can be observed that the convergence trend can be obtained for mixed BC as well.

Figure 26: Relative error in displacement and traction fields for pure dirichlet boundary condition and mixed boundary condition

In this section, the surface crack modeled by the phantom element method is checked by the edge crack example. Fig. 27a illustrates the numerical model of the edge crack, where the bottom face is fixed, and the top face is applied with uniform traction in the normal direction. The material constants E = 1000 and v = 0.3. Fig. 27b gives the deformation result of the edge crack.

This example present the possibility to model the surface crack problem using IGABEM while the original parametrization can be preserved. Further verification in some fracture parameters such as stress intensity factors as well as the surface crack propagation for non-trivial industrial parts will be done in the future work.

In this work, several numerical aspects, such as the singular integration, trimmed NURBS and surface crack modeling, of the IGABEM are investigated in detail. The singular and nearly singular integration has been studied for NURBS elements with high element aspect ratio or elemental distortion, which are revealed to be important in the accuracy in trimmed NURBS implementation. The conclusions are

(1) the singular integration exhibits a sensitive behavior on the element shape (or parametrization quality), and the distorted element usually with bad parametrization shows a deterioration on the accuracy of the singular integral. The proposed ‘transformation + subdivision’ scheme is shown to be a remedy for this issue of IGABEM;

(2) the proposed IGABEM for trimmed NURBS can handle the closed trimming curve and multi-curves in a single patch;

(3) the proposed surface crack modeling allows the crack to split the boundary of the body geometry while the original parametrization is preserved.

Future work will focus on describing crack tip which can stop inside an element without inserting knot at the tip point on the body surface. Adaptive spline technique is also of interest since it is crucial to obtain an accurate stress field around the crack front.

Funding Statement: Haojie Lian thanks the financial support of National Natural Science Foundation of China (NSFC) under Grant (No. 51904202).

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Scott, M. A., Simpson, R. N., Evans, J. A., Lipton, S., Bordas, S. P. A. et al. (2013). Isogeometric boundary element analysis using unstructured T-splines. Computer Methods in Applied Mechanics and Engineering, 254, 197–221. DOI 10.1016/j.cma.2012.11.001. [Google Scholar] [CrossRef]

2. Chen, L., Lu, C., Lian, H., Liu, Z., Zhao, W. et al. (2020). Acoustic topology optimization of sound absorbing materials directly from subdivision surfaces with isogeometric boundary element methods. Computer Methods in Applied Mechanics and Engineering, 362(3–5), 112806. DOI 10.1016/j.cma.2019.112806. [Google Scholar] [CrossRef]

3. Peake, M. J., Trevelyan, J., Coates, G. (2013). Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems. Computer Methods in Applied Mechanics and Engineering, 259, 93–102. DOI 10.1016/j.cma.2013.03.016. [Google Scholar] [CrossRef]

4. Peake, M., Trevelyan, J., Coates, G. (2015). Extended isogeometric boundary element method (XIBEM) for three-dimensional medium-wave acoustic scattering problems. Computer Methods in Applied Mechanics and Engineering, 284, 762–780. DOI 10.1016/j.cma.2014.10.039. [Google Scholar] [CrossRef]

5. Beer, G., Marussig, B., Zechner, J. (2015). A simple approach to the numerical simulation with trimmed CAD surfaces. Computer Methods in Applied Mechanics and Engineering, 285(37–40), 776–790. DOI 10.1016/j.cma.2014.12.010. [Google Scholar] [CrossRef]

6. Wang, Y., Benson, D. J., Nagy, A. P. (2015). A multi-patch nonsingular isogeometric boundary element method using trimmed elements. Computational Mechanics, 56(1), 173–191. DOI 10.1007/s00466-015-1165-y. [Google Scholar] [CrossRef]

7. Marussig, B., Zechner, J., Beer, G., Fries, T. P. (2015). Fast isogeometric boundary element method based on independent field approximation. Computer Methods in Applied Mechanics and Engineering, 284(39–41), 458–488. DOI 10.1016/j.cma.2014.09.035. [Google Scholar] [CrossRef]

8. Sun, F. L., Dong, C. Y., Wu, Y. H., Gong, Y. P. (2019). Fast direct isogeometric boundary element method for 3D potential problems based on HODLR matrix. Applied Mathematics and Computation, 359(39–41), 17–33. DOI 10.1016/j.amc.2019.04.030. [Google Scholar] [CrossRef]

9. Feischl, M., Gantner, G., Praetorius, D. (2015). Reliable and efficient a posteriori error estimation for adaptive IGA boundary element methods for weakly-singular integral equations. Computer Methods in Applied Mechanics and Engineering, 290(39–41), 362–386. DOI 10.1016/j.cma.2015.03.013. [Google Scholar] [CrossRef]

10. Aimi, A., Diligenti, M., Sampoli, M. L., Sestini, A. (2016). Isogemetric analysis and symmetric Galerkin BEM: A 2D numerical study. Applied Mathematics and Computation, 272(11), 173–186. DOI 10.1016/j.amc.2015.08.097. [Google Scholar] [CrossRef]

11. Kim, H. J., Seo, Y. D., Youn, S. K. (2009). Isogeometric analysis for trimmed CAD surfaces. Computer Methods in Applied Mechanics and Engineering, 198(37–40), 2982–2995. DOI 10.1016/j.cma.2009.05.004. [Google Scholar] [CrossRef]

12. Kim, H. J., Seo, Y. D., Youn, S. K. (2010). Isogeometric analysis with trimming technique for problems of arbitrary complex topology. Computer Methods in Applied Mechanics and Engineering, 199(45–48), 2796–2812. DOI 10.1016/j.cma.2010.04.015. [Google Scholar] [CrossRef]

13. Schmidt, R., Wüchner, R., Bletzinger, K. U. (2012). Isogeometric analysis of trimmed NURBS geometries. Computer Methods in Applied Mechanics and Engineering, 241–244(5), 93–111. DOI 10.1016/j.cma.2012.05.021. [Google Scholar] [CrossRef]

14. Marussig, B., Zechner, J., Beer, G., Fries, T. P. (2017). Stable isogeometric analysis of trimmed geometries. Computer Methods in Applied Mechanics and Engineering, 316(2), 497–521. DOI 10.1016/j.cma.2016.07.040. [Google Scholar] [CrossRef]

15. Sevilla, R., Fernández-Méndez, S., Huerta, A. (2008). NURBS-enhanced finite element method (NEFEM). International Journal for Numerical Methods in Engineering, 76(1), 56–83. DOI 10.1002/nme.2311. [Google Scholar] [CrossRef]

16. Mi, Y., Aliabadi, M. H. (1992). Dual boundary element method for three-dimensional fracture mechanics analysis. Engineering Analysis with Boundary Elements, 10(2), 161–171. DOI 10.1016/0955-7997(92)90047-B. [Google Scholar] [CrossRef]

17. Cisilino, A. P., Aliabadi, M. H. (2004). Dual boundary element assessment of three-dimensional fatigue crack growth. Engineering Analysis with Boundary Elements, 28(9), 1157–1173. DOI 10.1016/j.enganabound.2004.01.005. [Google Scholar] [CrossRef]

18. Moës, N., Dolbow, J., Belytschko, T. (1999). A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1), 131–150. DOI 10.1002/(ISSN)1097-0207. [Google Scholar] [CrossRef]

19. De Luycker, E., Benson, D. J., Belytschko, T., Bazilevs, Y., Hsu, M. C. (2011). X-FEM in isogeometric analysis for linear fracture mechanics. International Journal for Numerical Methods in Engineering, 87(6), 541–565. DOI 10.1002/nme.3121. [Google Scholar] [CrossRef]

20. Ghorashi, S., Valizadeh, N., Mohammadi, S. (2012). Extended isogeometric analysis for simulation of stationary and propagating cracks. International Journal for Numerical Methods in Engineering, 89(9), 1069–1101. DOI 10.1002/nme.3277. [Google Scholar] [CrossRef]

21. Nguyen, V. P., Anitescu, C., Bordas, S. P. A., Rabczuk, T. (2015). Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 117, 89–116. DOI 10.1016/j.matcom.2015.05.008. [Google Scholar] [CrossRef]

22. Verhoosel, C. V., Scott, M. A., de Borst, R., Hughes, T. J. R. (2011). An isogeometric approach to cohesive zone modeling. International Journal for Numerical Methods in Engineering, 87(1–5), 336–360. DOI 10.1002/nme.3061. [Google Scholar] [CrossRef]

23. Piegl, L., Tiller, W. (1995). The NURBS book. New York, USA: Springer. [Google Scholar]

24. Auricchio, F., Veiga, L. B. D., Hughes, T. J. R., Reali, A., Sangalli, G. (2010). Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 20(11), 2075–2107. DOI 10.1142/S0218202510004878. [Google Scholar] [CrossRef]

25. Rabczuk, T., Zi, G., Gerstenberger, A., Wall, W. A. (2008). A new crack tip element for the phantom node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 75(5), 577–599. DOI 10.1002/nme.2273. [Google Scholar] [CrossRef]

26. Hansbo, A., Hansbo, P. (2004). A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 193(33–35), 3523–3540. DOI 10.1016/j.cma.2003.12.041. [Google Scholar] [CrossRef]

27. Mergheim, J., Kuhl, E., Steinmann, P. (2005). A finite element method for the computational modelling of cohesive cracks. International Journal for Numerical Methods in Engineering, 63(2), 276–289. DOI 10.1002/(ISSN)1097-0207. [Google Scholar] [CrossRef]

28. Chau-Dinh, T., Zi, G., Lee, P. S., Rabczuk, T., Song, J. H. (2012). Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 92–93(11), 242–256. DOI 10.1016/j.compstruc.2011.10.021. [Google Scholar] [CrossRef]

29. Peng, X., Atroshchenko, E., Kerfriden, P., Bordas, S. (2017). Isogeometric boundary element methods for three dimensional static fracture and fatigue crack growth. Computer Methods in Applied Mechanics and Engineering, 316(23–24), 151–185. DOI 10.1016/j.cma.2016.05.038. [Google Scholar] [CrossRef]

30. Gong, Y., Dong, C., Qin, F., Hattori, G., Trevelyan, J. (2020). Hybrid nearly singular integration for three-dimensional isogeometric boundary element analysis of coatings and other thin structures. Computer Methods in Applied Mechanics and Engineering, 367(5566), 113099. DOI 10.1016/j.cma.2020.113099. [Google Scholar] [CrossRef]

31. Han, Z., Huang, Y., Cheng, C., Liang, Y., Hu, Z. et al. (2020). The semianalytical analysis of nearly singular integrals in 2D potential problem by isogeometric boundary element method. International Journal for Numerical Methods in Engineering, 121(16), 3560–3583. DOI 10.1002/nme.6370. [Google Scholar] [CrossRef]

32. Ye, W. (2008). A new transformation technique for evaluating nearly singular integrals. Computational Mechanics, 42(3), 457–466. DOI 10.1007/s00466-008-0262-6. [Google Scholar] [CrossRef]

33. Hayami, K., Matsumoto, H. (1994). A numerical quadrature for nearly singular boundary element integrals. Engineering Analysis with Boundary Elements, 13(2), 143–154. DOI 10.1016/0955-7997(94)90017-5. [Google Scholar] [CrossRef]

34. Hayami, K. (2005). Variable transformations for nearly singular integrals in the boundary element method. Publications of the Research Institute for Mathematical Sciences, 41(4), 821–842. DOI 10.2977/prims/1145474596. [Google Scholar] [CrossRef]