Numerical investigations in non-watertight models based on a surface-independent discretization boundary element method

Abstract

It is well known that boundary integral equations are exact mathematical representations of the governing differential equations of a boundary value problem when the integrals are written over a closed-shape boundary representation (B-representation) of the domain, usually reffered to as a watertight B-representation. However, practical geometric design technics (namely, NURBS surfaces) often do not render a watertight B-representation. Non-watertight geometric models with small gaps and overlaps are often generated in the design stage of projects. Based on a proposed surface-independent discretization approach, the present study investigates how unsought gaps affect the response of boundary element models of linear elasticity problems. The developed surface-independent discretization is applied to discretize multiple-patches NURBS B-representation geometries. Linear triangular and quadrilateral elements are adopted to discretize the independent surfaces. Generalized discontinuous elements at the edges of the visible areas of the NURBS parametric spaces are detected by a Level Set function. An offset collocation strategy is adopted for the nodes at the edges of the visible part of the parametric spaces. Thus, singularities and near singularities due to collocation are avoided in the BEM equations. The influence of gaps in the convergence of the L2-norm of boundary displacement error is verified in a 3D example with an available analytical solution. A second example with available numerical solution is analyzed with a non-watertight BEM discretization for qualitative boundary field validation. Finally, a non-watertight B-representation geometry of a crane hook is analyzed. The obtained results have pointed out that, as long as the gaps (and overlaps) are small enough, BEM models built up from non-watertight geometries may produce valuable solutions for practical purposes.

Keywords:
Non-watertight B-representations; NURBS surfaces discretization; boundary elements; convergence tests

Graphical Abstract

1 INTRODUCTION

Linear structural analyses of mechanical components used for various purposes are vital to ensure the structural safety in the design phase. To improve safety, engineers spend a significant amount of time designing, modeling, and analyzing each component. The geometric models of the mechanical components are first developed by designers in computer-aided design (CAD) software. Then, numerical methods, such as the finite or boundary element methods (FEM or BEM) (Zienkiewicz and Taylor, 1967; Brebbia et al., 1984), are employed to build discrete analysis models from suitable CAD models. Solvers are then employed to perform the numerical analysis, and the results are visualized in a pos-processing stage. A suitable geometric model must be provided from CAD software to be able to fast generate the discrete analysis models (Cottrell et al., 2009, Bazilevs et al., 2010). The most common solid geometric models are built as a set of independent parametric surfaces, which describes the boundary of the domain. Such models are known in the literature as boundary-representation (B-representation) models (Sederberg and Parry, 1986). Usually, B-representation models must be converted into analysis-suitable models. This involves several intermediate steps (Coll et al., 2014). The first step is to obtain a B-representation without gaps and overlaps in the intersections of the independent surfaces. Then, a closed-shape surface mesh with conforming surface-to-surface discretization must be obtained. For domain-based methods, i.e., finite elements, the last step required is to generate a standard volumetric mesh for finite element analysis purposes.

Problems often arise in ensuring no gaps or overlaps in the surface intersections. With the standard CAD technologies, namely Non-Uniform Rational B-splines (NURBS) surfaces (Piegl and Tiller, 1995; Roger, 2001), small gaps are generally inevitable (Sederberg et al., 2003, Bazilevs et al., 2010). If small gaps exist in a B-representation model, it is called a non-watertight geometric model. Non-watertight models are usually avoided in the design-through-analysis process, since for the standard analysis technology, namely the finite element method (Zienkiewicz and Taylor, 1967) and other domain and meshless methods (Belytschko et al., 1994; Strouboulis et al., 2001), the existence of gaps prevents the standard volumetric mesh (or regularly distributed internal points) generation. This issue is by far a nontrivial task for complex geometric engineering designs, as Cottrell et al. (2009) pointed out. The conversion of geometric models to analysis suitable models have been studied by Ted Blacker, Manager of Simulation Sciences, Sandia National Laboratories, and usually takes around 80% of the whole design-through-analysis process time (Hugues et al., 2005). Linear differential equations governing a given boundary value problem can be rewritten as boundary integral equations when the fundamental solution is known (Hsiao and Wendland, 2008). The boundary integral equations of linear elasticity are well-established nowadays and can be properly solved by the Boundary Element Method (BEM) (Aliabadi, 2002). The main characteristic of BEM is that domain meshes are avoided and the discretization is required only at the boundaries of the solids. This feature makes the construction of discrete BEM models from CAD geometries largely attractive because both CAD and BEM are based on the B-representation of the solids. Based on the so-called isogeometric concept (Cottrell et al., 2009), the construction of BEM models directly from CAD geometries for numerical analysis purposes has been reported in several works in the literature (Simpson et al., 2012, Simpson et al., 2013, Lian, et al., 2013; Scott et al., 2013, Marussig et al. 2015, Peng et al. 2017, Oliveira et al. 2020). The boundary integral equations are exact mathematical representations of the problem if they are written over closed-shape B-representations (watertight B-representations). If only a non-watertight B-representation is available, the boundary integral equations are only approximately true in the best case. Since integrals can be understood as infinity summations, the contributions that will not appear or appear spuriously in the integral equations due to the gaps or overlaps will introduce errors in the previously mathematically exact boundary integral equations. However, as reported by Daumas et al. (2022) for 2D analysis, if the gaps (and overlaps) in the non-watertight B-representations are small enough, such errors will be small, and the boundary integral equations written over the non-watertight B-representations can still be solved by the BEM.

In the present work, a surface-independent discretization scheme is developed for the discrete BEM model generation from NURBS-based B-representation geometries. Generalized discontinuous elements are adopted at the edges of the NURBS parametric spaces to allow the study of discrete models with very small gaps between surfaces without resulting singularities due to close-collocation points. The influence of the gaps in the numerical results of linear elasticity problems are verified in three examples reported.

2 NURBS BASED B-REPRESENTATION GEOMETRIC MODELS

The B-representation geometric models assume that the boundary $\mathbf{\Gamma }$ of a given domain can be described by the union of a set of independent curves ${\mathit{C}}_i$ (for 2D geometries) or surfaces ${\mathit{S}}_i$ (for 3D geometries):

$\mathbf{\Gamma }=\bigcup _{i=1}^I{\mathit{C}}_i\mathrm{}\mathrm{o}\mathrm{r}\mathrm{}\mathbf{\Gamma }=\bigcup _{i=1}^I{\mathit{S}}_i$ (1)

Figure 1a illustrates a B-representation 2D geometric model of a spanner while Figure 1b illustrates a B-representation 3D geometric model of a propeller.

Figure 1
B-representation geometric models: a) 2D model (Simpson et al., 2012). b) 3D model (Scott et al., 2013).

In general, the boundary of a 3D geometric object cannot be exactly described by a set of surfaces using the standard CAD technologies, namely NURBS surfaces. Gaps and overlaps in the surfaces intersections may occur for B-representations of topologically complex objects (Marussig & Hugues, 2018). For non-closed surface models, i.e., non-watertight B-representations, the union of the parametric surfaces is just an approximation $\tilde{\mathbf{\Gamma }}$ of the real boundary $\mathbf{\Gamma }$

$\mathbf{\Gamma }=\left(\tilde{\mathbf{\Gamma }}\bigcup {\mathbf{\Gamma }}_{gap}\right)\bigcap {\mathbf{\Gamma }}_{over}$ (2)

where ${\mathbf{\Gamma }}_{gap}$ are missing boundaries due to gaps and ${\mathbf{\Gamma }}_{over}$ the spurious boundary due to overlaps. Even though CAD technologies such as B-splines and NURBS surfaces are the standard CAD technologies and represent billions of dollars in investment, Non-watertight NURBS B-representations are inevitable in general CAD geometric modeling, limiting the design-through-analysis process, as domain mesh generation for FEM analysis requires watertight geometric models (Bazilevs et al., 2010).

2.1 Geometric mapping of B-splines curves

The geometric mappings of B-splines curves from their parametric spaces into the Cartesian space, i.e., $\mathbb{R}⟶{\mathbb{R}}^2$, are performed by linear combinations of basis functions and a polygonal of control points in ${\mathbb{R}}^2$. The univariate basis functions $N_{i,p}\left(\xi \right)$ are constructed from a vector $\mathrm{\Xi }=\left\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_{n+p+1}\right\}$ of parametric coordinates ${\xi }_i$ in $\mathbb{R}$. This vector is usually referred to as the knot vector, in which $n$ is the number of basis functions of order $p$ that can be constructed from the knot vector $\mathrm{\Sigma }$. For $p=0$, the B-spline basis functions are

$N_{i,0}\left(\xi \right)=\left\{\begin{array}{c}1\\ 0\end{array}\right.\mathrm{}\begin{array}{c}\mathrm{}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{}{\xi }_i\le \xi \le {\xi }_{i+1}\\ \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\mathrm{}$ (3)

For polynomial orders higher than zero, the basis functions are defined recursively from the Cox-de-Boor formula

$N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\mathrm{}\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}\left(\xi \right)$ (4)

The knot vectors are either uniform, i.e. equally spaced into the parametric coordinates, or non-uniform, in which the parametric coordinates are not equally spaced. The B-splines constructed from non-uniform knot vectors are referred to as non-uniform B-splines. For non-uniform knot vectors, repeated knot values may occur, i.e. multiplicity of knots. The multiplicity decreases the continuity of the geometric entities. Therefore, non-closed NURBS geometric entities adopt the multiplicity in their knot vectors to ensure the ends of the geometric entity. These entities are the most common case in CAD models. They are constructed from open knot vectors, i.e., knot vectors with multiplicity $k=p+1$ at the first and at the last values.

The B-spline curve $\mathit{C}\left(\xi \right)$ is constructed from the definition of the knot vector $\mathrm{\Xi }=\left\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_{n+p+1}\right\}$, a polynomial order $p$ for the basis and control points ${\mathit{B}}_i$, $i=1,\cdots ,n$ in ${\mathbb{R}}^2$. The B-spline curve mapping $\mathit{C}\left(\xi \right)$ is defined by the linear combination of the basis functions $N_{i,p}\left(\xi \right)$ and the control points ${\mathit{B}}_i$

$\mathit{C}\left(\xi \right)=\sum _{i=1}^n{N}_{i,p}\left(\xi \right){\mathit{B}}_i$ (5)

2.2 Geometric mapping of B-splines surfaces

The geometric mappings of B-splines surfaces from their parametric spaces into the Cartesian space, i.e., ${\mathbb{R}}^2⟶{\mathbb{R}}^3$, are performed by linear combinations of bivariate basis functions and a net of control points in ${\mathbb{R}}^3$. The bivariate basis functions $N_{ij,pq}\left(\xi ,\eta \right)=N_{i,p}\left(\xi \right)N_{j,q}\left(\eta \right)$ are obtained from the tensor product of the univariate basis functions $N_{i,p}\left(\xi \right)$ and $N_{j,q}\left(\eta \right)$ of orders $p$ and $q$, constructed from knot vectors $\mathrm{\Xi }=\left\{{\xi }_1,{\xi }_2,\cdots ,{\xi }_{n+p+1}\right\}$ and $\mathrm{{\rm H}}=\left\{{\eta }_1,{\eta }_2,\cdots ,{\eta }_{m+q+1}\right\}$, respectively. The parametric space of the surface is defined by the tensor product $\mathrm{\Xi }\otimes \mathrm{{\rm H}}$. The regular B-spline surface mapping $\mathit{S}\left(\xi \right)$ is defined by the linear combination of the basis functions $N_{ij,pq}\left(\xi ,\eta \right)$ and a net of control points ${\mathit{B}}_{ij}$, $i=1,\cdots ,n$, $j=1,\cdots ,m$ in ${\mathbb{R}}^3$.

$\mathit{S}\left(\xi \right)=\sum _{i=1}^n\sum _{j=1}^m{N}_{i,p}\left(\xi \right)N_{j,q}\left(\eta \right){\mathit{B}}_{ij}=\sum _{i=1}^n\sum _{j=1}^m{N}_{ij,pq}\left(\xi ,\eta \right){\mathit{B}}_{ij}$ (6)

2.3 NURBS curves and surfaces

The Non-Uniform Rational B-Splines (NURBS) curves and surfaces in ${\mathbb{R}}^n$ are achieved by a projective transformation of primitive B-splines curves and surfaces, respectively, in ${\mathbb{R}}^{n+1}$. From the algebraic point of view, the NURBS curves and surfaces may be presented in a similar way as presented for the B-Splines entities. Consequently, the NURBS curves are obtained by the linear combination of univariate rational basis functions and control points coordinates ${\mathit{B}}_i$, whereas the NURBS surfaces are defined by the linear combination of the bivariate rational basis functions and a control net ${\mathit{B}}_{ij}$

$\mathit{C}\left(\xi \right)=\sum _{i=1}^n{R}_i\left(\xi \right){\mathit{B}}_i$ (7)

$\mathit{S}\left(\xi ,\eta \right)=\sum _{i=1}^n\sum _{j=1}^m{R}_{ij}\left(\xi ,\eta \right){\mathit{B}}_{ij}$ (8)

in which $R_i$ and $R_{ij}$ are, respectively, the NURBS univariate and bivariate rational basis functions. These functions are defined by a rational weighting of the B-Splines basis functions $N_{i,p}$ (for curves) and $N_{ij,pq}$ (for surfaces). The $d+1$ components of the primitive B-spline control points are the weights, i.e., ${\left({\mathit{B}}_i\right)}_{\mathrm{d}+1}={\mathrm{\omega }}_i$ for curves and ${\left({\mathit{B}}_{ij}\right)}_{\mathrm{d}+1}={\mathrm{\omega }}_{ij}$ for surfaces. The NURBS basis functions are defined as

$R_i\left(\xi \right)=\frac{N_{i,p}\left(\xi \right){\mathrm{\omega }}_i}{\sum _{\widehat{i}=1}^n{N}_{\widehat{i},p}\left(\xi \right){\mathrm{\omega }}_{\widehat{i}}}$ (9)

$R_{ij}\left(\xi ,\eta \right)=\frac{N_{i,p}\left(\xi \right){M_{j,q}\left(\eta \right)\mathrm{\omega }}_{ij}}{\sum _{\widehat{i}=1}^n\sum _{\widehat{j}=1}^m{N}_{\widehat{i},p}\left(\xi \right){M_{\widehat{j},q}\left(\eta \right)\mathrm{\omega }}_{\widehat{i}\widehat{j}}}$ (10)

2.4 Triming curves and trimmed NURBS surfaces

The NURBS curves are also largely utilized together with the NURBS surfaces for the Boundary-representation of three-dimensional CAD models. These curves are utilized for the intersection definition between NURBS surfaces. To represent arbitrary boundaries of NURBS surfaces, trimming procedures must be applied. A trimmed NURBS surface is defined by a regular NURBS surfaces plus a set of trimming curves, ${\mathit{C}}^t\left(u\right)$, defined in the NURBS parametric space in ${\mathbb{R}}^2$, which result into closed loops that may or may not account for the boundary of the parametric space. The trimming curves are generally NURBS or B-splines entities, defined as

${\mathit{C}}^t\left(u\right)=\mathrm{}\left\{\begin{array}{c}{\xi }^t\left(u\right)\\ {\eta }^t\left(u\right)\end{array}\right\}=\sum _{i=1}^n{R}_i\left(u\right){\mathit{b}}_i$ (11)

where $u\in \mathbb{R}$ is the coordinates of the parametric space of the trimming curves and ${\mathit{b}}_i\in {\mathbb{R}}^2$ are the control points defined in the parametric space of the NURBS surface $\mathit{S}\left(\xi ,\eta \right)$. The loops defined by the sets of trimming curves define the visible parts of the parametric space $\mathrm{\Xi }\otimes \mathrm{{\rm H}}$ of the regular NURBS, which are mapped into the Cartesian space ${\mathbb{R}}^3$ and displayed in the CAD visualization tools. The resulting surface is the trimmed NURBS surface. Figure 2 illustrates the processes of converting a regular B-spline surface into a trimmed one by defining the visible area ${\mathcal{A}}^v$ in the parameter space of the regular surface through a trimming curve ${\mathit{C}}^t$.

Figure 2
Trimmed tensor product surface: a) regular surface. b) trimmed parameter space where a loop of trimming curves (thick line) specifies the visible aera ${\mathcal{A}}^v$. c) resulting trimmed surface. (Marussig & Hugues, 2018).

3 BOUNDARY ELEMENT MODELS

Unlikely finite element models, boundary element models provide numerical solutions for boundary value problems with a reduced-dimension for spatial discretization, i.e., only the surface boundaries must be discretized. This is achieved by rewriting the governing domain differential equations into Boundary Integral Equations (BIEs). The boundary integral equations can be deduced from the knowledge of the fundamental solution of the differential operator and integration by parts operations starting from a weighted integral form of the domain differential equations (Brebbia et al. 1984). The displacement BIE of linear elasticity is presented in the following, together with the numerical procedure for its solution based on the collocation version of the boundary element method.

3.1 Boundary integral equations of linear elasticity

The displacement balance equations of a given linear elastic solid with domain$\mathbf{}\mathbf{\Omega }$ and boundary $\mathbf{\Gamma }$, such as the one illustrated in Figure 3a, can be represented in terms of the BIEs. For general mixed boundary value problems, the boundary is decomposed into two complementary parts $\mathbf{\Gamma }={\mathbf{\Gamma }}_{\mathcal{D}}\bigcup {\mathbf{\Gamma }}_{\mathcal{N}}$, ${\mathbf{\Gamma }}_{\mathcal{D}}\bigcap {\mathbf{\Gamma }}_{\mathcal{N}}=\mathit{\varnothing }$, where ${\mathbf{\Gamma }}_{\mathcal{D}}$ is the part where Dirichlet boundary conditions are prescribed, i.e., the displacements are known, and ${\mathbf{\Gamma }}_{\mathcal{N}}$ is the part where Neumann boundary conditions are prescribed, i.e., the tractions are known.

Figure 3
Problem definition: a) boundary value problem. b) boundary source and field points.

The displacement BIE can be derived from the Navier-Cauchy equations of elasticity and, in the absence of body forces, results

$\mathbf{C}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\mathbf{u}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)={\int }_{\mathbf{\Gamma }}^{}\mathbf{U}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x})\mathbf{t}\left(\mathbf{x}\right)\mathrm{d}\mathrm{\Gamma }-\mathcal{C}{\int }_{\mathbf{\Gamma }}^{}\mathbf{T}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)\mathbf{u}\left(\mathbf{x}\right)\mathrm{d}\mathrm{\Gamma },$ (12)

where ${\mathbf{X}}^{\mathrm{\text{'}}}$ and $\mathbf{X}$ indicates the source and field points, respectively (see Fig. 3b). $\mathbf{u}\left({\mathbf{X}}^{\mathrm{\text{'}}}\right)$ is the displacement at ${\mathbf{X}}^{\mathrm{\text{'}}}$, while $\mathbf{u}\left(\mathbf{X}\right)$ and $\mathbf{t}\left(\mathbf{X}\right)$ are the displacements and tractions boundary fields, i.e., for boundary field points $\mathbf{X}\in \mathbf{\Gamma }$.

The Kelvin fundamental solution for displacements $\mathbf{U}({\mathbf{X}}^{\mathrm{\text{'}}},\mathbf{X})$ and tractions $\mathbf{T}({\mathbf{X}}^{\mathrm{\text{'}}},\mathbf{X})$ are given by

$\mathbf{U}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)={\mathrm{U}}_{ij}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right){\mathbf{e}}_i\otimes {\mathbf{e}}_j\mathbf{T}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)={\mathrm{T}}_{ij}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right){\mathbf{e}}_i\otimes {\mathbf{e}}_j,$ (13)

in which the components ${\mathrm{U}}_{ij}$ and ${\mathrm{T}}_{ij}$ of the fundamental solution for the Navier-Cauchy equations in ${\mathbb{R}}^3$ are

${\mathrm{U}}_{ij}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)=\frac{1}{16\pi \mu \left(1-\nu \right)r}\left[\left(3-4\nu \right){\delta }_{ij}+\frac{\partial r}{\partial x_i}\frac{\partial r}{\partial x_j}\right]$ (14)

${\mathrm{T}}_{ij}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)=\frac{-1}{8\pi \left(1-\nu \right)r^2}\left\{\frac{\partial r}{\partial n}\left[\left(1-2\nu \right){\delta }_{ij}+3\frac{\partial r}{\partial x_i}\frac{\partial r}{\partial x_j}\right]-\left(1-2\nu \right)\left(n_j\frac{\partial r}{\partial x_i}-n_i\frac{\partial r}{\partial x_j}\right)\right\}.$ (15)

In Eqs (14) and (15), $r=‖\mathbf{r}‖=‖\mathbf{x}-{\mathbf{x}}^{\mathrm{\text{'}}}‖$ is the distance between the source and the field point, $n_i$ stands for the components of the normal outward vector $\mathit{n}\left(\mathbf{x}\right)$ and $\partial r/\partial n=\nabla r\bullet \mathbf{n}=\partial r/\partial x_i{n}_i$ is the normal derivative of $r$. $\mathcal{C}\int$ Indicates that the integral is evaluated in the Cauchy Principal Value sense and $\mathbf{C}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\mathbf{u}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)$ is the free term arising from the singular integral (Aliabadi, 2002). For a non-watertight B-representation $\tilde{\mathbf{\Gamma }}$, the BIE (12) is not an exact mathematical integral representation of the problem. However,

$\mathbf{C}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\mathbf{u}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\approx {\int }_{\tilde{\mathbf{\Gamma }}}^{}\mathbf{U}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x})\mathbf{t}\left(\mathbf{x}\right)\mathrm{d}\mathrm{\Gamma }-\mathcal{C}{\int }_{\tilde{\mathbf{\Gamma }}}^{}\mathbf{T}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\right)\mathbf{u}\left(\mathbf{x}\right)\mathrm{d}\mathrm{\Gamma },$ (16)

as long as the missing boundaries ${\mathbf{\Gamma }}_{gap}$ due to gaps and spurious boundaries ${\mathbf{\Gamma }}_{over}$ due to overlaps are small enough compared to $\tilde{\mathbf{\Gamma }}$, i.e., $‖{\mathbf{\Gamma }}_{gap}‖\ll \mathrm{}‖\tilde{\mathbf{\Gamma }}‖$ and $‖{\mathbf{\Gamma }}_{over}‖\ll \mathrm{}‖\tilde{\mathbf{\Gamma }}‖$. Thus, enforcing the equality in (16) will not produce significant errors, as long as the continuity of the boundary fields is not required by the governing boundary integral equations.

3.2 Numerical approximations and solution by collocation point

Approximations are required for obtaining numerical solutions of the previous presented BIE by the classical weighted residue collocation method. Thus, the boundary $\tilde{\mathbf{\Gamma }}=\bigcup _{i=1}^I{\mathit{S}}_i$ must be discretized. Each NURBS surface ${\mathit{S}}_i$ are then discritezed into ${NE}_{\left(i\right)}$ linear (triangular or quadrilateral) boundary elements, over which displacements, tractions, and the boundary geometry are approximated regarding the isoparametric concept. Introducing the numerical approximations into Eq. (12) results

$\mathbf{r}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)=\mathbf{C}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\mathbf{u}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)+\sum _{i=1}^I{\int }_{{\mathbf{S}}_i}^{}\mathbf{T}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\left(\xi ,\eta \right))\mathbf{u}\left(\mathbf{x}\left(\xi ,\eta \right)\right)\mathrm{d}{S}_i-\sum _{i=1}^I{\int }_{{\mathbf{S}}_i}^{}\mathbf{U}\left({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\left(\xi ,\eta \right)\right)\mathbf{t}\left(\mathbf{x}\left(\xi ,\eta \right)\right)\mathrm{d}{S}_i$ (17)

$\mathbf{r}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)=\mathbf{C}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)\mathbf{u}\left({\mathbf{x}}^{\mathrm{\text{'}}}\right)+\sum _{i=1}^I\sum _{e=1}^{{NE}_{\left(i\right)}}\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\mathbf{T}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\left({\xi }_1,{\xi }_2\right))N_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right){\mathrm{J}}^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)\mathrm{d}{\xi }_1\mathrm{d}{\xi }_2{\mathbf{u}}_{\alpha }^{e\left(i\right)}-\sum _{i=1}^I\sum _{e=1}^{{NE}_{\left(i\right)}}\underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\mathbf{U}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x}\left({\xi }_1,{\xi }_2\right))N_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right){\mathrm{J}}^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)\mathrm{d}{\xi }_1\mathrm{d}{\xi }_2{\mathbf{t}}_{\alpha }^{e\left(i\right)}$ (18)

where $\mathbf{r}\left({\mathbf{X}}^{\mathrm{\text{'}}}\right)$ is the boundary integral equation residue, $I$ indicates the number of NURBS surfaces ${\mathit{S}}_i$ (NURBS patches) describing the boundary $\tilde{\mathbf{\Gamma }}$ and ${NE}_{\left(i\right)}$ the number of boundary elements utilized for the discretization of ${\mathit{S}}_i$. The mapping ${\mathit{S}}_i\left(\xi ,\eta \right)$, given by (8), is adopted for the generation of the nodal coordinates of the conventional linear elements, as detailed in the next subsection. The basis functions $N_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)$ of the element $e\left(i\right)\in {\mathit{S}}_i$ are linear Lagrange polynomials, which are applied for both the geometric mapping $\mathbf{x}\left({\xi }_1,{\xi }_2\right)$ and the fields approximations ${\mathbf{u}}^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)=N_{\alpha }^{e\left(i\right)}{\mathbf{u}}_{\alpha }^{e\left(i\right)}$ and ${\mathbf{t}}^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)=N_{\alpha }^{e\left(i\right)}{\mathbf{t}}_{\alpha }^{e\left(i\right)}$, in which ${\mathbf{u}}_{\alpha }^{e\left(i\right)}$ and ${\mathbf{t}}_{\alpha }^{e\left(i\right)}$ are, respectively, the coefficients of the displacement and tractions approximations for the element $e\left(i\right)\in {\mathit{S}}_i$. ${\mathrm{J}}^{e\left(i\right)}\left(\xi ,\eta \right)$ is the Jacobian of the Gaussian-to-physical space mapping of the linear isoparametric elements. Notice that ${\xi }_1$ and ${\xi }_2$ are the coordinated of the Gaussian space of the boundary elements and must not be confused with the first and secont knot values of $\mathrm{\Xi }$. The usual BEM system of equations

$\mathbf{H}\mathbf{u}=\mathbf{G}\mathbf{T}$ (19)

is thus obtained by imposing $\mathbf{r}\left({\mathbf{X}}_n^{\mathrm{\text{'}}}\right)=0$ for a set of collocation points ${\mathbf{X}}_n^{\mathrm{\text{'}}}$ , $n=1,\cdots ,N$ , where $N$ is the number of unknown coefficents of the numerical approximations. The dense (full-populated) matrices $\mathbf{H}$, $\mathbf{G}$ arise from the integral of the kernels $\mathbf{T}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x})$ and $\mathbf{U}({\mathbf{x}}^{\mathrm{\text{'}}},\mathbf{x})$, respectively, and $\mathbf{u}$, $\mathbf{t}$ are vectors collecting the coefficients ${\mathbf{u}}_{\alpha }^{e\left(i\right)}$, ${\mathbf{t}}_{\alpha }^{e\left(i\right)}$ of all boundary elements, accounting for their connectivity. After the boundary conditions imposition, the unknown coeficients of $\mathbf{u}$ and $\mathbf{t}$ are collected into a vector $\mathbf{x}$ and the final system $\mathbf{A}\mathbf{x}=\mathbf{b}$ can thus be solved. The weak-singular integrals in (14) are integrated in polar coordinates, and the singularity cancels with the Jacobian of the transformation. The strong-singular integrals on the other hand are regularized by the Guiggiani Method (Guiggiani and Gigante, 1990). The remaining regular integrals are performed with adaptative Gaussian quadrature and sub-elementation integration (Cordeiro, 2018). Discontinuous boundary elements, such as those presented in Aliabadi (2002), are adopted in regions of geometric or traction discontinuities, and also together with a special collocation strategy to ensure an adequate collocation points distance, when using the surface-independent discretization for the NURBS surfaces ${\mathit{S}}_i$ , $i=1,\cdots ,I$.

4 Patch independent discretizations

The discretization in the proposed numerical framework accounts for the following details: (i) Independent discretization of each NURBS surface into linear elements, based on the independent parametrizations ${\mathit{S}}_i\left(\xi ,\eta \right)$; (ii) Generalized discontinuous elements at the edges of the visible areas of the NURBS parametric spaces are detected by a Level Set function where the offset collocation strategy is adopted; (iii) The discontinuities in the displacement field between surfaces are small enought due to the fact that the integral kernels are non-null at all points of the ${\mathbb{R}}^3$ Cartesian space, resulting full-populated matrices that does not allow the independent motion of the surfaces; (iv) The offset collocation strategy avoids singularities and near-singularities due to coincident or small distance collocation points; (v) Regular NURBS surfaces are discretized into linear quadrilateral elements while trimmed surfaces can be discretized with both linear triangular and quadrilateral elements; (vi) The basis functions of the generalized discontinuous elements are built by imposing the Kronecker property at the offseted parametric coordinates of the collocation points, retaining the physical meaning of the coefficients of the elemental field’s aproximations.

The first step toward the surface discretization generation is to build ${\mathbb{R}}^2$ meshes for the parametric spaces of the NURBS surfaces. It is worth to emphasize that the discretization of each parametric space is performed independently. The nodes at the edges of the parametric spaces are identified by a level-set function and the respective elements are transformed into generalized discontinuous boundary elements where the offset collocation strategy is employed. Finally, based on the NURBS surface geometric mapping, Eq. (8), the nodal coordinates in the ${\mathbb{R}}^3$ Cartesian space are obtained.

4.1 Discretization of regular NURBS patches

The discretization of the parametric space of a regular NURBS surface is trivial. In the proposed numerical framework, the discretization is performed with linear quadrilateral boundary elements, as illustrated in Fig. 4.

Figure 4
Discretization of regular NURBS patches: a) NURBS parameter space. b) discretization of the parameter space. c) resulting discretization in ${\mathbb{R}}^3$.

4.2 Discretization of trimmed NURBS patches

The discretization at the visible areas of the parametric space of a trimmed NURBS surface requires the knowledge of the number of trimming loops $l=1,\cdots ,L$, the number o trimming curves $c=1,\cdots ,C\left(l\right)$ per loop $l$ and their mathematical description, i.e., the mapping ${\mathit{C}}_{c\left(l\right)}^t\left(u\right)$. The trimming curves and the edges of the parametric space are discretized into a finite number of line segments. The trimmed NURBS surfaces are classified into two types for discretization purposes: the inner and the outer trimmed NURBS surfaces. In the former case, the meshed region is delimited only by the line segments of the trimming curves. In the latter case, the line segments of both trimming curves and the parametric space edges delimit the meshed region. Based on this geometric definition, the discretization of the visible areas ${\mathcal{A}}^v$ of the parametric space into triangular or quadrilateral linear elements is performed based on bidimensional unstructured mesh generation algorithms. Fig. 5a and 5b illustrate the discretization of the parametric spaces of both inner and outer trimmed NURBS surfaces, respectively. The visible areas ${\mathcal{A}}^v$ are presented in grey. ${\mathit{C}}_{c\left(l\right)}^t$ indicates the trimming curve number $c\in l$. Fig. 5 also illustrates the ${\mathbb{R}}^3$ discretization of the surface in the physical space.

Figure 5
Discretization of trimmed NURBS patches: a) inner discretization. b) outer discretization.

4.3 Generalized discontinuous elements and the offset collocation

Collocation points in the same geometric position or small distance collocation between points leads to numerical singularities in any collocation based method, such as the collocation boundary integral equation method. The offset collocation strategy adopted herein avoids collocation points closely positioned. The collocation points are obtained by an offset of the edge nodes of the parametric mesh into the interior of the visible area of the NURBS parametric space, avoiding collocation at the edges or trimming curves. The edge nodes are identified by a signed distance function $\mathrm{\varphi }\left({\mathit{\xi }}_n,{\mathbf{\Gamma }}^v\right)$ (or, level set function) defined in the visible area ${\mathcal{A}}^v$ of the parametric spaces $\mathrm{\Xi }\otimes \mathrm{{\rm H}}$ as

$\mathrm{\varphi }\left({\mathit{\xi }}_n,{\mathbf{\Gamma }}^v\right)=\left\{\begin{array}{c}\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}\left({\mathit{\xi }}_n,{\mathbf{\Gamma }}^v\right),\\ 0\\ -\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}\left({\mathit{\xi }}_n,{\mathbf{\Gamma }}^v\right),\end{array}\right.\mathbf{}\begin{array}{c}\mathrm{}\mathrm{i}\mathrm{f}\mathrm{}{\mathit{\xi }}_n\in {\mathcal{A}}^v\\ \mathrm{}\mathrm{i}\mathrm{f}\mathrm{}{\mathit{\xi }}_n\in {\mathbf{\Gamma }}^v\\ \mathrm{}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$ (20)

where ${\mathit{\xi }}_n={\left\{{\xi }_n,{\eta }_n\right\}}^{\mathrm{T}}$ are the parametric coordinates of the $n=1,\cdots ,N\left(i\right)$ nodes adopted in the discretization of ${\mathit{S}}_i$. The boundary of the visible area ${\mathcal{A}}^v$ is ${\mathbf{\Gamma }}^v$ and $\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}\left({\mathit{\xi }}_n,{\mathbf{\Gamma }}^v\right)$ is the Eucliadian distance function between the boundary ${\mathbf{\Gamma }}^v$ and the node coordinate ${\mathit{\xi }}_n$ in the parametric space $\mathrm{\Xi }\otimes \mathrm{{\rm H}}$, which can be numerically computed regarding the parametrizations of the trimming curves ${\mathit{C}}_{c\left(l\right)}^t\left(u\right)$ and small parametric increments $\mathrm{\Delta }u$. The elements that intercept the boundary of the visible areas ${\mathcal{A}}^v$ are set as generalized discontinuous boundary elements. The triangular and quadrilateral linear generalized discontinuous boundary elements adoped herein are illustrated in Fig. 6a, while Figs. 6b illustrates their use for discretization of the visible areas of the parametric spaces of both regular and trimmed NURBS surfaces.

Figure 6
Generalized discontinuous elements: a) collocation points in the Gaussian space. b) collocation points in the parametric space of a regular NURBS. c) collocation points in the parametric space of a trimmed NURBS.

The basis functions $N_{\alpha }^{e\left(i\right)}$ of the generalized discontinuous elements are built imposing the Kronecker delta property for the Lagrange polynomials at the same parametric coordinates of the offseted collocation points. This ensure that the coefficients ${\mathbf{u}}_{\alpha }^{e\left(i\right)}$ and ${\mathbf{t}}_{\alpha }^{e\left(i\right)}$ of the fields approximations have physical meaning. The term “generalized” was introduced in the present work to emphasize the fact that the collocation points and anchors for the kronecker imposition can be chosen arbitraly inside the domain of the element (see Fig. 6a). The geometric mapping of the generalized discontinuous elements is defined by $\mathbf{x}\left({\xi }_1,{\xi }_2\right)=S_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right){\mathbf{x}}_{\alpha }^{e\left(i\right)}$, where $S_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)$ is the original basis with nodal Kronecker imposition and ${\mathbf{x}}_{\alpha }^{e\left(i\right)}$ the nodal coordinates of the element $e\in {\mathit{S}}_i$. On the other hand, the basis functions $N_{\alpha }^{e\left(i\right)}\left({\xi }_1,{\xi }_2\right)$ for the fields of the generalized elements are computed by the imposition of the Kronecker property:

$N_{\alpha }\left({\xi }_1,{\xi }_2\right)=c_1^{\alpha }+c_2^{\alpha }{\xi }_1+c_3^{\alpha }{\xi }_2\mathbf{}{N}_{\alpha }\left({{\xi }_1}^{\beta },{{\xi }_2}^{\beta }\right)=\left\{\begin{array}{c}1\\ 0\end{array}\begin{array}{c}\mathrm{}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{\alpha }=\mathrm{\beta }\\ \mathrm{}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{\alpha }\ne \mathrm{\beta }\end{array}\right.\mathbf{}\Rightarrow c_j^{\alpha },\mathbf{}\mathrm{\alpha }=\mathrm{\beta }=\mathrm{j}=1,\cdots ,3$ (21)

for linear triangular elements and

$N_{\alpha }\left({\xi }_1,{\xi }_2\right)=c_1^{\alpha }+c_2^{\alpha }{\xi }_1+c_3^{\alpha }{\xi }_2+c_4^{\alpha }{\xi }_1{\xi }_2\mathbf{}{N}_{\alpha }\left({{\xi }_1}^{\beta },{{\xi }_2}^{\beta }\right)=\left\{\begin{array}{c}1\\ 0\end{array}\begin{array}{c}\mathrm{}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{\alpha }=\mathrm{\beta }\\ \mathrm{}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{}\mathrm{\alpha }\ne \mathrm{\beta }\end{array}\right.\mathbf{}\Rightarrow c_j^{\alpha },\mathbf{}\mathrm{\alpha }=\mathrm{\beta }=\mathrm{j}=1,\cdots ,4$ (22)

for linear quadrilateral elements.

5 RESULTS

Three examples are presented to investigate the performance of non-watertight boundary element models, built from B-representations regarding the surface-independent discretization approach, in linear elasticity three-dimensional problems. The first example deals with a cube cut by a cylindrical surface with available analytical solution. A convergence study regarding small gaps in non-watertight models was performed. The second and the third examples deal with the numerical analysis of models built from complex B-representation CAD models. In the final example, the analysis model is built from a non-watertight geometry, which is not an option for finite element models.

5.1 Cube cut by a cilyndrical surface - convergence test

The first example corresponds to a cube of length$L=1$ (length unity), cut by a concentric cylindrical surface of radius $r=0.15$ (length unity) as illustrated in Fig. 7. The displacement components normal to the faces $x=0$, $y=0$ and $z=0$ are restricted, and a unity traction $t_z=1$ (stress unity) acts on the face $z=1$.

Figure 7
Cube cut by a cylindrical surface (red). The unspecified boundary conditions are zero tractions (Peng, 2016).

The analytical displacement solution for such a problem is given by

${\mathrm{u}}_{\mathrm{x}}\left(\mathbf{x}\right)=-\frac{\nu x}{E}\mathit{},\mathit{}{\mathrm{u}}_{\mathrm{y}}\left(\mathbf{x}\right)=-\frac{\nu y}{E},\mathit{}{\mathrm{u}}_{\mathrm{z}}\left(\mathbf{x}\right)=\frac{z}{E}$ (23)

The elastic constants adopted for the problem were $E=1000$ (stress unity) and $\nu =0.3$. From the seven surfaces that constitute the B-representation of the problem, the two trimmed faces $z=0$ and $z=1$ were discretized with triangular elements, while the remaining ones were discretized with quadrilateral elements. Three levels of discretizations were adopted for a convergence study: Level I, Level II and Level III. The Levels I, II and III discretizations comprise 704, 1328, 2924 elements and 772, 1446, 2995 nodes, respectively. Six models were then generated: three with the trimmed surfaces matching the cylindrical surface and three with intentioned gaps at the intersection of the trimmed and cylindrical surfaces. Fig. 8 illustrates a frontal view of the watertight and non-watertight meshes for the most refined level.

Figure 8
Discrete BEM models: a) watertight model. b) non-watertight model.

Two convergence analyses were performed, one for the watertight models and another for the non-watertight models. The influence of the gaps can be evaluated regarding the convergence in terms of the L2-norm of the boundary displacement error. The L2-norm of boundary displacement error can be computed as

${‖{\mathbf{e}}_{\mathrm{u}}‖}_{L2}=\sqrt{\frac{{\int }_{\tilde{\mathbf{\Gamma }}}^{}\left(\tilde{\mathbf{u}}-\mathbf{u}\right)\bullet \left(\left(\tilde{\mathbf{u}}-\mathbf{u}\right)\right)\mathrm{d}\mathbf{\Gamma }}{{\int }_{\tilde{\mathbf{\Gamma }}}^{}\mathbf{u}\bullet \mathbf{u}\mathrm{d}\mathbf{\Gamma }}}$ (24)

in which $\tilde{\mathbf{\Gamma }}=\bigcup _{i=1}^I{\mathit{S}}_i$ is the approximated boundary, $\tilde{\mathbf{u}}$ the approximated displacement solution and $\mathbf{u}$ the analytical displacement solution. The convergence results in terms of the L2-norm of boundary displacement error are presented in Fig. 9a concerning the number of collocation points of the models. Besides, Fig. 9b shows the exponential increase that occurs in ${‖{\mathbf{e}}_{\mathrm{u}}‖}_{L2}$ as the ratio of the gap areas per the total boundary area increases.

Figure 9
${‖{\mathbf{e}}_{\mathrm{u}}‖}_{L2}$ convergence results: a) watertight vs. non-watertight models. b) error vs. gap sizes.

As expected, the error reduces when the number of collocation points increases in both convergence analyses. However, the gaps do affect the error ${‖{\mathbf{e}}_{\mathrm{u}}‖}_{L2}$ in one order of magnitude, i.e., from ${10}^{-4}$ to ${10}^{-3}$. Even thought, global errors with magnitudes ${10}^{-3}$ are still acceptable for pratical engineering analysis purposes.

To illustrate that the numerical results are satisfatory even in the presence of small gaps, Fig. 10 brings the results of the norm of the displacement vector for both watertight and non-watertight models with refinement Level III. Note that the results are identical considering the decimal digits presented.

Figure 10
Displacement boundary field results: a) watertight model. b) non-watertight model.

This example is considered to be a benchmark for the validation of non-watertight discrete BEM models of linear elasticity. In the two next examples, non-watertight discrete BEM models are built from complex B-representation with the surface-independent discretization approach. The B-representation may even be non-watertight B-representations.

5.2 Base-plate - Comparative analysis

The second example deals with the analysis of a base-plate for structural applications. The multiple-path B-representation NURBS model was obtained from an Initial Graphics Exchange Specification (IGES) file, available in the GrabCAD library. Fig. 11 illustrates the visualization of the B-representation geometry obtained from the IGES Tollbox in MATLAB and the most refined discrete BEM model generated with the surface-independent discretization algorithm.

Figure 11
B-representation models: a) B-representation geometry. b) non-watertight discrete BEM model.

The B-representation model is composed of 36 NURBS surfaces. There are 30 regular NURBS surfaces of polynomial orders $p=1$ and $q=2$, and 6 trimmed NURBS surfaces of polynomial orders $p=1$ and $q=1$ . From the 6 trimmed surfaces, 4 are trimmed by 2 circular curves ($p=2$), and the other 2 are trimmed by 8 circular curves ($p=2$). Two levels of discretization were adopted for the discrete BEM models: Level I and Level II. The Levels I and II discretizations are composed of 4586 and 12367 elements, and 5586 and 13874 nodes, respectively. The BEM models were built without any concern about conforming meshes at the intersections of the surfaces, thus resulting non-watertight discrete models. Since the B-representation is a watertight geometric model, a volumetric Finite Element Model (FEM) of the problem was easily created in the commercial software Ansys. The volumetric mesh from Ansys is composed of 15322 tetrahedral quadratic elements of the type Solid186, resulting 27526 nodes. The elastic constants adopted for the analysis were $E=1000$ (stress unity) and $\nu =0.3$. The boundary conditions imposed in the model were null prescribed displacements,$\mathrm{}\stackrel{-}{\mathbf{u}}=0$, in the left-half surface of the eight circular holes and prescribed tractions $\stackrel{-}{\mathbf{t}}={\left\{\begin{array}{ccc}10& 0& 0\end{array}\right\}}^{\mathrm{T}}$ (stress unity) for the right-half surface of the central axial hole. Fig. 12 illustates the boundary conditions of the problem.

Figure 12
Discretization and boundary conditions: a) prescribed $\stackrel{-}{\mathbf{u}}=0$. b) prescribed $\stackrel{-}{\mathbf{t}}={\left\{\begin{array}{ccc}10& 0& 0\end{array}\right\}}^{\mathrm{T}}$.

The structural BEM and FEM analysis produced similar displacement fields. For the level I BEM model, the maximum values of displacement components were ${\mathrm{u}}_{\mathrm{x}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=4.72828$, ${\mathrm{u}}_{\mathrm{y}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=2.05836$ and ${\mathrm{u}}_{\mathrm{z}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=0.50842$ (length unities). The level II BEM model on the other hand produced the following maximum values of displacement components ${\mathrm{u}}_{\mathrm{x}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=4.86869$, ${\mathrm{u}}_{\mathrm{y}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=2.13497$ and ${\mathrm{u}}_{\mathrm{z}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=0.53051$ (length unities). The differences of theses and the FEM results are ${\mathrm{d}\mathrm{u}}_{\mathrm{x}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=1.86\%$, d${\mathrm{u}}_{\mathrm{y}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=8.56\%$ and d${\mathrm{u}}_{\mathrm{z}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=9.65\%$ for the level I BEM model and ${\mathrm{d}\mathrm{u}}_{\mathrm{x}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=1.05\%$, d${\mathrm{u}}_{\mathrm{y}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=12.60\%$ and d${\mathrm{u}}_{\mathrm{z}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=12.30\%$ for the level II BEM model. The undeformed and deformed configurations for the level II non-watertight BEM and Ansys FEM models in scale 1:1 are illustrated in Figure 13.

Figure 13
Undeformed/deformed mesh configurations: a) non-watertight BEM model. b) watertight Ansys FEM model.

Fig. 14 presents a qualitative displacement boundary fields evaluation comparing the most refined non-watertight BEM model with the Ansys FEM model.

Figure 14
Displacement boundary field results: a) non-watertight BEM model. b) watertight Ansys FEM model.

5.3 Crane hook

The last example deals with the analysis of a crane hook. The multiple-path B-representation model was obtained from an Initial Graphics Exchange Specification (IGES) file available in the GrabCAD library. The B-representation of the problem is a complex non-watertight one, which introduces severe difficulties for volumetric mesh generation. A surface mesh with non-conforming meshes at the surfaces intersections can be easily created with the surface-independent discretization algorithm on the other hand. Fig. 15a illustrates the visualization of the B-representation geometry obtained with the IGES Tollbox in MATLAB and a zoom in the non-watertight part of the B-representation is presented in Fig. 15b.

Figure 15
B-representation geometry: a) full B-representation geometry. b) non-watertight part of the geometry.

The B-representation model of the hook is composed of 69 NURBS surfaces. From the 69 surfaces, 32 are regular NURBS surfaces of polynomial orders $p=5$ and $q=5$, 8 are regular NURBS surfaces of polynomial orders $p=5$ and $q=4$, and 4 are regular NURBS surfaces of polynomial orders $p=1$ and $q=2$. From the remaining 25 NURBS surfaces, 20 are trimmed NURBS surfaces of polynomial orders $p=5$ and $q=5$, 2 are trimmed NURBS surfaces of polynomial orders $p=5$ and $q=4$, and 3 are trimmed NURBS surfaces of polynomial orders $p=1$ and $q=1$ . In the parameter space of the 25 trimmed surfaces, a total of 93 trimming curves of polynomial orders $p=1$, $p=2$ and $p=5$ are defined. A single refined discrete BEM model is generated with the surface-independent discretization algorithm. The discretization is composed of 9623 elements and 11614 nodes. The BEM model was built without any concern about conforming meshes at the intersections of the surfaces, thus avoiding complications due to the non-watertight nature of the geometric B-representation. Fig. 16a illustrates the discrete BEM model and a zoom in the non-watertight part of the discretization is presented in Fig. 16b.

Figure 16
Discrete BEM model: a) full BEM model. b) non-watertight part of the model.

The elastic constants adopted for the analysis were $E=100000$ (stress unity) and $\nu =0.3$. The boundary conditions imposed in the model were null prescribed displacements,$\mathrm{}\stackrel{-}{\mathbf{u}}=0$, in the bottom surface of the hook suport and prescribed tractions $\stackrel{-}{\mathbf{t}}={\left\{\begin{array}{ccc}0& -10& 0\end{array}\right\}}^{\mathrm{T}}$ (stress unity) in the central part of the body of the hook, as illustrated in Fig. 17.

Figure 17
Discretization and boundary condition: a) prescribed $\stackrel{-}{\mathbf{u}}=0$ (in red). b) prescribed $\stackrel{-}{\mathbf{t}}={\left\{\begin{array}{ccc}0& -10& 0\end{array}\right\}}^T$(in blue).

Fig. 18 presents the results of the norm of the boundary displacement field and the undeformed/deformed configurations in scale 1:50.

Figure 18
Crane-hook results: a) norm of the boundary displacement field. b) undeformed/deformed configurations.

Since the B-representation of the problem is a complex non-watertight one, the volumetric mesh generation needed for developing a finite element (FE) model would require CAD repaire algotithms to close the gaps of the geometry (Coll et. al. 2014). Therefore, no FE model was created for this example. Note that even though there is a visiable gap in the model, a reasonable solution has been achieved with the BEM model.

6 CONCLUSION

The present study investigates how unsought gaps and overlaps in non-watertight boundary element models affect the numerical response of linear elasticity problems. A surface-independent discretization approach was developed to achieve flexibility in model generation without any concern about surface-to-surface mesh compatibility. Linear quadrilateral elements were adopted for the discretization of regular surfaces, while both triangular and quadrilateral elements can be adopted for the discretization of trimmed surfaces. Generalized discontinuous elements at the edges of the visible areas of the NURBS parametric spaces are detected by a Level Set function, and an offset collocation strategy was adopted for the nodes at the edges of the visible areas of the NURBS parametric spaces, avoiding singularities in the collocation BEM equations. The influence of the gaps in the convergence of the L2-norm of the boundary displacement error was verified in a 3D example with an available analytical solution. It is also shown that the L2-norm of the boundary displacement error increases exponentially with the increase in the size of the gaps. A second example with an available finite element solution was analyzed regarding a non-watertight BEM model for qualitative boundary field validation. Finally, a complex non-watertight B-representation of a crane hook consisting of 44 regular NURBS surfaces and 25 trimmed NURBS surfaces was analyzed. The obtained results have pointed out that, as long as the gaps and overlaps are small enough, BEM models built up from non-watertight geometries may produce valuable solutions for practical purposes. Since the surface gaps and overlaps, inherent to practical geometric design technics, are usually very small, one can adopt the BEM as a plausible choice for the structural analysis, avoiding redesigning the solid in CAD software. Future investigations may include 3D path-discontinuous isogeometric convergence and sensitivity analysis for non-watertight models.

References

• Aliabadi, M.H.F. (2002). The Boundary Element Method, volume 2: Appliecation in solids and structures, John Wiley Sons (Chichester).
• Bazilevs, Y., Calo, V.M., Cottrell, J.A., Evans, J.A., Hugues, T.J.R., Lipton, S., Scott, M.A., Sedergerg, T.W. (2010). Isogeometric analysis using T-splines, Computer Methods in Applied Mechanics and Engineering 199: 229-263.
• Belytschko, T., Lu, Y.Y., Gu, L. (1994). Element-free Galerkin methods, International journal for numerical methods in engineering 37: 229-256.
• Brebbia, C.A., Telles, J.C.F., Wrobel, L.C. (1984). Boundary element techniques: theory and applications in engineering, Springer (New York).
• Coll, A., Dadvand, P., Oñate, E. (2014). Robust volume mesh generation for non-watertight geometries, Monograph CIMNE.
• Cordeiro, S.G.F. (2018). Contributions to fracture and fatigue analysis of three-dimensional components by the dual boundary element method, Ph.D. Thesis (in Portuguese), University of Sao Paulo, Brazil.
• Cottrell, A.J., Hugues, T.J.R., Bazilevs, Y. (2009). Isogeometric analysis: Toward Integration of CAD and FEA, John Wiley Sons (Chichester).
• Daumas, G.O., Lins, R.M., Cordeiro, S.G.F., Monteiro, F.A.C. (2022). Non-watertight geometries and boundary element models, Proceedings of the 8th International Symposium on Solid Mechanics, Campinas-SP, Brazil.
• Guiggiani, M., Gigante, A. (1990). A general algorithm for multidimensional Cauchy principal valeu integrals in the boundary element method, Journal of Applied Mechanics 57(4): 906-915.
• Hsiao, G.C., Wendland, W.L. (2008). Boundary integral equations, Springer-Verlag (Berlin).
• Hugues, T.J.R., Cottrell, J.A., Bazilevs, Y. (2005). Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194(39-41) : 4135-4195.
• Lian, H., Simpson, R.N., Bordas, S.P.A. (2013). Stress analysis without meshing: isogeometric boundary element method, Proceedings of the institution of civil engineers: Engineering and computational mechanics 166: 88-99.
• 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: 458-488.
• Marussig, B., Hugues, T.J.R. (2018). A review of trimming in isogeometric analysis: challenges, data exchange and simulation aspects, Archives of Computational Methods in Engineering 25(4): 1059-1127.
• Oliveira, H.L., Andrade, H.C., Leonel, E.D. (2020). An isogeometric boundary element approach for topology optimization using the level set method, Applied Mathematical Modelling 84: 536-553.
• Peng, X. (2016). Isogeometric boundary element methods for linear elastic fracture mechanicsContributions to fracture, Ph.D. Thesis, Cardiff University, United Kingdom.
• Peng, X., Atroshchenko, E. Kerfriden, P., Bordas, S.P.A. (2017). Linear elastic fracture simulation directly from CAD: 2D NURBS-based implementation and rode of tip enrichment, International Journal of Fracture 204: 55-78.
• Piegl, L., Tiller, W. (1995). The NURBS book, Springer (Berlin).
• Roger, D.F. (2001). An introduction to NURBS with historical perspective, Morgan Kaufmann Publishers (Burlington).
• Scott, M.A., Simpson, R.N, Evans, J.A., Lipton, S., Bordas, S.P.A., Hugues, T.J.R., Sederberg, T.W. (2013). Isogeometric boundary element analysis using unstructured T-splines, Computer Methods in Applied Mechanics and Engineering 254: 197-221.
• Sederberg, T.W., Parry, S.R. (1986). Free-form deformation of solid geometric models, ACM SIGGRAPH Computer Graphics 20: 151-160.
• Sederberg, T.W., Zheng, J., Bakenov, A., Nasri, A. (2003). T-splines and T-NURCS, ACM Transaction on Graphics 22: 477-484.
• Simpson, R.N, Bordas, S.P.A., Trevelyan, J., Rabczuk, T. (2012). A Two-dimensional isogeometric boundary element method for elastostatic analysis, Computer Methods in Applied Mechanics and Engineering 209: 87-100.
• Simpson, R.N., Bordas, S.P.A., Haojie, L., Trevelyan, J. (2013). An isogeometric boundary element method for elastostatic analysis: 2D implementation aspects, Computers & Structures 118: 2-12.
• Strouboulis, T., Copps, K., Babuška, I. (2001). The generalized finite element method, Computer Methods in Applied Mechanics and Engineering 190: 4081-4193.
• Zienkiewicz, O.C., Taylor, R.L. (1967). The finite element method: solid mechanics, Butterworth-heinemann (Oxônia).

Edited by

Editor: Marco L. Bittencourt and Josué Labaki

Publication Dates

• Publication in this collection
  04 Aug 2023
• Date of issue
  2023

History

• Received
  15 Mar 2023
• Reviewed
  04 May 2023
• Accepted
  26 June 2023
• Published
  30 June 2023