Isogeometric Analysis

The exact solution of an engineering problem is often impossible to be found. Consequently, the engineer is only able to approximate it. The major difference between engineers and mathematicians is that engineers don’t seek for the exact solution, but for an accurate enough one to satisfy certain criteria and regulations. The real challenge is to balance between accuracy and time. The more advanced tools and techniques an engineer has in his quiver, the more optimized his projects will be.

Isogeometric analysis is the powerful generalization of the classical finite element analysis (FEA). It is a recently developed computational approach offering the possibility of fully integrating finite element analysis into conventional NURBS-based CAD tools. To date, it is necessary to convert data between CAD and FEA software packages to analyse new designs during development, a demanding task as the computational geometric approach for each is different. This new technique employs complex NURBS geometry (most CAD packages are based on it) in the FEA application directly. This lets models to be designed, examined and adjusted in one effort, using a common data set.

Geomiso is a powerful tool for modern engineers, as it makes them capable of efficiently dealing with demanding projects. It does not simply offer the new isogeometric analysis, but a fully integrated CAD/ CAE solution.

The examined structures are designed as NURBS models in physical space (Cartesian system). All calculations take place in the parameter space, already known from the isoparametric concept of the Finite Element Method (FEM), where all 3D models are represented as rectangular parallelepipeds regardless their geometry. Isogeometric analysis takes advantage of an extra space, the so-called index space, which plays a significant role for some kinds of SPLines, such as T-SPLines. It is worth mentioning that index space is only auxiliary for NURBS.

​​​Index Space

Control points are defined as the center of their support. The support of each control point is its domain of influence consisting of p+1 knot value spans for 1D, (p+1)rectangles for 2D and (p+1)cubes for 3D. Index space indicates the:

• interconnection between the control points
• support of a basis function
• contribution of a knot value to the basis
• overlapping between the finite elements

1D, 2D, 3D models are represented as lines, rectangles and rectangular parallelepipeds respectively, while knots are equally spaced regardless their value.

Figure 1. 2D model in index space​.​

Parameter Space

Parameter space is the representation of the model with respect to knots. SPLine entities are represented as lines, rectangles and rectangular parallelepipeds, which are transformed into real geometries in physical space through the Jacobian transformation. This mapping is already known from FEM. In the 1D case, each knot corresponds to the starting point of a specific knot span and to the ending point of the next one. Support is the domain, in which the basis function has non-zero value. Basis functions with intersected supports overlap in parameter space and control a shared domain of the model in physical space.

​Figure 2. 2D model in parameter space​.

Physical Space

Physical space, the well-known Cartesian space, is the one where the real geometry of the examined structure is represented. The lines, rectangles and rectangular parallelepipeds of the parameter space are transformed into real geometries in physical space. Both the control point variables (physical coordinates) and the basis play a significant role in this mapping. It is worth mentioning that for a given set of control points, only a single set of basis functions leaves the geometry intact.

Control points are often located outside the model. This is why they are not material points and don’t belong to the model in contrast to FEM’s nodes, which always belong to the mesh. Due to this, isogeometric analysis is able to utilize the existing finite element mesh of the NURBS model, something impossible in the case of the classical FEA.

1D

x(xi)=sum_(i=1)^N{R_i(xi)*X_i}={R_i(xi)}_((1xN))^T*{X}_((Nx1))
y(xi)=sum_(i=1)^N{R_i(xi)*Y_i}={R_i(xi)}_((1xN))^T*{Y}_((Nx1))
z(xi)=sum_(i=1)^N{R_i(xi)*Z_i}={R_i(xi)}_((1xN))^T*{Z}_((Nx1))

2D

x(xi,eta)=sum_(i=1)^N{R_i(xi,eta)*X_i}={R_i(xi,eta)}_((1xN))^T*{X}_((Nx1))
y(xi,eta)=sum_(i=1)^N{R_i(xi,eta)*Y_i}={R_i(xi,eta)}_((1xN))^T*{Y}_((Nx1))
z(xi,eta)=sum_(i=1)^N{R_i(xi,eta)*Z_i}={R_i(xi,eta)}_((1xN))^T*{Z}_((Nx1))

3D

x(xi,eta,zeta)=sum_(i=1)^N{R_i(xi,eta,zeta)*X_i}={R_i(xi,eta,zeta)}_((1xN))^T*{X}_((Nx1))
y(xi,eta,zeta)=sum_(i=1)^N{R_i(xi,eta,zeta)*Y_i}={R_i(xi,eta,zeta)}_((1xN))^T*{Y}_((Nx1))
z(xi,eta,zeta)=sum_(i=1)^N{R_i(xi,eta,zeta)*Z_i}={R_i(xi,eta,zeta)}_((1xN))^T*{Z}_((Nx1))

Figure 3. 2D model in physical space.

BSPLine Basis Functions

Given a sequence of non-decreasing numbers Ξ= {ξ_1,ξ_2,…ξ_(n+p),ξ_(n+p+1)}, the value of basis functions at xiin[xi_1,xi_(n+p+1)] are calculated using the Cox-de Boor recursive formula [1]:

p=0

N_(i,0)(ξ)={(1,if ξ_i<=ξ<ξ_(i+1)),(0, otherwise):}

N_(n+p,0)(ξ)={(1,if ξ_(n+p)<=ξ<=ξ_(n+p+1)),(0, otherwise):}

p=1,2,…

N_(i,p)(ξ)=(ξ-ξ_i)/(ξ_(i+p)-ξ_i)*N_(i,p-1)(ξ)+(ξ_(i+p+1)-ξ)/(ξ_(i+p+1)-ξ_(i+1))*N_(i+1,p-1)(ξ)

with the assumption of 0/0-=0.

NURBS Shape Function

W(ξ) is the Z-coordinate of the projective B-SPLine curve. Projective transformation is applied by dividing the other two coordinates of the B-SPLine curve with the Z-coordinate. NURBS shape functions are calculated as follows:

1D

R_i^p(xi)=(N_(i,p)(xi)*W_i)/(W(ξ))=(N_(i,p)(xi)*W_i)/(sum_(i’=1)^n{N_(i’,p)(xi)*W_(i’)})

W(xi)=sum_(i=1)^n{N_(i,p)(xi)*W_i}

2D

R_(i,j)^(p,q)(ξ,η)=(N_(i,p)(xi)*M_(j,q)(η)*W_(ij))/(sum_(i’=1)^n sum_(j’=1)^m{N_(i’,p)(xi)*M_(j’,q)(η)*W_(i’j’)})

W(ξ,η)=sum_(i’=1)^n sum_(j’=1)^m{N_(i’,p)(xi)*M_(j’,q)(η)*W_(i’j’)}

3D

R_(i,j,k)^(p,q,r)(ξ,η,ζ)=(N_(i,p)(xi)*M_(j,q)(η)*L_(k,r)(ζ)*W_(ijk))/(sum_(i’=1)^n sum_(j’=1)^m sum_(k’=1)^1{N_(i’,p)(xi)*M_(j’,q)(η)*L_(k’,r)(ζ)*W_(i’j’k’)})

W(ξ,η,ζ)=sum_(i’=1)^n sum_(j’=1)^m sum_(k’=1)^1{N_(i’,p)(xi)*M_(j’,q)(η)*L_(k’,r)(ζ)*W_(i’j’k’)}

Elasticity Matrix

1D

[E]_((1×1))=E

2D

[E]_((3×3))=E/(1-v^2)*[[1,v,0],[v,1,0],[0,0,(1-v)/2]] (Plane Stress)

[E]_((3×3))=E/((1-v)*(1-2v))*[[1-v,v,0],[v,1-v,0],[0,0,(1-2v)/2]] (Plane Strain)

3D

[E]_((6×6))=E/((1-v)*(1-2v))*[[1-v,v,v,0,0,0],[v,1-v,v,0,0,0],[v,v,1-v,0,0,0],[0,0,0,(1-2v)/2,0,0],[0,0,0,0,(1-2v)/2,0],[0,0,0,0,0,(1-2v)/2]]

Jacobian Matrix

1D

[J(ξ)]_((1×1))=[[R_(1,ξ)(ξ),R_(2,ξ)(ξ),…,…,…,R_(n,ξ)(ξ)]]_((1xn))*[[X_1],[X_2],[.],[.],[.],[X_n]]_((nx1))

[J(ξ)]_((1×1))^(-1)=1/[J(ξ)]_(1×1)

2D

[J(ξ,η)]_((2×2))=[[R_(1,ξ)(ξ,η),R_(2,ξ)(ξ,η),…,…,…,R_(N,ξ)(ξ,η)],[R_(1,η)(ξ,η),R_(2,η)(ξ,η),…,…,…,R_(N,η)(ξ,η)]]*[[X_1,Y_1],[X_2,Y_2],[.,.],[.,.],[.,.],[X_N,Y_N]]

[J(ξ,η)]_((2×2))^(-1)=[[J_11^(**)(ξ,η),J_12^(**)(ξ,η)],[J_21^(**)(ξ,η),J_22^(**)(ξ,η)]]=1/det[J(ξ,η)]*[[J_22(ξ,η),-J_12(ξ,η)],[-J_21(ξ,η),J_11(ξ,η)]]

det[J(ξ,η)]=J_11(ξ,η)*J_22(ξ,η)-J_21(ξ,η)*J_12(ξ,η)

3D

[J(ξ,η,ζ)]_((3×3))=[[R_(1,ξ)(ξ,η,ζ),R_(2,ξ)(ξ,η,ζ),…,…,…,R_(N,ξ)(ξ,η,ζ)],[R_(1,η)(ξ,η,ζ),R_(2,η)(ξ,η,ζ),…,…,…,R_(N,η)(ξ,η,ζ)],[R_(1,ζ)(ξ,η,ζ),R_(2,ζ)(ξ,η,ζ),…,…,…,R_(N,ζ)(ξ,η,ζ)]]*[[X_1,Y_1,Z_1],[X_2,Y_2,Z_2],[.,.,.],[.,.,.],[.,.,.],[X_N,Y_N,Z_N]]

[J(ξ,η,ζ)]_((3×3))^(-1)=[[J_11^(**)(ξ,η,ζ),J_12^(**)(ξ,η,ζ),J_13^(**)(ξ,η,ζ)],[J_21^(**)(ξ,η,ζ),J_22^(**)(ξ,η,ζ),J_23^(**)(ξ,η,ζ)],[J_31^(**)(ξ,η,ζ),J_32^(**)(ξ,η,ζ),J_33^(**)(ξ,η,ζ)]]

[[(delφ)/(delξ)],[(delφ)/(delη)],[(delφ)/(delζ)]]=[[(delx)/(delξ),(dely)/(delξ),(delz)/(delξ)],[(delx)/(delη),(dely)/(delη),(delz)/(delη)],[(delx)/(delζ),(dely)/(delζ),(delz)/(delζ)]]*[[(delφ)/(delx)],[(delφ)/(dely)],[(delφ)/(delz)]]rArr [[(delφ)/(delξ)],[(delφ)/(delη)],[(delφ)/(delζ)]]=[J(ξ,η,ζ)]_((3×3))*[[(delφ)/(delx)],[(delφ)/(dely)],[(delφ)/(delz)]]

[[(delφ)/(delx)],[(delφ)/(dely)],[(delφ)/(delz)]]=[J(ξ,η,ζ)]_((3×3))^(-1)*[[(delφ)/(delξ)],[(delφ)/(delη)],[(delφ)/(delζ)]]

where:

• R_(i,ξ)(ξ)=del/(delξ)R_i(ξ).
• R_(i,η)(η)=del/(delη)R_i(η).
• R_(i,ζ)(ζ)=del/(delζ)R_i(ζ).

Deformation Matrix

1D

{ε}_((1×1))=[(delu)/(delx)]=[J]^(-1)*[(delu)/(delξ)]hArr{ε}_((1×1))=[B_1]_((1×1))*{u_ξ}_((1×1))

[B_1(ξ)]_((1×1))=[1/(J_11)]

[(delu)/(delξ)]=[[R_(1,ξ)(ξ),R_(2,ξ)(ξ),…,…,…,R_(n,ξ)(ξ)]]*[[u_1],[u_2],[.],[.],[.],[u_n]]hArr{u}_((1×1))=[B_2]_((1xn))*{d}_((nx1))

[B_2(ξ)]_((1xn))=[[R_(1,ξ)(ξ),R_(2,ξ)(ξ),…,…,…,R_(n,ξ)(ξ)]]

[B(ξ)]_((1xn))=[B_1(ξ)]_((1×1))*[B_2(ξ)]_((1xn))

{ε(ξ)}_((1×1))=[B(ξ)]_((1xn))*{d}_((nx1))

2D

{ε}_((3×1))=[[(delu)/(delx)],[(delv)/(dely)],[(delu)/(dely)+(delv)/(delx)]]=1/det[J]*[[J_22,-J_12,0,0],[0,0,-J_21,J_11],[-J_21,J_11,J_22,-J_12]]*[[(delu)/(delξ)],[(delu)/(delη)],[(delv)/(delξ)],[(delv)/(delη)]]

[B_1(ξ,η)]_((3×4))=1/det[J]*[[J_22,-J_12,0,0],[0,0,-J_21,J_11],[-J_21,J_11,J_22,-J_12]]

[[(delu)/(delξ)],[(delu)/(delη)],[(delv)/(delξ)],[(delv)/(delη)]]=[[R_(1,ξ),0,R_(2,ξ),0,…,…,…,R_(N,ξ),0],[R_(1,η),0,R_(2,η),0,…,…,…,R_(N,η),0],[0,R_(1,ξ),0,R_(2,ξ),…,…,…,0,R_(N,ξ)],[0,R_(1,η),0,R_(2,η),…,…,…,0,R_(N,η)]]*[[u_1],[V_1],[u_2],[V_2],[.],[.],[.],[u_N],[V_N]]

[B_2(ξ,η)]_((4x2N))=[[R_(1,ξ),0,R_(2,ξ),0,…,…,…,R_(N,ξ),0],[R_(1,η),0,R_(2,η),0,…,…,…,R_(N,η),0],[0,R_(1,ξ),0,R_(2,ξ),…,…,…,0,R_(N,ξ)],[0,R_(1,η),0,R_(2,η),…,…,…,0,R_(N,η)]]

[B(ξ,η)]_((3x2N))=[B_1(ξ,η)]_((3x4N))*[B_2(ξ,η)]_((4x2N))

{ε(ξ)}_((3×1))=[B(ξ)]_((3xn))*{d}_((nx1))

3D

[B_1(ξ,η,ζ)]_((6×9))=[[J_11^(**),J_12^(**),J_13^(**),0,0,0,0,0,0],[0,0,0,J_21^(**),J_22^(**),J_23^(**),0,0,0],[0,0,0,0,0,0,J_31^(**),J_32^(**),J_33^(**)],[J_21^(**),J_22^(**),J_23^(**),J_11^(**),J_12^(**),J_13^(**),0,0,0],[0,0,0,J_31^(**),J_32^(**),J_33^(**),J_21^(**),J_22^(**),J_23^(**)],[J_31^(**),J_32^(**),J_33^(**),0,0,0,J_11^(**),J_12^(**),J_13^(**)]]

[B_2(ξ,η,ζ)]_((9x3N))=[[R_(1,ξ),0,0,R_(2,ξ),0,0,…,…,…,R_(N,ξ),0,0],[R_(1,η),0,0,R_(2,η),0,0,…,…,…,R_(N,η),0,0],[R_(1,ζ),0,0,R_(2,ζ),0,0,…,…,…,R_(N,ζ),0,0],[0,R_(1,ξ),0,0,R_(2,ξ),0,…,…,…,0,R_(N,ξ),0],[0,R_(1,η),0,0,R_(2,η),0,…,…,…,0,R_(N,η),0],[0,R_(1,ζ),0,0,R_(2,ζ),0,…,…,…,0,R_(N,ζ),0],[…,…,…,…,…,…,…,…,…,…,…,…],[…,…,…,…,…,…,…,…,…,…,…,…],[…,…,…,…,…,…,…,…,…,…,…,…],[0,0,R_(1,ξ),0,0,R_(2,ξ),…,…,…,0,0,R_(N,ξ)],[0,0,R_(1,η),0,0,R_(2,η),…,…,…,0,0,R_(N,η)],[0,0,R_(1,ζ),0,0,R_(2,ζ),…,…,…,0,0,R_(N,ζ)]]

[B(ξ,η,ζ)]_((6x3N))=[B_1(ξ,η,ζ)]_((6x9N))*[B_2(ξ,η,ζ)]_((9x3N))

{ε(ξ)}_((6×1))=[B(ξ)]_((6xn))*{d}_((nx1))

Stiffness Matrix

1D

[K]_((nxn))=int{[B(ξ)]_((nx1))^T*[E]_((1×1))*[B(ξ)]_((1xn))*A*det[J]}dξ

[K]_((nxn))=sum_(i=1)^(GPξ){[B(ξ_i)]_((nx1))^T*[E]_((1×1))*[B(ξ_i)]_((1xn))*A*det[J]*W_i^(GPξ)}

2D

[K]_((2Nx2N))=int int{[B(ξ,η)]_((2Nx3))^T*[E]_((3×3))*[B(ξ,η)]_((3x2N))*t*det[J]}dξdη

[K]_((2Nx2N))=sum_(i=1)^(GPξ) sum_(j=1)^(GPη){[B(ξ_i,η_j)]_((2Nx3))^T*[E]_((3×3))*[B(ξ_i,η_j)]_((3x2N))*t*det[J]*W_i^(GPξ)*W_j^(GPη)}

3D

[K]_((3Nx3N))=int int int{[B(ξ,η,ζ)]_((3Nx6))^T*[E]_((6×6))*[B(ξ,η,ζ)]_((6x3N))*det[J]}dξdηdζ

[K]_((3Nx3N))=sum_(i=1)^(GPξ) sum_(j=1)^(GPη) sum_(k=1)^(GPζ){[B(ξ_i,η_j,ζ_k)]_((3Nx6))^T*[E]_((6×6))*[B(ξ_i,η_j,ζ_k)]_((6x3N))*det[J]*W_i^(GPξ)*W_j^(GPη)*W_k^(GPζ)}

where:

• A: the area of the cross-section
• t: the thickness of the cross-section
• GPξ: the total number of Gauss points per ξ for the specific patch
• GPη: the total number of Gauss points per η for the specific patch
• GPζ: the total number of Gauss points per ζ for the specific patch
• ξ_i, η_j, ζ_k: the coordinates of the tensor product Gauss point ijk
• W_i^(GPξ), W_j^(GPη), W_k^(GPζ): the corresponding weights of Gauss points

Figure 4. Stiffness matrix.

Boundary Conditions

The following types of boundary conditions are taken into consideration:

• Dirichlet
• Neumann
• Robin​​

1D

{F}_((Nx1))= int{{R(xi)}_((Nx1))*f(xi)_((1×1))*det[J]}dξ

2D

{F}_((Nx1))= int int{{R(xi,eta)}_((Nx1))*f(xi,eta)_((1×1))*det[J]}dηdξ

3D

{F}_((Nx1))= int int int{{R(xi,eta,zeta)}_((Nx1))*f(xi,eta,zeta)_((1×1))*det[J]}dζdηdξ

Equilibrium

{F_f}=[K_(ff)]*{D_f}

Pseudo Displacements

{D_f}=[K_(ff)]^(-1)*{F_f}

Displacement Field

The equation’s solution corresponds to the displacements of the control points, the so-called pseudo-displacements. Unlike in classical FEM, control points don’t belong to the model. The displacement field is calculated from the pseudo-displacements.

1D

d(xi)=sum_(i=1)^N{R_i(xi)*D_i}={R_i(xi)}_((1xN))^T*{D}_((Nx1))

2D

d(xi,eta)=sum_(i=1)^N{R_i(xi,eta)*D_i}={R_i(xi,eta)}_((1xN))^T*{D}_((Nx1))

3D

d(xi,eta,zeta)=sum_(i=1)^N{R_i(xi,eta,zeta)*D_i}={R_i(xi,eta,zeta)}_((1xN))^T*{D}_((Nx1))

where:

• R_i is the shape function i
• D_i is the displacement of the corresponding control point i
• N is the number of control points

If a control point c is interpolatory to the curve at (ξ_c, η_c, ζ_c),  it follows that:

d(xi_c,eta_c,zeta_c)=sum_(i=1)^N{R_i(xi_c,eta_c,zeta_c)*D_i}=1*D_c=D_c

Strain Field

The multiplication of the deformation matrix by the pseudo-displacement vector results in the strain vector.

1D

{epsilon(xi)}_((1×1))= [B(xi)]_((1xN))*{D}_((Nx1))

2D

{epsilon(xi,eta)}_((3×1))= [B(xi,eta)]_((3x2N))*{D}_((2Nx1))

3D

{epsilon(xi,eta,zeta)}_((6×1))= [B(xi,eta,zeta)]_((6x3N))*{D}_((3Nx1))

Stress Field

Hooke’s constitutive law leads to the following equations.​

1D

{sigma(xi)}_((1×1))=[E]_((1×1))* {epsilon(xi)}_((1×1))=[E]_((1×1))*[B(xi)]_((1xN))*{D}_((Nx1))

2D

{sigma(xi,eta)}_((3×1))=[E]_((3×3))* {epsilon(xi,eta)}_((3×1))=[E]_((3×3))*[B(xi,eta)]_((3x2N))*{D}_((2Nx1))

3D

{sigma(xi,eta,zeta)}_((6×1))=[E]_((6×6))* {epsilon(xi,eta,zeta)}_((6×1))=[E]_((6×6))* [B(xi,eta,zeta)]_((6x3N))*{D}_((3Nx1))