Isogeometric analysis of 3D solids in boundary representation for problems in nonlinear solid mechanics and structural dynamics

This contribution presents an isogeometric approach in three dimensions for nonlinear problems in solid mechanics and structural dynamics. The proposed approach aligns with the boundary representation modeling technique in CAD. Isogeometric analysis is combined with the parameterization of the scaled boundary finite element method. In this way, the boundary description of 3D solids in CAD can be directly utilized for the analysis in an isogeometric framework. The approximation of the solution on the boundary is based on bi‐variate NURBS. We also approximate the interior of the solid with uni‐variate B‐splines to facilitate the solution of nonlinear problems. The nonlinear case is extended to three dimensions in this present work. The method is further extended to dynamic problems. The mass and damping matrices are derived using the same basis functions. The solution of the entire solid is obtained with the Galerkin method. We study several nonlinear and dynamic problems with simple geometries and arbitrary number of boundaries. Moreover, a fiber reinforced composite serves as demonstration for complex shapes.

a collocation approach 9 or by using the deformed mesh of a simple geometry. 10 Moreover, NURBS-enhanced finite elements have been proposed to exploit the exact boundary representation. 11,12 This preserves the smoothness of the boundary, whereas a standard polynomial basis is used to approximate the solution. Alternatively, a hybrid representation of tri-variate NURBS on the boundary and unstructured tetrahedral elements in the interior of the structure can be constructed. 13 Adopting higher-order Bézier tetrahedra in the interior improves the continuity of the volume representation. 14 Another approach is based on Boolean operations to obtain volumetric T-splines and thereafter tri-variate Bézier elements. 15 One can also obtain suitable tri-variate models by decomposing the boundary surfaces into cuboids that preserve the topology of the boundary and provide a tri-variate spline representation for isogeometric analysis. 16,17 Moving from the geometry to analysis, researchers have focused on analysis procedures in boundary representation. 18 This is sufficient in case the geometry is represented by a surface, as for example in shell analysis. 19 However, the boundary description is not sufficient to analyze solids. To this end, the boundary representation of CAD can be exploited in combination with a NURBS approximation in the interior of the solid. 20 This approach combines isogeometric analysis with the scaled boundary finite element method (SBFEM). The latter constitutes a semianalytical method employing polar coordinates for the parameterization in two dimensions, which was introduced for linear elasticity problems 21 and further extended to structural dynamics. 22 Its relation to the finite element method has been explored based on a virtual work derivation. 23 Scaled boundary polygon formulations have been derived to use SBFEM for elasto-plastic problems 24 and linear elastodynamics. 25 Furthermore, special shape functions derived from the analytical solution of SBFEM can be used for nonlinear problems. 26 Most recently, SBFEM has been enhanced to solve elasto-plasticity in three dimensions as well as dynamic problems. [27][28][29][30] The main feature of these studies is that they employ scaled boundary polygonal or polyhedral elements combined with quadtree decomposition for automatic mesh generation. This approach is based on a Lagrange approximation on the boundary, whereas other types of shape functions are in principle possible. [31][32][33] At this point, it should be stressed that in particular vibration analysis can benefit from the isogeometric approach. 34 Non-uniform rational B-splines (NURBS) offer higher smoothness, which can be combined with a higher-order refinement procedure (k-refinement) for superior accuracy in structural vibration problems. 35 Their advantage over higher-order Lagrange interpolations has been shown for example, the frequencies of higher modes. 36 SBFEM has been combined with B-spline approximations in the context of IGA. 37,38 However, these studies focus on linear problems.

Contributions and outline
This contribution deals with an isogeometric framework for the analysis of 3D solids in boundary representation. In this context, a nonlinear formulation is derived in three dimensions. Unlike existing 3D approaches for elastostatics, 20 the derived formulation exploits a general solution form based on the principle of virtual work instead of the scaled boundary finite element equation. This enables the solution of nonlinear problems as shown for 2D static problems in our previous works. 39,40 Here, we extend the nonlinear case to three dimensions and dynamic problems for the first time. To this end, the mass and damping matrices are derived using a consistent approach. The formulation exploits the parameterization of SBFEM in order to fit the boundary representation modeling technique in CAD. We adopt the bi-variate NURBS basis functions from CAD to approximate the displacements on the boundary surfaces. The approximation in the interior of the solid is based on uni-variate B-splines. In other words, the basis functions are decomposed into two-dimensional NURBS on the boundary that are a function of the B-splines in the interior of the solid. This differs from IGA, where a tri-variate tensor-product structure is used to approximate the displacement solution of the entire solid with NURBS. The weak form of equilibrium is employed to obtain the solution of the volumetric domain. The main features of the proposed approach are summarized: • Direct analysis of volumes based on the original surface representation from CAD.
• B-spline approximation in the interior of the solid allows for physical and geometrical nonlinearities.
• Enhanced to account for dynamic problems.
• Feasible to model arbitrary volumetric domains based only on the boundary description. Flexibility in meshing without restrictions of a single tri-variate tensor-product structure for the entire solid.
In view of the above contributions, this article is structured as follows: Section 2 presents the parameterization for the 3D case and Section 3 the approximation with NURBS. Section 4 provides a brief overview of the equilibrium equations and the discretized weak form. Section 5 is concerned with several numerical problems in structural vibrations and nonlinear mechanics. For comparison, we employ computations with isogeometric analysis and finite elements. Section 6 summarizes the main conclusions that can be drawn from this study and provides an outlook for future work.

PARAMETERIZATION OF 3D SOLIDS
Hereafter, the parameterization is presented considering nonlinear kinematics for the three-dimensional case. The reference configuration of the domain is denoted as Ω 0 . Based on previous work in two dimensions, the parameterization is inspired by the scaled boundary method. 39,40 The geometry is described by using a radial scaling parameter 0 ≤ ≤ 1 to define the interior of the solid. The boundary of the solid consists of two-dimensional surfaces that are scaled with respect to the scaling center C as depicted in Figure 1. In relation to the scaling center, the domain is partitioned into sections Ω 0,s , which are bounded by the boundary surfaces Ω 0,s . These are parameterized in the circumferential direction with the boundary parameters and , where 0 ≤ ≤ 1 and 0 ≤ ≤ 1.
The coordinates of the scaling center C are contained in the position vectorX = (X 1 ,X 2 ,X 3 ). The coordinates of a point on the boundary surfaces and in the interior of the solid are contained in the position vectorsX = (X 1 ,X 2 ,X 3 ) and X = (X 1 , X 2 , X 3 ), respectively. Any point on the boundary or in the domain can then be expressed as where N b ( , ) contains the NURBS basis functions that describe the boundary surfaces. Note that the exact geometry of the boundary surfaces designed in CAD is exploited within the presented approach. The vector X s contains the coordinates of the control points on the boundary. The position vectors on the section Ω 0,s are depicted in Figure 1.
Considering Equation (1), the Jacobian J of the domain can be expressed in matrix notation as .
(2) This results in a multiplicative decomposition of the determinant det J = 2 det J = 2 J. The transformation of a volume element from the physical space to the parameter space reads The inverse of the Jacobian is necessary to obtain the partial derivatives in Cartesian coordinates. With this in mind, the differential operator reads in matrix notation Employing the scaled boundary parameterization, the differential operator can be rewritten in tensor notation as where the matrices b 1 , b 2 , and b 3 are only a function of the boundary parameters and (see Appendix). For geometrically linear theory, the position vectors are replaced by their counterparts in the current configuration and the linear differential operator can be derived as shown in previous work. 39 One limitation is that the scaling requirement needs to be fulfilled in order to employ the presented parameterization, in other words the scaling center should be visible from the entire boundary. In case of complex geometries, subdivision into star-shaped subdomains is required. 38,39 Some algorithms for the subdivision are assessed in the literature. 41

APPROXIMATION WITH NURBS
In order to align with the idea of isogeometric analysis, the displacements u at the boundary surfaces are approximated with the same NURBS basis functions that describe the geometry. These are directly adopted from the geometry model, which is designed with the boundary representation modeling technique in CAD. Therefore, the approximation of the displacements can be expressed as where U s,I contains the displacements of the control point with index I and n bc denotes the total number of control points on the boundary surface Ω 0,s . The NURBS basis functions employed to describe the boundary surfaces are termed as boundary NURBS and are independent of . They can be computed as described in the literature. 42 Note that these are the same basis functions as in Equation (1). The associated control points are termed as boundary control points. It is remarked that the index I represents a global node number. After rearranging the displacement vectors U s,I of all control points in the vector U s , Equation (6) can be rewritten in matrix notation as The dimension of the displacement vector U s,I for I = 1, .., n bc is three, as three nodal degrees of freedom are assumed for each boundary control point. In scaling direction, B-splines are employed for the approximation of the displacements to allow for nonlinear problems. 39,40 These are denoted with R r j , where r is the polynomial degree in scaling direction. Note that the scaling lines constitute straight lines and the weighting factors are equal along each line. Approximating with B-splines in scaling direction, the displacement vector U s,I becomes only a function of the radial scaling parameter and reads Here, n cp is the number of control points along each radial scaling line and N s,I the uni-variate B-spline basis functions approximating a radial scaling line. It is remarked that the displacement vector U r,I for I = 1, .., n bc refers to the control points along a radial scaling line and its dimension is 3 ⋅ n cp , as three nodal degrees of freedom are attributed to each control point in scaling direction. Equation (8) can be rewritten in matrix notation as The displacement vectors U r,I associated with each scaling line can be rearranged into the vector U containing the displacements of all control points in the scaling direction of Ω 0,s . The basis functions N s in scaling direction are then rearranged in the matrix The dimension of N s is (3 ⋅ n bc , 3 ⋅ n bs ), where n bs = n cp ⋅ n bc is the total number of control points in the section. Inserting Equation (10) into Equation (6) leads to a multiplicative decomposition of the displacements as well as the virtual displacements, where it holds The vectors U and U contain all nodal degrees of freedom in the scaling direction of a section Ω 0,s . Moreover, N( , , ) can be written in matrix notation as with dimension (3, 3 ⋅ n bs ).  be found in previous studies. 43 These indicate that the singularity leads to condition numbers of higher order regarding the stiffness matrix, which grow fast with mesh refinement. Note that we constrain the solution at the scaling center by linking the respective degrees of freedom, which reduces this effect. The domain Ω 0 consists of several sections and is a multipatch structure. Instead of employing a tri-variate tensor-product for the approximation of each section, the displacement response is derived here by decomposing into a bi-variate tensor-product on the boundary that is a function of the uni-variate B-splines in the interior of the domain. In this way, the boundary representation from CAD is sufficient to define the volumetric domain. Here, adjacent boundary surfaces are assumed to have the same number of control points and polynomial degree. Furthermore, all radial scaling lines are assumed to be identical. In general, boundary surfaces might not be connected to each other in the geometry model. Also, an efficient discretization of the domain might necessitate local refinement with nonmatching B-spline approximations in scaling direction. The presented approach is, in principle, flexible and can be used for nonconforming discretizations by adopting an appropriate coupling method. 44,45

EQUILIBRIUM EQUATIONS
In this section, the definitions of the stress, strain and material tensors of the three-dimensional continuum theory are adopted. For each time instant t and material point X = (X 1 , X 2 , X 3 ), the unknown time-dependent displacement field is u = u(X, t). The weak form of the initial boundary value problem is considered as where u ∈ {H 1 | u = 0 on u Ω 0 } is the virtual displacement and P is the first Piola-Kirchhoff stress tensor in the reference configuration of the domain Ω 0 . The density is denoted by , the body force by b 0 and the acceleration vector byü. The Dirichlet and Neuman boundary conditions are given as PN = t 0 on t Ω 0 and u = u on u Ω 0 respectively. Note that N is the normal vector perpendicular to the boundary Ω 0 . The initial boundary conditions are expressed by u(t = 0) = u 0 andu(t = 0) =u 0 for t ∈ (0, T), where T denotes the end of the dynamic process. Equation (13) is rewritten by applying integration by parts as where S is the second Piola-Kirchhoff stress and E denotes the Green-Lagrange strain. These can be represented in . Considering the components of the deformation gradient F ik = (X i +u i ) X k for i, j, k = 1, … , 3, the first variation of the Green-Lagrange strain in index notation reads E ij = 1 2 ( F ki F kj + F ki F kj ), where the summation convention over repeated indices holds. With respect to Voigt notation and symmetry of E, it is rewritten as Here,F is a 6 × 9 matrix containing the components of the deformation gradient. F is the gradient of the virtual displacements u and is represented as a vector of dimension nine. According to the interpolation of the virtual displacements in Equation (11) and applying the differential operator from Equation (5), the interpolation of F reads Considering Equation (3), the discretized weak form of equilibrium is written for one section Ω 0,s as F int s , M s , F ext s denote the internal force vector, the mass matrix and the external force vector for one section Ω 0,s , respectively. The same NURBS basis functions are used to approximate the accelerations and virtual displacements, thus it holds N = N b ( , )N s ( ) according to Equations (11) and (12). It is remarked that a consistent mass matrix is adopted here, as the lumping approach has performed poorly in eigenvalue analysis of higher order NURBS elements. 3 For the numerical integration we employ standard, full Gauss quadrature with (p + 1) × (q + 1) × (r + 1) quadrature points in each parametric direction. It is remarked that the study of quadrature rules tailored to the proposed approach is out of the scope of this paper and is an aspect of future work. Here, a constant damping is considered, which results in the damping matrix on section level d denotes the modal damping ratio of the first mode and 1 the first natural frequency. Assembling over all elements leads to the residuum vector The vectors V,V,V represent the total number of degrees of freedom after assembly over all sections for the displacement, velocity and acceleration, respectively. For the discretization in time the Newmark method 46 is employed. Let Δt = t n+1 − t n be the time step size and n the discrete step number. The current displacement and velocity are given as For an implicit time integration the parameters are defined as = 1∕2 and = 1∕4. In the discrete time setting the residuum vector from Equation (19) is rewritten as In case of geometric or material nonlinearity, the internal force vector F int n+1 is a nonlinear function of the current displacement vector V n+1 . The Newton-Raphson scheme is employed to solve Equation (21) at time t n+1 . It necessitates a linearization of the internal force vector F int n+1 . On section level the Gâteaux derivative in the direction of the displacement increment ΔU = U k+1 n+1 − U k n+1 is employed, where the superscript k denotes the iteration step. Since S andF depend on the actual displacement-state, the chain rule is applied For the sake of a compact notation, we will skip the index n + 1 and remark that all quantities in the latter equation are evaluated at the actual time step t n+1 . The quantity ΔE is a function of the displacement increment ΔU. Considering Equations (15), (16) and substituting the first variation by the second one it follows The quantity ΔF in Equation (22) is a 6 × 9 matrix containing the components of the gradient of the displacement =ŜΔF. (24) Replacing ΔF = BΔU the linearization of the internal force vector on section level reads where C T = S E denotes the consistent tangent modulus and K s the tangential stiffness matrix for the section Ω 0,s . It is remarked that the matricesF andŜ contain several zero entries. From a programming viewpoint, the use of sparse matrix form can be used for savings in memory storage.

NUMERICAL EXAMPLES
In this section, 3D benchmarks demonstrate the ability of the proposed approach to analyze volumetric domains designed with the boundary representation technique. We focus on problems with physical and geometrical nonlinearities as well as dynamic analysis. The first example is a simple geometry of a cantilever beam with a known analytical solution for the linear case to demonstrate the application to structural dynamics problems. For better illustration, the eigenfrequency and vibration analysis is evaluated by comparison to IGA and FEM. We employ different discretization schemes and polynomial degrees of NURBS in order to evaluate the accuracy of the proposed method. The second example is a hollow thick cylinder subjected to uniform internal pressure. Elasto-plastic material behavior is considered and the convergence for different discretization schemes is evaluated. Moreover, a comparison with IGA and FEM is provided in order to evaluate the overall deformation and stress behavior for both static and dynamic analysis. The third example concerns a block and L-shaped geometry under compression, which aim to verify the proposed approach for geometrically nonlinear problems. To this end, the results are compared with FEM. The fourth example is a free-form solid, which is chosen to illustrate the capability of the presented approach to handle solids with arbitrary geometry. Nonlinear computations are performed with the Neo-Hookean material law and compared to FEM. The last example is a fiber reinforced composite, which illustrates the application to complex, nonstar-shaped domains. All models are developed using Rhinoceros 3D. Finite element computations are denoted with FEM 1st or FEM 2nd in the legend, depending on the order of the shape functions. Standard isogeometric analysis is denoted as IGA in this section. The presented approach has been implemented into a customized isogeometric version of the finite element analysis program (FEAP). 47 All computations are performed in FEAP unless stated otherwise. Table 1 summarizes the declarations for the description of the computational models. For the computations with the presented formulation, the polynomial degree is the same in all parametric directions, thus p = p b = p c unless stated otherwise. The same holds also for the number of elements n = n b = n c .

Cantilever beam
The first example to be considered here is a cantilever beam. The aim of this benchmark is to validate the proposed approach for dynamic problems. In particular, the dynamic behavior of simple geometries is investigated with the proposed discretization scheme ( Figure 3). For simplicity, a one-dimensional problem is considered where the cantilever is discretized as a solid. The beam is modeled only by its boundary surfaces. The geometry and boundary conditions are depicted in Figure 3. The dimensions are the length L = 20 m, the width B = 2 m and the height H = 1 m. All displacements are fixed at the left side as shown in Figure 3. Plane strain conditions are assumed and displacements are fixed in z-direction as shown in Figure 3. An analytical solution exists for this problem assuming linear theory. For time-dependent analysis, the beam is subjected to a triangular force excitation at the free end as shown in Figure 4, which induces an initial displacement and velocity in y-direction. Free vibration of the system occurs in y-direction after t = 1 s. It is remarked that the resultant force of the applied line load is depicted in Figure 4. Linear kinematics and elastic material behavior are assumed here. The material parameters are the Young's modulus E = 10 9 N/m 2 , the Poisson's ratio = 0.0 and the density = 1000 kg/m 3 . The domain is bounded by six surfaces, which results in six sections with respect to the scaling center. For the computations, the polynomial degree of NURBS basis functions is chosen as p = 2, 3, 4. The order elevation is handled by p-refinement. The number of elements in all directions in initially set to n = 2 and gradually refined to n = 4, 8, 10, 16 with k-refinement. First, eigenfrequency analysis is carried out. To study the convergence behavior, a reference solution is computed with isogeometric analysis. For the purpose of assessing the quality of the computed reference solution, we consider the analytical solution of an Euler-Bernoulli beam. An analytical solution for the natural frequencies i = 1, 2, 3 … of the employed system is given in the literature 48  where i = i L represents the solution of the characteristic frequency equation Alternative computations with standard isogeometric analysis for p = 4 and n = 64 result in 1 = 5.0419738 rad/s and 2 = 30.447204 rad/s. Note that a small approximation error of 0.8% is observed in comparison to the analytical solution. The error is quite low and the computed reference solution is chosen for better comparison with the numerical results of the proposed approach. The low natural frequencies are commonly most predominant for the dynamic excitation of the system. Thus, the relative error of the first and second natural frequency is evaluated here for the proposed approach with different polynomial degrees. In Figure 5, the error of the approximated natural frequencies h i obtained by the computations is calculated as ( h i − i )∕ i . The solution converges for all polynomial degrees and the error decreases with mesh refinement. An error magnitude of 10 −6 is reached for 1 and 10 −5 for 2 with p = 4. For 1 , the convergence rate is almost the same for all polynomial degrees, whereas p = 3 converges faster in the fine range of discretizations. The situation for 2 is similar, whereas the convergence rate is higher for p = 4. The results are compared with finite element computations with quadratic solid elements. The number of elements is set to n = 4, 8, 16, 32 per side for FEM. Note that a lumped mass approach is adopted, whereas the proposed formulation is based on a consistent mass approach. It is observed, that the proposed approach for p = 2 is less accurate than quadratic FEM. This might be due to the differences in the structure of the employed discretization scheme and the choice of discretization parameters in scaling direction. Superior accuracy of the eigenfrequencies is obtained in comparison to quadratic elements by increasing the polynomial degree to p = 3 and p = 4. Note that the convergence rate is p-1, which is measured with respect to the degrees of freedom. We remark that the error of eigenfrequencies for the Euler-Bernoulli beam converges as h 2(p−1) under optimal conditions in IGA. 34 The convergence of higher polynomial degrees may be affected by overintegration. Optimal quadrature rules can be used instead of standard Gauss quadrature to improve this. 49 Furthermore, dynamic analysis is carried out. The purpose here is to compare the transient response to standard IGA and FEM. First, we study the influence of h-refinement to the undamped free vibration. To this end, the displacement u y at point A is monitored (see Figure 3). The time step is set to Δt = 0.005 s for the analysis. The values of Newmark parameters are chosen as = 0.25 and = 0.5, therefore the employed implicit time integration scheme is unconditionally stable with no numerical damping. free vibration is distinguished in the curve. Thus, the natural period of vibration is T = 2 ∕ 1 = 1.25 s, which matches the previously employed analytical solution. Figure 6 (left) shows that by increasing the number of elements the results become nearly identical with the reference solution. Moreover, the free vibration of the system with damping is considered and the influence of p-refinement is investigated. The damping ratio is chosen as d = 0.05. Figure 6 (right) depicts the displacement-time curve for different polynomial degrees with n = 2. As expected, the displacement amplitude decreases with each vibration cycle. The exponential decay can be described by the envelope curve ± exp − d 1 t , where it holds 50 Note that damping decreases the natural frequency to D = 1 The period increases marginally with the chosen damping ratio. This can be directly observed in the curves of Figure 6 (right), where one cycle of free vibration can be distinguished at nearly the same time instant as in Figure 6 (left). Note that the imposed load is applied slowly and acts as a static load until t = 1.0, which causes an initial displacement in the direction of the first eigenmode. For the response of the cantilever, the first eigenmode is predominant. Therefore, only the fundamental mode is considered in Equation 28. Note that the solution with p = 2 slightly underestimates the damping, whereas the solution with p = 3 overestimates the damping. By increasing the polynomial degree to p = 4, more accurate results are obtained and the time-displacement curve touches the envelope curve around the peak values. It should be remarked, that a total of only 48 quartic elements are required to accurately capture the damped vibration. To sum up, the results validate the proposed approach for the simplified 1D dynamic case. Based on this problem, the results indicate that the accuracy can be increased by employing higher polynomial degrees and refinement strategies such as k-refinement.

Hollow cylinder subjected to uniform internal pressure
The aim of the second example is to assess the proposed formulation for elasto-plastic analysis of static and dynamic problems. For this purpose, a curved geometry of a hollow thick cylinder subjected to uniform internal pressure is studied. Note that the setup of this benchmark is adopted from the literature. 1 The capability of the proposed formulation is evaluated here for the 3D case by comparison to IGA, FEM and an analytical solution. 51 Moreover the influence of h-refinement and order elevation in scaling direction is investigated. The geometry and boundary conditions are depicted in Figure 7 (left). Small strain theory is assumed. Although the employed system constitutes an academic example, all dimensions can be assumed in meter for simplicity. The dimensions of the cylinder are R min = 5 and R max = 10 m for the inner and outer radii and the height is given by h = 5 m. An elastic-perfectly plastic material without hardening is used, where the elasticity modulus E = 100 N/m 2 , the Poisson's ratio = 0.0, the yield stress y = 50 N/m 2 and the density = 1 kg/m 3 are chosen. For the static case the magnitude of the inner pressure is q = 10 N/m 2 . The symmetry of the model is exploited such that only one quarter of the system is modeled and symmetry boundary conditions are employed. Plane strain conditions are assumed in z-direction for kinematic stability, therefore the test can be simplified to a two-dimensional problem. For the purpose of evaluating the capability of the proposed approach, the system is modeled here in three dimensions. The scaling center C is chosen at x = y = 6.9 m, which corresponds to a radius R c = 9.758 m. Note that the average of all controls on the boundary would lead to a scaling center outside of the domain. With respect to the scaling center, the domain is partitioned into six sections bounded by the surface representation designed in CAD (see also Figure 7 (left)). Note that only the depicted boundary surfaces are modeled in Rhino3D and adopted for the analysis model.
The boundary surfaces and the interior of each section are described by NURBS basis functions with polynomial degrees p = 2, 3 for the analysis. For conformity of adjacent patches, the number of elements in all three directions of each section is initially chosen as n = 2 and then refined to n = 4, 8, 16 with k-refinement. Elasto-plastic analysis is performed for the static case and the inner pressure is applied incrementally in a total number of 12 load steps to obtain a load factor = 3.0. The relative error of the displacement u y,A at point A in the elastic regime and u y,B at point B in the plastic regime are monitored (see Figure 7 (right)). A two-dimensional isogeometric computation using p = 6 and 4096 elements is used as a reference, which yields u y,A = 2.20445 m and u y,B = 2.89313 m. The analytical solution 51 also assumes plane strain conditions. Therefore, the strain in thickness direction is zero. It is assumed that elastic and plastic strain both vanish separately, where actually only the sum of the two should vanish. Because of this assumption, the analytical solution entails an approximation and is not used here for comparison of the displacements.
To investigate the influence of order elevation and k-refinement, the relative error of displacements u y,A and u y,B is compared for the same polynomial degrees on the boundary and in scaling direction p = p b = p c over the ratio of element number n c ∕n b , see Figure 8. The number of elements in scaling direction n c is increased from 0.5n b to 2n b . Similar as in Reference 52 a kink of the convergence lines can be observed for n b = n c on some of the curves, which leads to the similar recommendation of using n c > n b for reliable results. Generally speaking it can be said that for increasing the number of elements on the boundary n b , as well as for increasing the number of elements in scaling direction n c the accuracy improves. Order elevation, while keeping p = p c = p b , generally improves the accuracy, however the influence in the plastic and elastic regime is different. While in the elastic regime (point A) the influence is greater for smaller number of boundary control points n b , in the plastic regime (point B) the order elevation is more significant for larger number of n b . Because h-refinement in scaling direction might not be optimal in terms of efficiency, 40 the polynomial degree p c in scaling direction is increased and convergence is observed, see Figure 9. The number of control points on the boundary are refined to n b = 2, 4, 8, 16. Generally, the order elevation in scaling direction leads to more efficient results than h-refinement. However, it is difficult to draw generalized conclusions because many factors influence the convergence behavior. For more detailed information about the influence of the polynomial degrees p c , p b and the number of control points on the convergence behavior, one can refer to Reference 52. Moreover, Figure 10 depicts the comparison of the load deformation curve for u y,A between two-dimensional quadratic IGA and the proposed approach. The mesh with p = 2 and n = 8 is employed for the comparison. The resulting proposed mesh with a total of 3072 elements is depicted in Figure 7 (right). In addition, the minimum and maximum pressures, for which parts of the cylinder deform plastically are depicted in Figure 10. The range of pressures is defined analytically 51 as As mentioned previously, this solution is only approximate. However, the nonlinear part of the curve fits well into these bounds. The results of the proposed formulation are nearly identical with IGA. This verifies the obtained results for elasto-plastic material behavior.  Figure 11 depicts the distribution of the radial displacements. A very good agreement with FEM can be observed. Moreover, the von Mises stress distribution is depicted in Figure 12. It can be observed that the red yield zone propagates from the inner surface to approximately r = 7.5, which agrees well with the finite element solution. The proposed mesh consists of 6000 control points and requires overall 12,779 equations for the solution. Finally, the time-dependent F I G U R E 13 Cylinder subjected to uniform internal pressure: displacement-time curve at point B and comparison of proposed formulation with quadratic FEM deformation behavior including material nonlinearity is investigated. The elastic-perfectly plastic material behavior with the same material parameters as for the static case is adopted. For the dynamic case, the inner pressure is applied as a harmonic function of time q = 30 sin(2 t). Note that the harmonic load starts as zero at t = 0 and continues periodically until t = 3 sec with an amplitude of A = 30. This corresponds to a load factor of = 3.0 in order to assess the behavior in the nonlinear range (see static case, Figure 10) . The time step is chosen as Δt = T∕2000 = 0.005 s. It is remarked that the employed Newmark scheme is unconditionally stable with parameters = 0.25 and = 0.5. The displacement at point B (see Figure 7 (right)) is monitored over time. Figure 13 depicts the displacement-time curve for p = 2 and different levels of h-refinement. The results are compared to finite elements computations with quadratic solid elements in ANSYS Workbench. 53 Due to the employed harmonic pressure function, the displacements follow a similar path over time. It can be observed, that the displacement amplitude and thus the accuracy depends on the employed level of refinement. For the proposed discretization requiring 2139 equations, the results are almost identical with FEM using 3852 equations. It can be concluded that for the simplified 2D problem the obtained results with the proposed approach are as good as FEM for both static and dynamic computations with material nonlinearity. Regarding the choice of discretization parameters, some tendencies become visible by investigating the influence of h-refinement and order elevation. As the convergence behavior is influenced by many factors, the reader is further referred to the literature for a thorough investigation. 52

Compression test
The third example is concerned with a compression test. For this purpose, a 3D block is considered first. This test aims to validate the proposed formulation for geometrically nonlinear analysis with large deformations and strains. The second geometry is a nonconvex domain of an L-profile, thus it differs from the previous examples. For simplicity, the static case is considered.

Block
The first geometry is concerned with a 3D block. The geometry and boundary conditions are adopted from the literature, where they served as a benchmark for enhanced element formulations under severe locking. 54,55 Here, strongly nonlinear behavior is examined in the scope of the presented displacement-based formulation. To this end, a comparison with displacement-based FEM is employed. The geometry and boundary conditions of the system are illustrated in Figure 14 (left). Symmetry is exploited here, therefore, only one quarter of the system is modeled. The length and width of the block are set to L = B = 1 mm and the height to H = 1 mm. The Dirichlet boundary conditions include the symmetry conditions on the side surfaces, fixed vertical displacements on the bottom and fixed horizontal displacements on the top as shown in Figure 14 (left). The 3D block is partially loaded on the top with a vertical surface load q = 4 MPa, where denotes the load factor (see Figure 14 (left)). The load is applied incrementally in 20, 40, 60, 80 load steps in order to investigate the convergence behavior for various load levels causing moderate to large deformations. It is remarked that dead load is considered during the structural deformation, thus the load is acting on the reference configuration. The material properties of the Neo-Hooke material model are the shear modulus = 80.194 MPa and Lamé parameter Λ = 120.291 Pa. This corresponds to a Poisson's ratio of = 0.3. Figure 14 (middle) depicts an exemplary distorted mesh at the final deformed configuration for the load factor = 60. The geometry is modeled with B-Splines. The scaling center is placed at the average of all control points on the boundary surfaces and partitions the domain into six sections. Each section can be regarded as a volume patch. The polynomial degree is chosen as p = 2. The number of elements is initially set to n = 2 and increases to 4, 8, 16 with k-refinement. It is remarked that the polynomial degree and number of elements is equal in all parametric directions. A representative mesh with n = 8 can be seen in Figure 14 (middle). For comparison, solid finite elements with quadratic shape functions are employed. The finite element mesh is discretized with n = 2, 4, 8, 16 elements per parametric direction for the comparison. Geometrically nonlinear analysis is performed. The percentage compression level is monitored at point A (see Figure 14 (left)).
The results are depicted in Figure 15. It can be observed, that apart from the coarsest discretization for increased load levels, the results of the proposed formulation are similar with FEM. Both solutions converge with mesh refinement. Note that due to the simple geometry with straight edges, FEM already provides a good approximation. Nevertheless, we focus here on the overall capability to handle large deformations and strains. Moreover, Figure 14 (right) illustrates the von Mises stress distribution on the distorted mesh, which is computed with the proposed formulation for the load factor = 60 and a mesh with n = 8. Note that the illustration refers to the stress distribution in the reference configuration. Table 2 shows the convergence of the Newton-Raphson scheme for the proposed formulation and FEM. For this purpose, an exemplary mesh is chosen as the finest mesh with 16 elements per parametric direction for FEM and eight elements per direction of each section for the proposed formulation. This results in a slightly finer finite element mesh. Note that for coarser meshes the convergence of the Newton-Raphson scheme is similar. It can be observed that both solutions converge nearly quadratically with the same number of iterations, which implies a comparable computational cost. Overall, accurate results could be obtained with the proposed formulation for nonlinear behavior distinguished by large deformations and strains as well as various load levels.

L-shaped geometry
The second geometry is an L-shaped profile. The geometry is depicted in Figure 16 (left). The dimensions of the system are L 1 = B 1 = 9 m, L 2 = B 2 = 1 m and H = 1 m as depicted in Figure 16 (left). All displacements are fixed on the left square surface. Displacement control is employed, where the displacement u y = 10 m is applied incrementally on the top square surface in 20 load steps. The shear modulus for the Neo-Hooke model is chosen as = 160 MPa and the Poisson's ratio is = 0.0. It is remarked that this results in a 3D stress state with vertical displacements significantly smaller than the in-plane displacements. Therefore, the stresses in the vertical direction could be neglected and the system would simplify to a 2D problem. Here, we model the problem in 3D to investigate the proposed approach for geometries modeled by multiple subdomains. It is remarked that although the geometry is star-shaped, the subdivision of the domain leads to better visibility of the scaling center from the boundary. Note that in the scaled boundary approach, low visibility tends to affect adversely the numerical results as shown in previous studies. 40,56 The geometry is divided into three subdomains in order to ensure a conforming discretization between adjacent surfaces. The borders of each star-shaped subdomain are visible in Figure 16 (left). For each subdomain, a scaling center is defined as the average of all control points on the boundary. For the discretization the polynomial degree is initially set to p = 2 for all parametric directions and then elevated to p = 3. The number of elements in scaling direction and on the boundary of the square surfaces is initially set to n = 1 and then gradually refined to n = 2, 3, 4 with k-refinement. The boundary along the long edges is discretized with n b = 8n = 8, 16, 24, 32 elements accordingly. Nonlinear analysis is carried out and the displacement is applied in 20 increments. Figure 16 Figure 16 (left)). For comparison, finite element computations with linear and quadratic shape functions are employed. The number of elements for the finite element mesh is refined as n = 1, 2, 3, 4, 5 per parametric direction. Figure 17 depicts the convergence of the displacement u y,A at point A for different polynomial degrees. The results are compared to finite element computations. It is observed that the proposed approach with quadratic NURBS converges slower than finite elements for coarse discretizations. It is noted that finite elements appear to converge from above, whereas the proposed approach converges from below. By elevating the polynomial degree to p = 3, nearly converged results are obtained also in the coarse range of discretizations. All computations converge to identical results with mesh refinement. Moreover, Figure 18 (left) shows the load deformation curve for the displacement u x,B at point B. Note that the load factor refers to the incrementally imposed displacement u y,B . The finest discretization at the final load step is chosen for all computations. It can be observed that the results are identical and the curves are overlapping. The equivalent von Mises stresses at point B are compared in Figure 18 (right) over the prescribed displacement u y,B for the finest discretization. For better comparison, isogeometric computations are carried out for degrees p = 2, 3 and the same boundary discretization as the proposed approach. It can be observed that the proposed approach slightly underestimates the stress, whereas increasing the polynomial degree to p = 3 leads to better agreement with quadratic IGA. It can be concluded that for both displacements and stresses, increasing the polynomial degree leads to more accurate results.

Free-form solid
The fourth example concerns a solid with free-form geometry. This benchmark aims to illustrate the capability of the presented formulation to analyze solids with arbitrary shape and curved surfaces. A common application of such geometries is, for example, suspension elements for industrial machines. This geometry is likewise inspired by rubber springs typically employed for vibration isolation. The geometry is modeled in CAD with cubic NURBS by using only two symmetric boundary surfaces. The boundary surfaces are then directly adopted from the CAD model and used for the analysis. The system is depicted in Figure 19 (left). A surface traction with a magnitude of q = 4 N/m 2 is applied on the entire area of both boundary surfaces as depicted in Figure 19 (left). The Dirichlet boundary conditions are applied by fixing the displacements at the symmetry axis. For the analysis, the Neo-Hookean material model is employed with shear modulus = 38.4615 Pa and Lamé parameter Λ = 57.6923 Pa. This corresponds to a Poisson's ratio of = 0.3. The scaling center coincides with the centroid and is placed at the symmetry axis. This partitions the solid into two sections. For the analysis, the model is discretized with a total of 2000 elements. The polynomial degree in scaling direction is first set to p c = 1 and then refined to p c = 2, 3 with order elevation, whereas the total number of elements is fixed. A representative deformed mesh is depicted in Figure 19 (middle). Note that the finest mesh consists of 4394 control points and requires 9576 equations for the solution. We perform static analysis allowing geometrical nonlinearities and apply the surface load incrementally in 10 load steps. To compare the solution, finite element computations are performed in ANSYS Workbench. 53 We use quadratic solid elements and refine by decreasing the element size (h-refinement). The finest finite element mesh consists of 2857 quadratic solid elements, 8828 nodes and requires 26,352 equations for the solution (see Figure 19 (right)). Figure 20 shows the comparison between the maximum equivalent von Mises stresses obtained with the proposed meshes and FEM.
By increasing the polynomial degree in scaling direction, the stresses converge faster with the proposed approach requiring less degrees of freedom than FEM. For a qualitative comparison, Figures 21 and 22 illustrate the horizontal displacements (y-direction) and von Mises stresses of the final deformed configuration for the finest mesh with both solutions. It is noted that the Cauchy stress distribution is considered here, which is computed in the current configuration. The comparison indicates that the proposed approach is in very good qualitative agreement with the finite Note that the proposed mesh with p = 3 leads to a better approximation of the curved surfaces than quadratic FEM requiring less degrees of freedom. The presented approach makes direct use of the surface representation and is thus perfectly suitable for geometries with curved surfaces and arbitrary number of boundaries designed with the boundary representation technique in CAD.

Fiber reinforced composite
The last example serves as a demonstration for complex geometries. For this purpose, a fiber reinforced composite is modeled. There are numerous applications of textile reinforcement in several engineering fields, such as aerospace, automotive or civil engineering. Additionally to the reinforcement with fibers, the use of laminar textile structures, such as woven fabrics or laid webs, has become more popular in the past years. With the proposed approach, the complex geometry of a woven fabric (see Figure 23 (left)) can now be directly analyzed using the original surface representation. As the geometry is not star-shaped, subdivisions are required to apply the scaled boundary parameterization. Note that these are defined manually here, as the algorithmic treatment of the star-shaped subdivison is out of the scope of this work. The fiber rovings are idealized here as circular with radius R = 1 m. Note that due to the manufacturing process they usually have an elliptical shape. The woven fabric is embedded in a 10 m × 10 m × 10 m matrix material, see Figure 23 (middle). It should be noted, that the dimensions are only exemplary. To obtain star-shaped geometries for the woven fabric and surrounding matrix material, each roving is partitioned into 6 × 4 subdomains. The partitioned geometry of the fiber rovings is visualized by black lines in Figure 23 (left). In total 108 subdomains and 1120 patches are needed to discretize the structure with star-shaped domains. For each patch the surface polynomial degree is set to p b = 5 for both circumferential directions. The polynomial degree p c in scaling direction is adapted with p-refinement.
To validate the discretization strategy for complex shapes a tension test is performed. As shown in Figure 23 (right), one surface is fully fixed, while the opposite surface is loaded with a surface traction p = 10 N/m 2 . The points A and B are introduced to evaluate the displacements, where A is a corner point and B is located at the center of a roving. Firstly, it is assumed that the roving and matrix material are identical to compare the displacement to the analytical solution. For simplicity, a linear elastic material behavior is assumed with elasticity modulus E = 100 N/m 2 and Poisson's ratio = 0.0. The solution of the proposed method is compared to the analytical solution. Additionally, the model from Rhinoceros 3D is imported into ANSYS Workbench 53 to obtain a comparable numerical solution. The results are shown in Figure 24. The employed finite element mesh can be seen in Figure 24 (right). For the proposed method the analytical displacement of u y = 1 m is exactly obtained, while for the standard FEM a small error of 0.1% is observed, presumably due to the approximation of the geometry. Note that FEM is modeled with quadratic solid elements. For the proposed approach, a linear approximation in scaling direction with p c = 1 is sufficient to obtain the exact result. This leads to 40,848 equations, which have to be solved with the proposed approach. For the finite element model, 950,715 equations are required for the solution.
In a next step, the roving is assumed to be 100 times stiffer compared to the matrix material (E rov = 10000 N/m 2 ). Note that no analytical solution is available for this case. For the proposed approach the degree of the NURBS functions in scaling direction is increased with p-refinement. The degree of the boundary NURBS is kept constant at p b = 5. Figure 25 shows the displacement of point A and point B in loading direction for increasing polynomial degree p c . It is observed, that the displacement converges with p-refinement. Note that FEM is not discretized in scaling direction, therefore is omitted here. In Table 3, the required number of equations and the resulting displacement values in loading direction for points A and B are summarized for the different models.  For the finest mesh with the proposed approach, 369,840 equations are required for the solution. For comparison, the finite element model with quadratic solid elements achieves comparable accuracy with 950,715 equations. Moreover, the displacement contour for the converged solution with the finest proposed mesh is displayed in Figure 26 and compared to the solution with FEM. It can be observed that the displacement of the rovings is remarkably lower due to the increased stiffness. Altogether the results indicate, that the proposed method is capable of analyzing complex geometries with fewer number of equations, while obtaining comparable results to the finite element models.

CONCLUSIONS
This contribution is concerned with an isogeometric analysis framework that aligns with the boundary representation modeling technique in CAD. A nonlinear formulation is derived in three dimensions. The method is further extended to dynamic problems. To fit the boundary representation in CAD, we employ the parameterization of the SBFEM. Thus, the solid is partitioned into sections in relation to its boundary surfaces. It is distinguished between circumferential parameters on the boundary and a radial scaling parameter in the interior of the solid. Following the isogeometric paradigm, the same NURBS basis functions describe the geometry and approximate the solution on the boundary. In the interior of the solid, we also use B-splines for the approximation to facilitate the solution of nonlinear problems. The mass and damping matrices are derived using the same basis functions. The Galerkin method forms the basis for the solution on the boundary and also in the interior of the solid. The derived linearized operator is employed to obtain the numerical solution based on the iterative Newton-Raphson scheme. Several examples of increasing complexity demonstrated the capabilities of this approach. For eigenfrequency and dynamic analysis, higher accuracy was achieved by order elevation and k-refinement. It was shown that the solution of nonlinear problems required less degrees of freedom for both static and dynamic cases. An exemplary free-form geometry with arbitrary number of boundaries showed the application to geometries that differ topologically from a cube. Furthermore, a fiber reinforced composite was used as demonstration of a complex, nonstar-shaped geometry. It turns out that these geometries can be analyzed based on the exact surface representation while requiring less degrees of freedom to achieve comparable accuracy with finite element computations. The main advantages of the presented approach are summarized: • Feasible meshing and analysis of volumetric domains in alignment with the boundary representation in CAD.
• Applicable to wide class of problems in solid mechanics including strongly nonlinear behavior and dynamic response.
• Results as good as other numerical methods while requiring less degrees of freedom.
In this article, we demonstrated an exemplary application to a complex geometry. However, a general algorithmic framework for the star-shaped decomposition of volumetric domains is still an open issue. This would allow the application of the approach to complex, real-world problems. In this work we focused only on the accuracy with respect to the degrees of freedom, as the computational cost depends to some extent on the implementation. One important aspect here is the nonlinear case, which is computationally demanding as the system of equations is solved multiple times. Besides the number of iterations, further aspects such as preconditioning, the cost of forming the system matrices and solving the equations are not treated here, as an elaborate study of the computational cost is out of the scope of this work. Moreover, as adjacent patches are only C 0 -continuous, future work might be concerned with preserving the global continuity of the solution. In particular, stress and vibration analysis could benefit from smooth solutions. Finally, the treatment of trimmed NURBS comprises an interesting direction for future research in order to extend the applicability to trimmed models designed with the boundary representation modeling technique in CAD.

APPENDIX. DEFINITION OF MATRICES
The submatrices of Equation (5) are given in matrix notation as where j (···) denotes the entries of the inverted Jacobian J −1 .