An Isogeometric Framework for the Modeling of Curvilinear Anisotropic Media

NURBS-based isogeometric analysis utilizes Non-Uniform Rational B-Spline (NURBS) basis functions defined by a knot vector:

The knot vector is a non-decreasing set of coordinates in the parameter space, where $i=1,\cdots,n+p+1$, and $n$ is the number of basis functions of the order $p$ used to construct the B-spline curve, which were first proposed by Gordon and Riesenfeld [27, 28] for computer-aided graphic design. Given the knot vector in the parameter space defined above, the $i^{th}$ B-spline basis function of order $p$, denoted by $N_{i,p}(\xi)$ can be computed following Cox-de Boor recursion formula [29, 30, 31] as shown in Eq. (2) and (3):
For $p=0:$

for $p=1,2,3,\cdots,$

By taking a linear combination of the B-spline basis $N_{i,p}(\xi)$ and the corresponding control points $P_{i}\in\mathbb{R}$, the B-spline curve is described as

Similarly, assigning a control net $P_{i,j}$, where $i=1,...,n$, $j=1,...,m$, and B-spline functions in $\xi$- and $\eta$-directions $N_{i,p}(\xi)$ and $M_{j,q}(\eta)$ with the order of $p$ and $q$ constructed by knot vectors $\boldsymbol{\Xi}=[\xi_{1},...,\xi_{n+p+1}]$, and $\boldsymbol{H}=[\eta_{1},...,\eta_{m+q+1}]$, respectively, the B-spline surface is also given by

Making a rational function divided by weight functions enables exact representations for a wide range of complex geometry, especially in conic sections (circles, ellipses, etc.), which are beyond an array of piece-wise polynomials. This rational function is referred to as Non-Uniform Rational B-Spline function or NURBS (See the references e.g. [32, 33, 34]). Given the knot vector $\boldsymbol{\Xi}$ defined above, the $p^{th}$-order NURBS function on the knot span $\xi\in[\xi_{i},\xi_{i+1})$ is given by

where $\omega_{i}$ is the $i^{th}$ weight and $N_{i,p}(\xi)$ is the $i^{th}$ B-spline function with order $p$. Then, the NURBS curve and surface based on the knot vectors $\boldsymbol{\Xi}$ and $\boldsymbol{H}$ can be generated using the NURBS basis functions as shown in Eq. (6), respectively:

and

where

$R^{p,q}_{i,j}(\xi,\eta)$ is a bi-variable NURBS basis function and $\omega_{i,j}$ is a weight for the corresponding control net $P_{i,j}$.

Defining the $p^{th}$-order B-spline function $N_{i,p}(\xi)$ by the knot vector $\boldsymbol{\Xi}=[\xi_{1},...,\xi_{n+p+1}]$, the derivative $N^{\prime}_{i,p}(\xi)$, is calculated recursively as:

This can be proven by a mathematical induction [35].

Starting from the derivative of B-Splines, the derivative of NURBS functions is simply derived using the quotient rule to Eq. (6) such that:

where

Further properties and techniques for B-Splines and NURBS can be found in [36, 37, 38, 39] and [35, 32, 33, 34], respectively. For detailed theoretical framework and the implementation in IGA, the reader is referred to Hughes et al. [20, 22].

Typically, the description of curvilinear fiber paths is performed assigning a constant tangent vector in each element and considering each element orientation as a variable of the parameter space  [8, 9, 11, 12]. While this approach can be very accurate, it also leads to a significant computational cost due to the large number of design variables (the number of design variables is equal to the number of elements). Since the main goal of the optimization study proposed here is to showcase the use of IGA rather than identifying the perfect fiber configuration, we opted for the simplified approach proposed in [13, 14]. Hence, we utilize the concept of level-set function [13, 14] to present a subset of possible fiber configurations in the chosen domain, which also ensures the following important characteristics in terms of manufacturability: 1) continuity of fiber paths, and 2) absence of intersections between fibers. Therefore, in the present work, we describe the curvilinear fiber paths by using a $5^{th}$-order complete polynomial as a level-set function:

where $C_{ij}$, with $i,j=0,...,5$, are design variables of the $5^{th}$-order complete polynomial function. The fiber orientations $\theta$ are also calculated by partially differentiating the level set function of Eq. (13) with respect to $x$ and $y$ as:

It is important to restate here that the optimization study only serves the purpose of demonstrating the IGA framework. The level-set functions can only represent a subset of the fiber path solutions and cannot be exhaustive. Ongoing work by the authors is focusing on robust implementation of a comprehensive optimization framework. However, its description goes beyond the scope of the present work.

The constitutive relation of CTI media depends on the local fiber orientation. As a consequence, the element stiffness matrix must be evaluated at each Gauss point. The constitutive relations for a linear elastic problem can be expressed as:

where $\{F\}$, $\{d\}$, and $[K]$ are the loading vector, the displacement vector, and the global stiffness matrix, respectively. Then, the element stiffness matrix in the physical domain $\Omega^{(e)}$ can be written in the form:

where $[B]$ is the strain-displacement matrix generated from the NURBS basis functions, Eq. (9), and

where $[T_{\sigma}^{-1}(\theta)]$ and $[T_{\epsilon}^{-1}(\theta)]$ are transformation matrices for the stress and strain, respectively, and $[C]$ is the in-plane stiffness coefficients matrix depending on three engineering moduli ($E_{1}$, $E_{2}$, and $G_{12}$), and the Poisson’s ratio ($\nu_{12}$) [40, 41]. Then, the element stiffness is calculated by Gauss quadrature such that:

where $N_{p}$ is the number of integration points, $w_{\tilde{\xi}_{m}}$ and $w_{\tilde{\eta}_{m}}$ are the corresponding weighs for the each direction, and $J_{m}$ is the determinant of the Jacobian which maps from physical space to parent space. Detail derivation and calculation of the element stiffness in the isogeometric framework can be found in [21].

As Fig. 1 shows, the analysis of CTI media in previous works [8, 9, 13, 14, 11, 12] was performed using a uniform constitutive relation for the integration points in each element. In other words, $\theta=\theta^{(n_{e})}$, with $n_{e}$ is an element number, i.e., $[K]^{(e)}=[K^{(e)}\{\theta^{(n_{e})}\}]^{(e)}$. In contrast, in the present work we let the constitutive relation depend on the location of integration points within the element during the computation of element stiffness matrix such as $\theta=\theta(x,y)$, where $x$ and $y$ are the coordinates of the integration points, i.e., $[K]^{(e)}=[K\{\theta(x,y)\}]^{(e)}$ [42]. This results in accounting for the effects of the gradient in fiber curvature in calculating the element stiffness.

For the design optimization study in this paper, we adopted Sequential Quadratic Programming (SQP) assuming that the object function is twice differentiable. The SQP is one of the most powerful optimization methods designed for Constrained Non-Linear Programming (NLP) problems [43, 44]. This enables to achieve an attractive rate of convergence using $2^{nd}$-order derivatives of Lagrange function.

In this paper, the SQP follows the Constrained Steepest-Descent (CSD) method via two major procedures such that 1) Quadratic Programming (QP) subproblem and 2) Quasi-Newton method. Let define an NLP problem of the form:

Minimize

subject to

In the SQP method, QP subproblem is solved for a search direction $\boldsymbol{d}$ in a CSD framework. Thus, Karush-Kuhn-Tucker necessary conditions (KKT conditions) [45, 46] can be applied utilizing the Newton-Raphson method, and a general constrained optimization problem $\bar{f}$ is derived as:

Minimize

subject to the constraints:

where $\boldsymbol{H}$ is the Hessian matrix of Lagrange function $\nabla^{2}{L}$, $\boldsymbol{n}^{(i)}$ and $\boldsymbol{a}^{(i)}$ are the gradients of the functions $h_{i}$ and $g_{i}$, respectively, $\boldsymbol{c}$ is the gradient of the objective function of Eq. (19), and $\boldsymbol{x}^{(k)}$ represents the current iteration point. Thus, the search direction $\boldsymbol{d}$ is computed by the gradient conditions of the KKT necessary conditions.

Once the QP subproblem has been derived, the Hessian matrix $\boldsymbol{H}$ must be approximated for each iteration of the CSD framework. The Hessian matrix $\boldsymbol{H}$ is approximated using Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [47, 48, 49, 50]. Define the gradient difference vector $\boldsymbol{y}^{(k)}$ of the Lagrange function at two iteration points and the design change vector $\boldsymbol{s}^{(k)}$ with the Lagrange multipliers $\boldsymbol{v}^{(k)}$ and $\boldsymbol{u}^{(k)}$ at the current point such that

and set the vector $\boldsymbol{r}^{(k)}$ with the scalar quantity $\theta$ as:

Then, the Hessian matrix can be updated as:

Finally, the solutions of Eq. (19) can be obtained by the use of above two methods 1) SQP subproblem and 2) Quasi-Newton method in the CSD framework. Further explanations including schematic algorithms on these methods can be found in [44, 43].

To demonstrate the IGA framework, the goal of our optimization study is obtaining optimal curvilinear or linear fiber paths with a minimum stress concentration or Tsai-Wu failure index. Thus, assuming that our objective functions $K_{t}$ for a minimum stress concentration factor and a minimum Tsai-Wu failure index $\Phi$ are both twice differentiable nonlinear functions, the optimal design problems can be defined as

Minimize

where $C_{ij}\in\mathbb{R}$ are the coefficients of the level-set function defined in Eq. (13). Thus, the objective functions are subject to the fiber orientations in Eq. (14). In this study, the SQP method was implemented in the MATLAB environment utilizing built-in function called fmincon with an option of sqp algorithm [51].

Stress analysis on an elastic plate weakened by a semi-circular hole notch is conducted using both NURBS-based IGA and classical FEA assuming a linear elastic problem illustrated in Fig. 2. The semi-circular notched plate is under a unit constant tension of $T_{y}=1$ MPa along the top edge and assumed as a plane-stress problem with a unit thickness. The plate is constrained along the bottom and the right-hand side edges in the $y$- and the $x$-directions, respectively. The elastic plate is modeled using quadratic NURBS basis functions in the IGA simulations with the two knot vectors $\boldsymbol{\Xi}=[0,0,0,0.25,0.5,0.5,0.75,1,1,1]$ and $\boldsymbol{H}=[0,0,0,1,1,1]$, and mesh refinement is implemented using knot insertion [20, 22]. In conventional FEA, simulations are conducted applying quadrilateral 9-node element (Q9) and, the exact same node coordinates used in IGA are employed for the sake of fair comparison in terms of convergence. The semi-circular notched plate with several fiber configurations is examined by the following properties:

1. 1.

   Stress concentration factor $K_{t}$:

   $$K_{t}=\sigma_{y}^{max}/T_{y},\quad\sigma_{y}^{max}\text{: Maximum longitudinal stress found in the body,}$$

   (29)

2. 2.

   Average in-plane stiffness along the top edge $\bar{K}$:

   $$\bar{K}=\frac{P_{tot}}{\bar{v}},\quad\text{with}\quad P_{tot}=\int_{0}^{L}\sigma_{yy}d\,x\quad\text{and}\quad\bar{v}=\frac{1}{L}\int_{0}^{L}vd\,x,$$

   (30)

   where $\sigma_{yy}$ and $v$ are the longitudinal stress and displacement along the top edge, respectively, and $L$ is the length of the top edge.

3. 3.

   Tsai-Wu failure index $\Phi$ assuming Transverse Isotropy under plane-stress condition [52]:

   $$\Phi=F_{1}\sigma_{11}+F_{2}\sigma_{22}+F_{11}\sigma_{11}^{2}+F_{22}\sigma_{22}^{2}+F_{66}\tau_{12}^{2}+2F_{12}\sigma_{11}\sigma_{22},$$

   (31)

   where

   $$\begin{split}F_{1}&=\frac{1}{\sigma_{1t}^{f}}-\frac{1}{\sigma_{1c}^{f}},\qquad F_{2}=\frac{1}{\sigma_{2t}^{f}}-\frac{1}{\sigma_{2c}^{f}},\qquad F_{66}=\frac{1}{\tau_{12}^{f}},\\ F_{11}&=\frac{1}{\sigma_{1t}^{f}\sigma_{1c}^{f}},\qquad F_{22}=\frac{1}{\sigma_{2t}^{f}\sigma_{2c}^{f}},\qquad F_{12}=-\frac{1}{2}\sqrt{F_{11}F_{22}}.\\ \end{split}$$

   (32)

The choice of this criterion is motivated by its effectiveness and accuracy in capturing the failure condition of composite laminates featuring smooth or slightly notched surfaces. However, for more complex structural configurations and loading conditions, the emergence of large Fracture Process Zones (FPZs) can lead to significant size effects. This is typically the case of structures featuring open and filled holes, sharp notches and other stress raisers. In such a case, a quasibrittle fracture mechanics framework should be preferred over a stress-based failure criterion [53, 54, 55, 56, 57]. All the required material properties such as longitudinal and transverse elastic moduli $E_{1}$ and $E_{2}$, shear modulus $G_{12}$, Poisson’s ratio $\nu_{12}$, and the ply strengths: $\sigma_{1t}^{f}$, $\sigma_{1c}^{f}$, which are tensile and compressive strengths in the fiber direction, $\sigma_{2t}^{f}$, $\sigma_{2c}^{f}$, which are tensile and compressive strengths perpendicular to the fiber directions, and $\tau_{12}^{f}$, that is shear strength are respectively assigned and tabulated in Table 1.

In the present study, four different fiber configurations as described in Fig. 4 were investigated: a) conventional longitudinal straight fiber paths, b) conventional transverse straight fiber paths, c) concentric fiber paths following a semi-circular notch, and d) curvilinear fibers following the holomorphic path defined by the conformal mapping presented in [19]. The four fiber orientations measured from the global $x$-axis are defined as a function of $(x,y)$ coordinate such that:

1. (a)

   Conventional longitudinal straight fibers: $\theta(x,y)=\pi/2,$

2. (b)

   Conventional transverse straight fibers: $\theta(x,y)=0,$

3. (c)

   Concentric fibers following the semi-circular hole:

   $$\theta(x,y)=\frac{\pi}{2}+\cos^{-1}{\frac{x}{a}},\quad a=\sqrt{x^{2}+y^{2}},$$

4. (d)

   Curvilinear fibers following the holomorphic function presented in  [19]:

   $$\theta(x,y)=\tan^{-1}\Bigg{[}\frac{1+\frac{R^{2}}{x^{2}+y^{2}}\cos\big{(}2\tan\frac{y}{x}\big{)}}{-\frac{R^{2}}{x^{2}+y^{2}}\sin\big{(}2\tan\frac{y}{x}\big{)}}\Bigg{]}.$$

Figure 5 shows the error computed in $L^{2}$-norm of stresses for each of the fiber configurations presented in Fig. 4. The error for each component is defined by the following:

where $V$ is the domain of the whole material, and the $\boldsymbol{\epsilon}_{\sigma}$ is the vector of the difference between the sufficiently converged stress field (change in the $K_{t}$ less than $0.5$%) and the current stress field. Then, the error $\epsilon$, the root-mean-value of the Eq. (33) is taken as:

where $n$ is the number of components for the stress field and “$\cdot$” represents the dot product here.

In addition, we also focused on the continuity of the basis functions for both IGA and FEA. In the present study, quadratic Lagrange polynomials are adopted as the basis function in the FEA framework, which leads to $C^{0}$-continuity between elements all the time. On the contrary, the basis function of NURBS-based IGA is decided by the quadratic NURBS basis function with $C^{1}$-continuity everywhere except along the $x$-axis due to the multiplicity of the knot values at $\xi=0.5$. The enhanced smoothness of the interpolating functions element by element results in capturing the field variables between elements in better numerical performance for a given number of degrees-of-freedom compared to classical FEA. Moreover, as can be seen in Fig. 5d, the difference between the errors for IGA and FEA is way greater than the ones of the conventional fiber configurations (Fig. 5a and 5b), which indicates that IGA especially plays a significant role in the CTI specimens.

In summary, IGA provides:

1. 1.

   time reduction

2. 2.

   computational efficiency

3. 3.

   more accurate description of the field variables for a given number of degrees-of-freedom, compared to classical FEA

These features make IGA a valuable tool for the simulation of structural and electrostatics problems in media featuring curvilinear anisotropy systems [19, 8, 9, 10, 11, 12, 13, 14, 15, 16].

As can be seen in Table 2, generally, conventional longitudinal fiber reinforced composites have high-level of in-plane stiffness with significantly high stress concentration around notches and other stress raisers. In contrast, reinforced composites featuring transverse straigth fibers have lower stress concentration and in-plane stiffness. Leveraging the IGA theoretical framework, here we explore how the different curvilinear fiber paths described in Fig. 4 affect stress concentration and stiffness. All the stress concentration factors and Tsai-Wu failure indexes obtained from simulations are tabulated in Table 2.

The optimization study defined in Eq. (26) and (27) is conducted in order to find optimal fiber paths which minimize the stress concentration factor and the Tsai-Wu failure index via SQP method on the elastic plate defined in Fig. 2 and Table 1. As initial condition, the parameters $C_{ij}$ of the level-set function describing the family of fiber paths, Eq. (13), are chosen to take unit values.

The optimal fiber paths identified via SQP are shown in Fig. 11a, and the results in terms of stress concentration, failure index, and overall stiffness are summarized in Table 3. The contours of the stresses are shown in Fig. 11b-d. The stress concentration for the optimized case 1.28 (Table 2a) which is remarkably 82 % lower than the longitudinal fiber case and even lower than the isotropic case. This result implies that, by harnessing curvilinear anisotropy, it is possible to significantly reduce the severity of the semi-circular notch. This cannot be achieved in isotropic materials unless the elastic properties are properly functionally-graded. It is interesting to note that the obtained stress concentration is even lower than in the transverse straight fiber case, which features the lowest stress concentration among the configurations previously investigated. In fact, the stress concentration for the optimal fiber path case (Table 2a) is lower by 46 %. In conclusion, harnessing curvilinear anisotropy can significantly reduce the stress concentration even compared to isotropic case. Factoring the increase in strength and stiffness provided by the embedded fibers, it is clear that media featuring CTI can surpass the structural performance of traditional composites. Future work will focus on translating these results to other notch configurations including e.g. sharp V-notches or U-notches.

The optimization study for the minimum Tsai-Wu index is also implemented under the same condition for the minimum stress concentration. The initial design variables $C_{ij}$ are decided as the longitudinal straight fiber path shown in Fig. 4a such that:

because this fiber path provides the minimum Tsai-Wu failure index based on the four different fiber configurations (Fig. 4a and Table 2).

As a result of the optimization, the design variables remained as same as the initial conditions indicating that the longitudinal straight fiber case is the one minimizing the failure index (Table 2). In terms of Tsai-Wu index, conventional longitudinal straight fiber path is the best option, however, the stress concentration is very significant compared to the optimal fiber case of the minimum stress concentration. Thus, it is interesting to make a comparison between the two focusing on Tsai-Wu index as discussed next.

Based on Table 3, the maximum Tsai-Wu failure index for the fiber path that provides the minimum stress concentration is found as $0.020$ and it is located at the very top-left node. This value is about five times higher than the one from the longitudinal straight fiber path, which is the optimal fiber path for minimum Tsai-Wu index. Table 4 includes the local stresses of the each component, and the local dominance of the ply strength at the points of the maximum Tsai-Wu failure index. The local ratios between stresses to the lamina strength such that:
For the longitudinal fiber direction:

For the transverse fiber direction:

while for the shear component the ratio, $\tau_{12}/\tau_{12}^{f}$ is independent of the sign of the local shear stress. Considering the case of the optimal fibers for the minimum stress concentration (Fig. 11a), the transverse stress component is the closest to the related strength value, showing that matrix failure by e.g. splitting at the very top-left of the plate (Fig. 12a) is the most likely damage mechanism. However, the resulting crack would probably grow along the curvilinear fiber path and eventually be arrested. Thus, even though the plate begins to experience the failure, it might not immediately be exposed to be an ultimate failure. In the medium optimized for minimum Tsai-Wu index, the transverse and the shear components of the strengths take the highest values compared to the respective strengths, which implies that splitting crack is again the dominant failure mode. However, as can be seen in Fig. 12b, the splitting crack will initiate close to the notch tip and propagate undisturbed following the longitudinal fiber path until reaching the end of the plate.

Thus, it is very difficult to predict which fiber configuration will attain the highest load capacity. Of course, the Tsai-Wu failure criterion is not equipped to predict stable crack growth and quasibrittle fracture mechanics models will be used in future works to investigate the highly nonlinear progressive damage following initiation [57, 54, 58, 59, 60].