Isogeometric shape design sensitivity analysis of stress intensity factors for curved crack problems

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496 www.elsevier.com/locate/cma

Isogeometric shape design sensitivity analysis of stress intensity factors for curved crack problems Myung-Jin Choi, Seonho Cho ∗ National Creative Research Initiatives (NCRI) Center for Isogeometric Optimal Design, Department of Naval Architecture and Ocean Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 151-744, Republic of Korea Received 8 September 2013; received in revised form 1 July 2014; accepted 5 July 2014 Available online 14 July 2014

Highlights • • • • •

Exact curvature and normal vector are embedded in the isogeometric shape sensitivity of SIFs. Exact geometry using the NURBS precisely determines the local coordinate system at crack-tip. Configuration design variation is naturally considered due to the local coordinate system at crack-tip. Design dependence of auxiliary field and weight is considered in the shape design sensitivity. The crack-face integral and interaction integral are essential for SIFs in curved crack problems.

Abstract An isogeometric shape design sensitivity analysis method is developed for the stress intensity factors (SIFs) in curved crack problems. In order to obtain an enhanced shape sensitivity of SIFs, exact higher-order geometric information of curvature and normal vector is embedded in the isogeometric shape sensitivity expressions. A direct differentiation method is employed for the design sensitivity analysis and the size and orientation of the crack are selected as design variables. In the isogeometric approach, the NURBS basis functions in CAD system are directly utilized in the response analysis, which enables a seamless incorporation of exact geometry and higher continuity into the computational framework. The precise evaluation of the crack-face integral as well as the interaction integral is essential for the computation of the stress intensity factors for the curved crack problems in a mixed-mode. The CAD-based exact representation of tangential and normal vectors enables us to exactly define a local coordinate system at the crack-tip, whose shape dependency naturally leads to configuration design variations that include the change of crack orientation. Compared with the conventional finite element approach, a higher continuity of stress and strain fields are expected in the interaction integral domain. The design dependencies of auxiliary field and q-function are considered in the design sensitivity formulation, reflecting the convection effects due to the non-constant design velocity in the interaction integral domain. Various numerical examples of curved crack problems are presented to verify the developed isogeometric sensitivity analysis method through the comparison with exact and finite element solutions. c 2014 Published by Elsevier B.V. ⃝

Keywords: Isogeometric analysis; Shape design sensitivity analysis; Stress intensity factor; Curved crack; Crack-face integral; Configuration design

∗ Corresponding author. Tel.: +82 2 880 7322; fax: +82 2 880 7322.

E-mail address: [email protected] (S. Cho). http://dx.doi.org/10.1016/j.cma.2014.07.002 c 2014 Published by Elsevier B.V. 0045-7825/⃝

470

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Nomenclature X I , (I = 1, 2), xi , (i = 1, 2) Global and local rectangular Cartesian coordinates Xc Position of crack-tip in global Cartesian coordinates N I , (I = 1, 2), n i , (i = 1, 2) Normal vectors in global and local Cartesian coordinates at crack-tip Γc+ , Γc− Upper and lower crack-faces θ c , κ c Tangential angle and curvature of crack-face curve at crack-tip 1. Introduction The identification of crack initiation and growth is crucial for evaluating structural safety. Recently, Nam et al. [1] applied the control of crack initiation, propagation, and termination to microscopic pattern generation. The controlling method overcame the problem of deteriorated resolution of patterns due to electron diffraction in the existing lithography method by an electron beam, together with the significant reduction of operation time for pattern generation. Through numerous experiments, this research shows the feasibility of the micro-notch design in obtaining a desired crack path. In linear elastic fracture mechanics, there are three linearly independent cracking modes of opening, sliding, and tearing. The stress intensity factors (SIFs) for those modes are used to predict the stress state or stress intensity near the crack tip and provide a failure criterion for brittle materials. The magnitude of SIFs depends on the geometry of specimen, the size and location of crack, and the magnitude and the modal distribution of loads on the material. Among various methods to compute the SIFs, we select the method of M-integral or interaction integral derived from the J -integral which is equivalent to the energy release rate in linear elastic material. To evaluate the J -integral, we introduce the “q-function” that is a weight function ranging from 1 at the crack tip to 0 on the boundary of J -integral. For more details, interested readers may see Refs. [2,3]. For the precise control of the crack formation and propagation, it is necessary to predict the initiation and stability of crack growth using the variation rate of SIFs with respect to the design parameters such as the size, shape, and orientation of cracks. There are some research works on design sensitivity analysis (DSA) of SIFs only for straight crack problems in linear fracture mechanics. However, for the curved crack problems, no publication is reported yet regarding the DSA of SIFs even though many research works were reported for SIFs analysis. Chen et al. derive the first order shape sensitivity of J -integral using a direct differentiation method to apply to the first-order reliability method for probabilistic fracture-mechanics (see [4]) and suggest a shape DSA method of interaction integral formulation using the direct differentiation method for mixed-mode straight crack problems (see [5]). However, they do not consider the shape dependency of the auxiliary field and the q-function. Rao and Rahman [6] extend the DSA method to functionally graded materials (FGMs), using a homogeneous or non-homogeneous auxiliary field. They derive the shape sensitivity of the q-function but do not consider the shape dependency of the auxiliary field. A constant design velocity field called “virtual rigid body translation” is used in the M-integral domain so that the shape dependency of only an actual field is considered. When it comes to curved crack problems, the design velocity field is no longer constant and should be defined according to the shape dependency. It is one of the contributions in this paper that suggests analytical DSA formulation for actual as well as auxiliary fields in interaction integral formulation to consider the non-constant design velocity field in the interaction integral domain, which is essential in the DSA of curved crack problems. Feijoo et al. [7] employ the shape DSA to obtain the first order derivative of total potential energy which is the strain energy release rate. Taroco [8] extends the research to the second order shape derivative of total potential energy which is the change of strain energy release rate as a crack grows. However, the second-order derivatives of stress and strain are required, which is not appropriate for the C 0 finite element methods. When a crack growth is considered as the shape variation of the domain, the concept of shape DSA is also used to calculate the energy release rate and the mixed-mode SIFs of isotropic FGMs (see [9]). In this paper, aforementioned researches about the analytical DSA of straight crack problems are extended to the curved crack problems using the isogeometric approach, which enables to accurately consider configuration design variation as the local coordinate system at the crack-tip changes with the curved variation of crack geometry. Although the aforementioned studies on the shape DSA of SIFs are limited to straight crack problems, there are many researches on the SIF analysis of the curved crack problems. The interaction integral formulation derived from

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

471

the J -integral by considering a superposition of actual and auxiliary fields (see [3]) is known to be the most accurate and efficient method for calculating SIFs under mixed-mode loading conditions. Gosz and Moran [10] adopted the interaction integral method to study 3D non-planar cracks in homogeneous materials. Using the interaction integral formulation, Walters et al. [11] compute the SIFs in mixed mode for 3D curved crack problems subjected to surface tractions. They point out that for the evaluation of the response in the auxiliary field, the local coordinates of the crack-tip should be highly accurate but could be inaccurate due to the geometric approximation errors if linear finite elements are used. To overcome the aforementioned difficulty, they employ a local orthogonal coordinate system with a highly refined mesh for precise computation. Yu et al. [12] extended the interaction integral method to solve the fracture problems in 3D nonhomogeneous materials with arbitrary interfaces in the integral domain. Nevertheless, the geometric approximation in numerical methods has been an important and persisting issue to overcome. Recently, taking advantage of the NURBS (Non-Uniform Rational B-Splines) basis functions, Hughes et al. [13] develop an isogeometric analysis (IGA) method that has many advantages over the standard finite element methods in various engineering applications like structural vibrations, fluid–structure interactions, turbulent flow simulations, and shell analyses (see, e.g. [14]). The major feature of the isogeometric method is to employ the same NURBS basis functions as used in the CAD systems, which gives a significant advantage of exact representation of geometry. The isogeometric analysis method has many features in common with the finite element method; it invokes the isoparametric concept where the dependent variables and the geometry share the same basis functions. Also, the method has some features in common with the meshfree methods; it is not interpolatory. In addition to the benefits of the isogeometric analysis, the isogeometric design sensitivity analysis also has the following advantages (see, e.g. [15,16]): first, it provides more accurate design sensitivity for complex geometries including higher order geometric information such as normal vector and curvature, which could be significantly helpful for the DSA of curved crack problems. The NURBS functions of higher continuity offer a much more compact representation of the response and sensitivity of systems than the standard finite element functions do, yielding more accurate results, even at the same order of basis functions. Second, it vastly simplifies the design modification without communicating with the CAD description of geometry if the developed DSA method for the curved crack problems is extended to design optimization. However, updating internal control points by maintaining the injectivity of domain parameterization during shape optimization process is still challenging. The internal control points can be updated using either an algebraic constraint (see, e.g. [17]) or an optimization scheme enforcing the Jacobian to be positive in the domain (see, e.g. [18,19]). Verhoosel et al. [20] employs the isogeometric method having higher order continuity for the approximation of higher order gradient damage, which overcomes the limitation of the conventional FEA having C 0 continuity. Verhoosel et al. [21] uses the NURBS basis function for the discretization of cohesive zone formulation, and represents the discontinuity through the insertion of knots. Also, the evolution of T-mesh is utilized to represent the growth of crack. To improve accuracy of the solution in the isogeometric method, De Luycker et al. [22] employ the Heaviside enrichment and the crack-tip enrichment using asymptotic solution to represent the discontinuity and singularity of the solution, respectively. Only a mode-I straight crack problem is considered so that there still is a limitation on the solution accuracy in curved crack problems due to geometry approximation even though a sub-triangle technique is used for the partial numerical integration of elements. Ghorashi et al. [23] also uses the enrichment and the subtriangle technique for partial numerical integration in the isogeometric analysis, and tried to solve a mixed mode straight crack. For general curved crack problems, the geometry is still approximated by straight lines and thus the accuracy of stress intensity factor is not thoroughly verified. In this paper, the curved geometry of the cracks is exactly represented by taking advantage of the isogeometric approach, and it is verified that the exact description of curved geometry of crack significantly improves the accuracy of stress intensity factors. This paper is organized as follows: in Section 2, we briefly summarize the NURBS basis function and explain how to represent the discontinuity for cracks in the CAD model. In Section 3, using the isogeometric analysis method, we explain the interaction integral formulation of curved crack problems to obtain the stress intensity factor. In Section 4, the DSA method using an interaction integral is formulated, including the shape dependency of auxiliary field, qfunction, and local coordinate system at the crack-tip. In Section 5, several numerical examples are presented to verify the advantage of developed IGA and DSA methods. Finally, the importance of developed shape DSA method and the superior points of the isogeometric method are discussed in the concluding remarks.

472

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

2. NURBS basis function 2.1. B-spline basis function In the isogeometric analysis, the solution space is represented in terms of the same basis functions used for representing the geometry. The isogeometric analysis has several advantages over the conventional finite element analysis such as geometric exactness and easy refinements due to the use of NURBS basis functions based on B-splines. Consider a set of knots ξ in one dimensional space, which consists of knots ξi in a parametric space.   ξ = ξ1 , ξ2 , . . . , ξn+ p+1 , (1) where p and n are the order of the basis function and the number of control points, respectively. The B-spline basis functions are recursively defined as  1 if ξi ≤ ξ < ξi+1 0 Ni (ξ ) = ( p = 0) (2) 0 otherwise, and p

Ni (ξ ) =

ξi+ p+1 − ξ ξ − ξi p−1 p−1 N (ξ ), N (ξ ) + ξi+ p − ξi i ξi+ p+1 − ξi+1 i+1

( p = 1, 2, 3, . . .) .

p

(3) p

Using the B-spline basis function Ni (ξ ) and weight wi , the NURBS basis function Ri (ξ ) is defined as p

p Ri (ξ )

N (ξ )wi = n i .  p N j (ξ )w j

(4)

j=1

Generally, the isogeometric approach using higher order basis functions offers higher regularity than the conventional p finite element approach. NURBS curves are obtained from the linear combination of basis functions Ri (ξ ) and the corresponding control points Bi . C(ξ ) =

n 

p

Ri (ξ ) Bi .

(5)

i=1

Similarly, NURBS surface S is defined as a tensor product of coordinates, S (ξ, η) =

l  m 

p

q

Ri (ξ )R j (η) Bi j .

(6)

i=1 j=1

For the brevity of expression, Eq. (6) can be rewritten as S(4) ≡ S (ξ, η) =

CP 

W I (4)B I

(7)

I

where CP and 4 are the number of control points and the parametric domain of surface, respectively. The constructed NURBS basis functions of order p possess the property of affine covariance and ( p − 1) continuous differentiability. If the knots are repeated k-times, the continuity of NURBS basis functions decreases by k as well. For more details of NURBS geometry, interested readers may see Rogers [24], Piegl and Tiller [25]. 2.2. Modeling of crack-tip To locate the crack-tip at the desired position, NURBS basis functions are modified to possess the Kronecker delta property by knot repetition. In the CAD model, the explicit representation of discontinuity can be achieved through the use of multiple patches. Fig. 1 illustrates how to represent discontinuity using two NURBS patches and repeated knots.

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) Physical space.

473

(b) Parametric space. Fig. 1. Representation of discontinuity in the CAD model.

In Fig. 1(a), X = 1.45 is the desired position of crack-tip in physical space, and the corresponding position in parametric space is ξc = 0.4833 which can be found through a point inversion algorithm using the Newton–Raphson iteration (see [26]). The control points B1 ∼ B5 correspond to the quadratic NURBS basis functions R12 ∼ R52 , respectively. In Fig. 1(a), to locate the crack-tip at X = 1.45, the displacement should be discontinuous in X ∈ [0, 1.45) and continuous in X ∈ [1.45, 3] across the patch boundary, which can be achieved by repeated knots at the desired position of crack-tip ξc = 0.4833 in parametric space. As in Fig. 1(b), the knot multiplicity is the same as the order of the basis function, which allows the Kronecker delta property at the desired position of crack-tip and enables the two sets of basis functions R12 ∼ R32 and R32 ∼ R52 to have a local support in two regions ξ ∈ [0, 0.4833] and ξ ∈ [0.4833, 1], respectively. The continuity across the crack-tip position ξc = 0.4833 is decreased to C 0 due to the repeated knots. However, in the other regions ξ ∈ [0, 0.4833) and ξ ∈ (0.4833, 1] of each NURBS patch, C 1 -continuity is still maintained. Even though a loss of continuity is present at the crack-tip in the ξ -direction, the other regions ξ ∈ [0, 0.4833) and ξ ∈ (0.4833, 1] maintain the given continuity. In the η-direction, there is no loss of continuity. Furthermore, J -integral is to be employed so that there is no need to evaluate the stress field at the crack-tip precisely. Thus, the scheme of repeated knots will be employed hereafter. It is indeed difficult to exactly represent the complicated cracks using the interface between NURBS patches, which could be an important and challenging research topic. For general crack propagation problems, the exact solutions are usually unavailable and the precise computation of SIF is very difficult. However, the development of more precise numerical method is necessary regardless of modeling accuracy. In this paper, we are not trying to solve the complicated crack problems but to develop more precise numerical methods for the design sensitivity analysis as well as the response analysis for fracture mechanics. 3. Isogeometric analysis of SIFs in curved crack problems 3.1. Review of isogeometric analysis Consider the plane elasticity problem shown in Fig. 2. The elastic body occupies an open domain Ω bounded by a closed boundary Γ . The boundaries are composed of a prescribed displacement boundary Γ D and a prescribed traction boundary Γ N and mutually disjointed as Γ = ΓD ∪ΓN

and

Γ D ∩ Γ N = φ.

(8)

The body is subjected to the body force intensity b and the prescribed traction t on Γ N . N is an outward unit vector normal to the boundary. Using the principle of virtual work, an equilibrium equation is expressed as a(z, z¯ ) = ℓ(¯z),

∀¯z ∈ Z¯ ,

(9)

where the bilinear strain energy and linear load forms are defined, respectively, as   a(z, z¯ ) = ci jkl εi j (z)εkl (¯z)dΩ = ci jkl z i, j z¯ k,l dΩ Ω

Ω

(10)

474

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Fig. 2. Plane elasticity problems.

and ℓ(¯z) =



 bi z¯ i dΩ + Ω

ΓN

ti z¯ i dΓ .

(11)

Note that (•),i represents the covariant differentiation with respect to the coordinate i. The variational and the solution spaces are respectively defined as Z¯ = {¯z ∈ [H 1 (Ω )]d : z¯ = 0 on Γ D }

(12)

Z = {z ∈ [H 1 (Ω )]d : z = zˆ on Γ D },

(13)

and

where d is the dimension of the given problem and zˆ is a prescribed displacement. Using the isoparametric mapping and the Galerkin method, the response and virtual response in discrete space are expressed, in terms of NURBS basis functions and control points, as zh (4) =

CP 

W I (4)y I

(14)

W I (4)¯y I ,

(15)

I

and z¯ h (4) =

CP  I

where zh ∈ Z h ⊂ Z and z¯ h ∈ Z¯ h ⊂ Z¯ . Note that the control point B I and the response coefficient y I at control point are defined in physical coordinates. For brevity of notation for discrete response and function space, z and Z will be used instead of zh and Z h , hereafter. Also, same notation is applicable for the virtual ones. Using the isogeometric discretization of Eqs. (14) and (15), the variational equation (9) can be rewritten as a(z, z¯ ) = ℓ(¯z),

∀ z¯ =

CP 

W J (4)¯y J ∈ Z¯ ,

(16)

J

where a(z, z¯ ) =

  CP Ω I,K

Ci jkl W I, j W K ,l yi I y¯k K dΩ

(17)

and ℓ(¯z) =

  CP Ω

I

 bi W I y¯i I dΩ +

ΓN

BCP 

ti W˜ I y¯i I dΓ .

(18)

I

The BCP and W˜ denote the number of control points on the boundary and the one dimensional NURBS basis function on the boundary, respectively. Note that W I, j represents the covariant differentiation of the I-NURBS basis function

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) Contour integral.

475

(b) Local coordinate system. Fig. 3. Contour integral around crack tip.

defined in the parametric space with respect to the coordinate j in the physical domain Ω . The convergence, accuracy, and computing cost of isogeometric analysis compared to classical finite element analysis were discussed in a very clear and complete manner in Hughes et al. [13]. 3.2. Fracture mechanics Interaction integral formulation for SIF The fracture parameter such as the direction of crack growth is represented in terms of the SIFs that can be obtained in several ways; virtual crack extension method, J -integral method, crack opening displacement method, and so on (see, e.g. [2,27]). In this paper, the SIFs are obtained using an interaction integral (M-integral) derived from the J -integral (see, e.g. [3]). For a combined problem of mode I (opening) and mode II (in-plane shear) based on linear elasticity, we consider two equilibrium states; (1) an actual state under given loading and boundary conditions and (2) an auxiliary state corresponding to unit SIF. Fig. 3(a) shows the domain Ω and boundary Γb of the interaction integral (M-integral). The closed boundary Γ1 around the crack-tip approaches zero as the perimeter of Γ1 is reduced. In order to precisely evaluate the auxiliary stress and displacement fields, it is crucial to exactly represent the geometry of the crack-face as well as the local coordinate system in Fig. 3(b). At a new equilibrium state (superscript s) obtained from the superposition of the actual state (superscript 1) and auxiliary state (superscript 2), the strain energy density W (s) at the state (s) is given by   1  (1) (2) (1) (2) W (s) ≡ W (1) + W (2) + W (1,2) = σi j + σi j εi j + εi j 2  1 (1) (1) 1 (2) (2) 1  (1) (2) (2) (1) (19) = σi j εi j + σi j εi j + σi j εi j + σi j εi j . 2 2 2 The J -integral can be defined as       ∂ z i(1) + z i(2) W (S) δ1 j − σ (1) + σ (2)  n j dΓ J (s) ≡ J (1) + J (2) + M (1,2) = ij ij ∂ x1 Γb 



W

= Γb

(1)

(1) (1) ∂z i δ1 j − σi j ∂ x1



 + W

(2)

(2) (2) ∂z i δ1 j − σi j ∂ x1



n j dΓ + M (1,2)

(20)

where M (1,2) is called the interaction integral (M-integral). Since the relation between the J -integral and the SIFs is given by      1  (1) (2) 2 (1) (2) 2 (s) J = KI + KI + K II + K II , (21) E˜

476

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

the M-integral can be derived as  2  (1) (2) (1) (2) K I K I + K II K II , M (1,2) = E˜

(22)

where E˜ takes E for plane stress problems and E/(1 − v 2 ) for plane strain ones. In Eqs. (21) and (22), the SIFs associated with the actual and auxiliary fields are denoted by the superscripts 1 and 2, respectively. The displacement and stress fields for the actual mixed mode state (superscript 1) can be obtained from either FEA or IGA. In the auxiliary state (superscript 2) of mode I, the displacement and the stress corresponding to a unit SIF are respectively given by Refs. [2,27].     θ  θ 2   (2)   cos  κ − 1 + 2 sin  z1 1 r  2 2    = (23) (2) θ  2µ 2π  z2 2 θ   sin  κ + 1 − 2 cos 2 2 and      θ 3θ   sin  1 − sin    (2)   2 2              σ11   1 θ θ 3θ  (2) , cos 1 + sin sin σ22 = √   2  2 2  2πr       (2)         σ12       sin θ cos 3θ 2 2

(24)

where κ = 3 − 4v for plane strain problems and κ = (3 − v)/(1 + v) for plane stress ones. µ is the shear modulus. Also, in the auxiliary state (superscript 2) of mode II, the displacement and the stress corresponding to a unit SIF are respectively given by      θ 2 θ    (2)     sin κ + 1 + 2 cos  z1 1 r  2 2    = (25) (2) θ  θ 2µ 2π  z2   cos  −κ + 1 + 2 sin2 2 2 and      θ   θ 3θ   sin 2 + cos cos    (2)      2 2 2   σ            11    1 θ 3θ θ (2) sin sin . cos σ22 = √    2 2 2 2πr   (2)               σ12   θ 3θ    cos θ  1 − sin sin 2 2 2 (2)

Taking K I

(1)

KI

(26)

(2)

= 1, K II = 0 in Eq. (22), the SIF for mode I can be obtained as

=

E˜ (1,2) M . 2 (2)

Similarly, taking K I

(27) (2)

= 0, K II = 1, the SIF for mode II can be obtained as

E˜ (1,2) M . 2 Generally, the M-integral in Eq. (20) is written as (1)

K II =

(1,2)

M (1,2) = MΩ

(1,2)

+ M Γc ,

(28)

(29)

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

477

(1)

where, under the traction-free condition σi j n j = 0 on the crack-face (see [12]), (1,2)

M Γc

 =



 (2) (1) (2) (1) σi j εi j n 1 − σi j z i,1 n j qdΓ



 (2) (1) (2) (1) (2) (1) (2) (1) σ21 z 1,2 n 1 + σ22 z 2,2 n 1 − σ12 z 1,1 n 2 − σ22 z 2,1 n 2 qdΓ .

Γc+ ∪Γc−

 = Γc+ ∪Γc−

(30)

n 1 and n 2 represent the x and y components of the normal vector in the local coordinate system at crack-tip, (2) (2) (1,2) respectively. MΓc vanishes for straight crack problems since n 1 , σ12 , and σ22 in the auxiliary state vanish on the crack-face but does not for general curved ones. The domain integral in the interaction integral formulation is written as    (2) (1) ∂q ∂q ∂q (1,2) (1) ∂z i (2) ∂z i MΩ = σi j + σi j − W (1,2) dΩ , (31) ∂ x1 ∂ x j ∂ x1 ∂ x j ∂ x1 Ω where q is a smooth weight function and takes the maximum of one around the crack tip and zero at the Γ2 boundary. In the isogeometric approach, using the NURBS in Fig. 4(a), a smooth q-function field can be constructed in the parametric space as shown in Fig. 4(b). q(4) =

CP 

W I (4)Q I ,

(32)

I

where W I and Q I are the NURBS basis function and the corresponding coefficient, respectively. If we take the coefficients Q I = 1 for the red-colored basis functions and Q I = 0 for the others as shown in Fig. 4(a), the resulting q-function field looks like a plateau in Fig. 4(b). In the M-integral domain, except the region of C 0 -continuity due to the knot repetition and the patch division, C 1 -continuity is obtained, which is advantageous to obtain a smooth q-function field. The domain integral of M-integral in Eq. (31) can be explicitly rewritten as    (2) (2) (2) (2) (1) ∂z 1 ∂q (1,2) (1) ∂z 1 ∂q (1) ∂z 2 ∂q (1) ∂z 2 ∂q σ11 MΩ = + σ12 + σ21 + σ22 dΩ ∂ x1 ∂ x1 ∂ x1 ∂ x2 ∂ x1 ∂ x1 ∂ x1 ∂ x2 Ω    (1) (1) (1) (1) (2) ∂z 1 ∂q (2) ∂z 2 ∂q (2) ∂z 1 ∂q (2) ∂z 2 ∂q σ12 + + σ22 − σ12 − σ22 dΩ ∂ x1 ∂ x2 ∂ x1 ∂ x2 ∂ x2 ∂ x1 ∂ x2 ∂ x1 Ω  (33) ≡ [h 1 + h 2 + h 3 + h 4 + h 5 + h 6 − h 7 − h 8 ] dΩ . Ω

Note that ∂q/∂ x j represents the covariant differentiation of the q-function defined in the parametric space with respect to the local coordinate j in the physical domain Ω . Combining Eqs. (30) and (33), the interaction integral of Eq. (29) is finally obtained, in an isogeometric discrete form, as M

(1,2)

=

  CP  Ω I,J

+

    (2) (2) (2) (2) z 1,1 C11kl + z 2,1 C21kl W J,1 + z 1,1 C12kl + z 2,1 C22kl W J,2 W I,l yk I Q J dΩ

  CP 

W I,1 W J,2 − W I,2 W J,1



Ω I,J

 +

BCP 

Γc+ ∪Γc− I,J

n 1 W˜ I,2 − n 2 W˜ I,1

 (2) (2) σ12 y1I + σ22 y2I Q J dΩ

  (2) (2) σ21 y1I + σ22 y2I W˜ J Q J dΓ .

Notice that the boundary integral and the normal vector can be precisely evaluated in the isogeometric method.

(34)

478

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) NURBS basis functions.

(b) Resulting q-function field.

Fig. 4. Representation of q-function using NURBS. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4. Isogeometric shape design sensitivity analysis of SIFs 4.1. Material derivatives Consider the variation of domain from an original domain Ω to a perturbed one Ωτ as shown in Fig. 5. Suppose that only one parameter τ defines a transformation T . The mapping T : X → Xτ , X ∈ Ω is given by Xτ ≡ T (X, τ )

(35)

Ωτ ≡ T (Ω , τ ).

(36)

and

A design velocity field that is equivalent to a mapping rate can be defined as V(Xτ , τ ) ≡

dXτ dT (X, τ ) ∂ T (X, τ ) = = . dτ dτ ∂τ

The material derivative of response z at X ∈ Ω is defined as   d z˙ ≡ zτ (X + τ V(X)) = z′ + ∇zT V, dτ τ =0

(37)

(38)

where z′ and ∇z are the partial derivative and the gradient of z, respectively. Consider the performance measures in domain and boundary integral forms,  ψ1 = h(X)dΩ (39) Ω

and ψ2 =

 Γ

g(X)dΓ .

The first order variations of Eqs. (39) and (40) with respect to the shape design parameter τ are derived as    ψ1′ ≡ h ′ (X) + ∇h(X)T V(X) + h(X)div V(X) dΩ Ω

(40)

(41)

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

479

Fig. 5. Variation of domain.

and ψ2′ ≡

  Γ

  g ′ (X) + g(X) div V(X) − N(X)T ∇V(X)N(X) dΓ ,

(42)

where N(X) denote a normal vector. For more detailed derivation of shape sensitivity expressions, see Haug et al. [28]. 4.2. Isogeometric shape DSA Taking the first order variation of the bilinear strain energy form of Eq. (10) and the linear load form of Eq. (11) with respect to the shape design parameter τ , we have the following: ′ [a(z, z¯ )]′ = aex (z, z¯ ) + a(˙z − ∇zT V, z¯ ) + a(z, z˙¯ − ∇ z¯ T V) ′ = aex (z, z¯ ) + a(˙z, z¯ ) − a(∇zT V, z¯ ) + a(z, z˙¯ ) − a(z, ∇ z¯ T V)

(43)

and [ℓ(¯z)]′ = ℓ′ex (¯z) + ℓ(z˙¯ − ∇ z¯ T V) = ℓ′ex (¯z) + ℓ(z˙¯ ) − ℓ(∇ z¯ T V),

(44)

′ (z, z ¯ ) and ℓ′ex (¯z) denote the explicit variation terms with the dependence of their arguments on the shape where aex design parameter suppressed. Also, the linearity of strain energy and load forms and the fact that z′ = z˙ − ∇zT V and z¯ ′ = z˙¯ − ∇ z¯ T V are used. Assuming that external loads are independent of shape variations, b′ = 0 and t′ = 0. Using the fact that [a(z, z¯ )]′ = [ℓ(¯z)]′ , ∀¯z ∈ Z¯ , the shape sensitivity equation is finally derived as

a(˙z, z¯ ) = ℓ′V (¯z) − aV′ (z, z¯ ) + ℓ(z˙¯ ) − a(z, z˙¯ ),

∀¯z ∈ Z¯ ,

(45)

where a(˙z, z¯ ) = ℓ′V (¯z) ≡



ci jkl z˙ i, j z¯ k,l dΩ , Ω ℓ′ex (¯z) − ℓ(∇ z¯ T V)

= Ω

(46)

(bi,m z¯ i Vm + bi z¯ i Vm,m )dΩ +



(ti,k z¯ i Nk + κti z¯ i )Vm Nm dΓ ,

(47)

ΓN T

′ T aV′ (z, z¯ ) ≡ aex (z, z¯ ) − a(∇z V, z¯ ) − a(z, ∇ z¯ V) = − (ci jkl z i,m Vm, j z¯ k,l + ci jkl z i, j z¯ k,m Vm,l − ci jkl z i, j z¯ k,l Vm,m )dΩ , Ω   ℓ(z˙¯ ) = bi z˙¯ i dΩ + ti z˙¯ i dΓ ,

(48) (49)

ΓN

Ω

and a(z, z˙¯ ) =

 Ω

ci jkl εi j (z)εkl (z˙¯ )dΩ =

 Ω

ci jkl z i, j z˙¯ k,l dΩ .

(50)

480

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

For flexibility of design variation, weights as well as positions of control points are employed as design variable. Nagy et al. [29] and Qian [30] conducted design sensitivity analysis in isogeometric shape design optimization through differentiating the discrete matrix equation with respect to position and weight of control points. A design velocity field can be defined as the first variation of position X and written in a discrete form, using the isogeometric discretization of Eq. (14).   CP  d  dXτ   = W I ((w I )τ ; ξ, η)(B I )τ  V(X) ≡   dτ τ =0 dτ I τ =0

=

CP 

CP 

W˙ I (w I ; ξ, η)B I +

I

W I (w I ; ξ, η)VBI ,

(51)

I

where the weight w I and the control point are perturbed as (w I )τ = w I + τ δw I = w I + τ VIw ,

(52)

(B I )τ = B I + τ δB I = B I + τ V IB .

(53)

and

The detailed derivation of W˙ I is found in Appendix D. Considering the design dependence of NURBS functions, the material derivatives of response and virtual response are obtained as     CP CP   · z˙ ≡ y˙ + y = W I (4)˙y I + W˙ I (4)y I (54) I

I

and z˙¯ ≡ y˙¯ + y¯ · =

 CP 

 W I (4)y˙¯ I

+

 CP 

I

 W˙ I (4)¯y I

.

(55)

I

Eqs. (49) and (50) are discretized, using the isogeometric way of Eqs. (54) and (55), as    BCP BCP   p p · ti W˜ I y˙¯ i I dΓ ℓ(z˙¯ ) ≡ L(y˙¯ ) + L(¯y ) = bi W I y˙¯ i I dΩ + Ω

+

 BCP  Ω

ΓN

I

p bi W˙ I y¯i I dΩ

I

 +

BCP 

ΓN

p ti W˙˜ I y¯i I dΓ

I

 (56)

I

and a(z, z˙¯ ) ≡ A(z, y˙¯ ) + A(z, y¯ ) = ·

+

 CP  Ω I,K

 CP  Ω I,K

 p p Ci jkl W I, j W K ,l yi I y˙¯ k K dΩ

 p p Ci jkl W I, j W˙ K ,l yi I y¯k K dΩ

.

(57)

Using the fact that A(¯z, y˙¯ ) = L(y˙¯ ), ∀ y˙¯ ∈ Z , the variational equation (45) is discretized as A(˙y, z¯ ) = L ′V (¯z) − A′V (z, z¯ ) − A(y·, z¯ ) + L(¯y·) − A(z, y¯ ·),

∀ z¯ =

CP 

W I (4)¯y I ∈ Z ,

(58)

I

where A(˙y, z¯ ) =

  CP Ω I,K

p

p

Ci jkl W I, j W K ,l y˙i I y¯k K dΩ ,

(59)

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

L ′V (¯z) =

CP  



Ω I,M

 (bi,m W I W M + bi W I W M,m ) y¯i I VmBM + (bi,m W I W˙ M + bi W I W˙ M,m ) y¯i I Bm M dΩ

BCP 

 +

Γ N I,M

A′V (z, z¯ ) = − ×



481

 (ti,k n k + κti )W˜ I W˜ M y¯i I VmBM Nm + (ti,k n k + κti )W˜ I W˙˜ M y¯i I Bm M Nm dΓ ,

CP 

Ω I,K ,M

(60)

Ci jkl (W I,m W M, j W K ,l + W I, j W K ,m W M,l − W I, j W K ,l W M,m ) CP 



yi I y¯k K VmBM dΩ

−

Ω I,K ,M

Ci jkl (W I,m W˙ M, j W K ,l + W I, j W K ,m W˙ M,l

− W I, j W K ,l W˙ M,m )yi I y¯k K Bm M dΩ ,

(61)

and A(y·, z¯ ) =

  CP Ω I,K

p

p

Ci jkl W˙ I, j W K ,l yi I y¯k K dΩ .

(62)

Higher order geometric information like the normal vector and curvature in the boundary integrals in Eq. (60) are crucial for the accurate evaluation of shape design sensitivity. Contrary to the finite element DSA that uses approximate geometry, the isogeometric DSA uses exact geometry, which eventually leads to precise results of shape DSA (see [16]). The isogeometric shape design sensitivity provides the identical results to the finite element one as long as all design boundaries have linear geometry and the domain is perturbed under linear velocity field. For general velocity fields, however, the isogeometric approach prevents the loss of higher order geometric information so that precise sensitivity is expected and could be used in shape design optimization. 4.3. Isogeometric shape DSA of SIFs The M-integrals in Eqs. (30) and (31) include the derivative quantities with respect to the local coordinates at the crack-tip. Thus, it is necessary to derive the first variations of local coordinates, q-function and its derivative, and tangential angle θ c before deriving the shape sensitivity of SIFs. Local coordinates The auxiliary fields of σ(2) (x, y) and z(2) (x, y) are defined on the local coordinate system at the crack-tip as shown in Fig. 3. The local coordinate system is also dependent on the shape variations. Using the tangential angle θ c and global coordinate Xc at the crack-tip, an arbitrary point x is written, in the local coordinates, as x = T(X − Xc ),

(63)

where X is the global coordinate of arbitrary point x and the transformation is given by   cos θ c sin θ c . T= − sin θ c cos θ c

(64)

Taking the material derivative of Eq. (63) leads to ˙ − Xc ) + T(X ˙ −X ˙ c) = V ˜ θ x + v − vc , x˙ = T(X

(65)

where  − sin θ c T˙ = θ˙ c − cos θ c

  cos θ c 0 c = − sin θ −θ˙ c

 θ˙ c ˜ θ T, T≡V 0

(66)

and ˙ −X ˙ c ). v − vc = T(X

(67)

482

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Fig. 6. Variation of the crack-tip local coordinate system.

Using Eq. (65) and the fact that x˙ = x′ + v, the first variation of local coordinates is obtained by ˜ θ x − vc . x′ = V

(68)

θc

Tangential angle A local Cartesian coordinate system can be defined using the tangential and normal directions of crack-face curve at the crack-tip. Also, all the physical quantities in the interaction integral are described in the local coordinate system whereas the elasticity analysis gives the solution in the global coordinate system. Therefore, the physical quantities in global coordinates should be transformed into those in the local coordinates. Also, it is crucial to consider the change of the local coordinate system in a perturbed design as shown in Fig. 6. Γ and Γτ represent the crack-face curves in original and perturbed designs, respectively. Actually, the crack orientation could drastically change during crack propagation. However, in this paper, we are not trying to consider such drastic shape variation of crack but to develop a precise numerical method to obtain the design sensitivity of SIFs in the curved shape variations of crack. Thus, it is assumed that the crack-face curve Γτ is regular at the position of original crack-tip Xc ; a sufficiently smooth shape variation is considered in the perturbed design. Then, Γτ can be locally represented by the graph of a regular function f , say X 2 = f (X 1 ). Also, the tangential angle θ = θ (X) does not depend on the shape parameter τ explicitly but on the position X = (X 1 , f (X 1 )). Its material derivative is obtained by θ (Xcτ ) − θ (Xc ) = θ,X 1 (Xc )V1 (Xc ), τ →0 τ

θ˙ c ≡ lim

(69)

where Xcτ = Xc + τ V(Xc ).

(70)

The material derivative of the tangential angle is also rewritten as (see Appendix A for details)   θ˙ c = κ c sec θ c V1 (Xc ) = κ c vc  , where ∥Vc ∥

∥vc ∥ is used. Vc

V(Xc ) is the global design velocity vector at the crack-tip Xc

(71) and vc

v(Xc ) is the

= ≡ ≡ local design velocity vector. For the exact evaluation of κ c , it is necessary to describe the exact geometry of the crack near the crack-tip. Due to the approximate representation of geometry in the finite element analysis, sufficient mesh refinement near the crack-tip is required to accurately calculate κ c . In the isogeometric DSA, however, it is possible to obtain exact κ c in a coarse mesh since the geometry is exactly described. Therefore, the isogeometric DSA can give accurate sensitivity of SIFs even in coarse mesh. Also, the calculation of the curvature requires at least quadratic basis function, which means that Eq. (71) vanishes if linear basis functions are used. In other words, the linear FEA is not able to consider the effects of configuration design variation.

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

483

q-function The derivative of the q-function with respect to the local coordinate system is present in the M-integral. The first variation of the q-function is obtained by ¯ T v, q ′ = q˙ − ∇q T V = q˙ − ∇q

(72)

where ∇¯ and ∇ denote the gradient operators for the differentiations with respect to the local and global coordinates, respectively. Using the coordinate transformation of Eq. (64), the derivative of the q-function is written as ¯ = T∇q. ∇q

(73)

Taking the material derivative of Eq. (73) and using Eq. (66) leads to ˙ ˜ θ ∇q ¯ · = T∇q ¯ + ∇¯ q˙ − ∇q ¯ T ∇v. ¯ (∇q) + T(∇q)· = V

(74)

The first variation of Eq. (73) is derived, using Eq. (72), as ¯ ′ = (∇q) ¯ · − ∇( ¯ ∇q) ¯ Tv (∇q) ˜ θ ∇q ¯ + ∇¯ q˙ − ∇q ¯ T ∇v ¯ − ∇( ¯ ∇q) ¯ T v. =V

(75)

Interaction integral The material derivative of the M-integral is derived as (1,2) (1,2) M˙ (1,2) = M˙ Ω + M˙ Γc   ′  h (X) + div (h(X)V(X)) dΩ = Ω     + g(X) ˙ + g(X) div V(X) − N(X)T ∇V(X)N(X) dΓ ,

(76)

h = h1 + h2 + h3 + h4 + h5 + h6 − h7 − h8

(77)

  (2) (1) (2) (1) (2) (1) (2) (1) g = σ21 z 1,2 n 1 + σ22 z 2,2 n 1 − σ12 z 1,1 n 2 − σ22 z 2,1 n 2 q.

(78)

Γc+ ∪Γc−

where

and

In general curved crack problems, h i′ = Hi + H¯ i + Hiθ ,

(i = 1–8),

(79)

where Hi is the sensitivity considering the shape dependency of actual field only, which is suggested by Chen et al. [5] for straight crack problems. H¯ i is the shape dependency of auxiliary field and q-function, which is necessary for the curved crack problems with a design velocity field varying with position in the M-integral domain. For the case of i = 1,   (2) (2) (2) ′ (1) (2) (1) . (80) H¯ 1 = σ11 (z 1,1 q,1 )′ = σ11 (z 1,11 x1′ + z 1,12 x2′ )q,1 + z 1,1 q,1 Also, Hiθ represents the effects of configuration design variation to the actual field.  (1) (2) E  ′ ′ ′ T1I T1J + T1I T1J + ν(T2I′ T2J + T2I T2J ) z I,J z 1,1 q,1 2 1−ν   E (1) (1) (2) = (z 1,2 + z 2,1 )z 1,1 q,1 κ c vc  . 1+ν

H1θ =

(81)

Detailed derivations for the other components are available in Appendix B. In the previous study [6], H¯ i is not considered in the straight crack problems since a constant velocity field is used in the M-integral domain. Also, Hiθ

484

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

does not appear in the straight crack problem as no configuration design variations occur. The curvature of crack-face curve κ c is included in the H¯ i and Hiθ (i = 1–8), which means that the exact representation of the crack geometry near the crack-tip is significant for accurate shape design sensitivity of the curved crack problem. The material derivative of the integrand function g in Eq. (76) is obtained by (2) (1)

(2) (1)

(2) (1)

(2) (1)

(2) (1)

(2) (1)

˙ g˙ = (σ˙ 21 z 1,2 + σ21 z˙ 1,2 + σ˙ 22 z 2,2 + σ22 z˙ 2,2 )n 1 q + (σ21 z 1,2 + σ22 z 2,2 )(n˙ 1 q + n 1 q) (2) (1)

(2) (1)

(2) (1)

(2) (1)

(2) (1)

(2) (1)

˙ − (σ˙ 12 z 1,1 + σ12 z˙ 1,1 + σ˙ 22 z 2,1 + σ22 z˙ 2,1 )n 2 q − (σ12 z 1,1 + σ22 z 2,1 )(n˙ 2 q + n 2 q),

(82)

where (2)

(2)′

σ˙ i j = σi j

(2)

(2)

(2)

(2)

+ σi j,k vk = σi j,1 x1′ + σi j,2 x2′ + σi j,k vk

(83)

and (1) (1) (1)′ (1) (1) z˙ i, j = T˙i I T j J z I,J + Ti I T˙ j J z I,J + Ti I T j J z I,J + z i, jk Vk (1) (1) (1) (1) (1) (1) = T˙i I T j J z I,J + Ti I T˙ j J z I,J + z˙ i, j − z i,k j Vk − z i,k Vk, j + z i, jk Vk (1) (1) (1) (1) = T˙i I z I, j + z i,J T˙ j J + z˙ i, j − z i,k Vk, j .

(84)

Using a coordinate transformation method, the material derivative of a normal vector in the local coordinate system at crack-tip is represented as   n˙ i = T˙i I N I + Ti I N˙ I = T˙i I N I + n j v j,k n k n i − vk,i n k , (85) where the following (see [28]) is used.     N˙ I = N J V J,K N K N I − VK ,I N K = n j v j,k n k N I − vk,I n k .

(86)

Detailed derivations for the other components are available in Appendix C. 5. Numerical examples 5.1. Inclined crack problems For the verification of the SIF analysis using the isogeometric approach, the analysis results for the inclined crack problem in Fig. 7(a) are compared with those by G. Chen et al. [5]. Consider a plate with an inclined crack under a remote tension of σ ∞ = 1, and the angle of inclination of the crack is β = 45◦ . The plate has a width (2W ) of 20 and a length (2L) of 20. Two cases are considered; crack sizes (a/W ) of 0.1 and 0.05. An h-refinement scheme is employed around the inclined crack for the numerical results. The exact solution of K I and K II for the inclined crack in an infinite plate (a/W → 0) is given as √ K I = σ ∞ sin2 β πa (87) and √ K II = σ ∞ sin β cos β πa.

(88)

In the isogeometric approach, quadratic NURBS basis functions are used. For the 45◦ inclined crack problem, modeI SIF (K I ) and mode-II SIF (K II ) are computed from the FEA and the IGA using quadratic basis functions. The computed SIFs are compared with exact ones in Table 1. The IGA (1308 DOFs) gives more accurate results than the conventional FEA (2800 DOFs), even with the model of less number of DOFs, for all ratios a/W . Fig. 7(b) and (c) show the contours of realistic design velocity fields for a/W = 0.1 and a/W = 0.05, respectively, obtained by solving a displacement–loaded elasticity problem. The amount of design perturbation δ takes 1% of initial crack length (δ = 0.01a). The exact sensitivities (a) are obtained by differentiating Eqs. (87) and (88) with respect to the design variable a, and the IGA results are compared with those from the literature [5] and finite differencing as varying the ratio a/W . The results of finite differencing are not available in the literature [5], and the IGA results

485

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) Design perturbation δ.

(b) a/W = 0.1.

(c) a/W = 0.05.

Fig. 7. Design velocity field in the inclined crack problem. Table 1 Stress intensity factors for 45◦ inclined crack problem. a/W = 0.1 Exact SIF (a)

Analytical SIF (b)

(b)/(a) × 100%

a/W = 0.05 Exact SIF (c)

Analytical SIF (d)

(d)/(c) × 100%

KI

Quadratic FEA [5] Quadratic IGA

0.8862

0.8729 0.8923

98.50 100.68

0.6267

0.6298 0.6246

100.50 99.67

K II

Quadratic FEA [5] Quadratic IGA

0.8862

0.8552 0.8881

96.50 100.21

0.6267

0.6091 0.6225

97.20 99.33

Table 2 Sensitivity of SIFs for 45◦ inclined crack problem (a/W = 0.10). (a) Exact sensitivity

(b) Analytical sensitivity

(c) FDM sensitivity

(b)/(a) × 100%

(b)/(c) × 100%

KI

Quadratic FEA [5] Quadratic IGA

0.4431

0.4350 0.4554

N/A 0.4541

98.20 102.76

N/A 100.27

K II

Quadratic FEA [5] Quadratic IGA

0.4431

0.4265 0.4459

N/A 0.4459

96.20 100.64

N/A 100.00

Table 3 Sensitivity of SIFs for 45◦ inclined crack problem (a/W = 0.05). (a) Exact sensitivity

(b) Analytical sensitivity

(c) FDM sensitivity

(b)/(a) × 100%

(b)/(c) × 100%

KI

Quadratic FEA [5] Quadratic IGA

0.6267

0.6297 0.6294

N/A 0.6273

100.49 100.43

N/A 100.33

K II

Quadratic FEA [5] Quadratic IGA

0.6267

0.6078 0.6133

N/A 0.6113

97.00 97.86

N/A 100.31

show good agreement with the results of finite differencing. Compared with the FEA results, the IGA gives better ones in spite of using less than 46% of DOFs in Table 2 and Table 3. The enhanced sensitivity results in the IGA are obtained by using more realistic design velocity field and considering the design dependency of auxiliary field and the q-function. Also, the enhancement comes from the smoother solution and the q-function using the NURBS basis functions.

486

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) Design perturbation.

(b) x-component.

(c) y-component.

Fig. 8. Design velocity field in the parabolic crack problem.

5.2. Parabolic crack Consider a parabolic crack under a remote tension σ ∞ = 1. The plate has a width (2W ) of 20 and a length (2L) of 20. The parabola equation is represented by y=

α (x + a) (x − a) , a

|x| ≤ a.

(89)

The value of α controls the shape of the parabola; as α increases, the parabola becomes sharper. Since the parabolic curve is represented as quadratic polynomials, both quadratic NURBS and Lagrange polynomials could describe the exact geometry of the crack-face curve. However, in general cases like the circular arc, the Lagrange polynomials are not able to describe the exact geometry. The SIF analysis is conducted for various values α = 0.1–1 with a fixed crack-tip position a. Fig. 9(a) shows circular arc is able to be represented by three control points B1 , B2 , B3 and the corresponding weights are w1 = w3 = 1, w2 = sin(γ ). If the weight w2 = 1, then the curve geometry becomes parabolic as in Fig. 9(b), which means the parabolic arc can be represented by a B-spline using the same control net as used in the circular arc. Since the exact solution of SIFs for the parabolic crack problem is unavailable, the computed SIFs are compared with those using a singular integral equation (see [31]). Fig. 10 shows the SIF analysis results of K I and K II for various α values. The numbers of DOFs for the FEA and the IGA are 1360 and 1308, respectively. Both the quadratic FEA and the IGA show good agreement with the reference solution. However, the linear FEA results in significant errors due to the linearly approximated representation of the crack geometry which leads to inaccurate evaluation of auxiliary stress and displacement fields. Detailed discussions are given in the following circular arc crack problem. A shape DSA for the SIFs in the parabolic crack problem is performed. Fig. 8(b) and (c) show the contours of componential design velocity field, respectively, obtained by solving a displacement–loaded elasticity problem. The design perturbation is made such that the shape of the crack maintains parabolic shape as shown in Fig. 8(a). The amount of design perturbation δ takes 0.01% of initial length a (δ = 10−4 a). Since the exact sensitivity of SIFs for the parabolic crack problem is not available in the literature, the IGA results are compared with finite differencing as varying the α. Table 4 shows the DSA results of SIFs when a linear FEA is employed. In the shape sensitivity expressions, the higher order geometric information like the curvature in Eq. (71) requires the second order derivatives and thus cannot be computed by using the linear FEA. In other words, the linear FEA is not able to consider the effects of configuration design variation, which might result in inaccurate sensitivity results.   θ˙ c = κ c sec θ c V1 (Xc ) = κ c vc  . Next, the sensitivity results of SIFs by the quadratic FEA and IGA are shown in Table 5, as varying the α value.

487

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Fig. 9. Circular and parabolic arc curves.

KII

KI

α

α

(a) Mode-I (K I ).

(b) Mode-II (K II ). Fig. 10. SIF analysis results for various α values.

Table 4 DSA results of SIF using linear FEA. α 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

KI (a) FDM variation

(b) Analytical variation

(b)/(a) × 100%

K II (a) FDM variation

(b) Analytical variation

(b)/(a) × 100%

5.3606E−05 3.5547E−05 1.3939E−05 −5.1232E−06 −1.9060E−05 −2.7959E−05 −3.3034E−05 −3.5543E−05 −3.6503E−05 −3.6593E−05

5.8665E−05 5.3230E−05 4.6404E−05 3.9855E−05 3.4330E−05 2.9860E−05 2.6168E−05 2.3011E−05 2.0184E−05 1.7573E−05

109.44 149.75 332.90 −777.93 −180.12 −106.80 −79.22 −64.74 −55.29 −48.02

3.0364E−05 5.3269E−05 6.5639E−05 6.8915E−05 6.6441E−05 6.1246E−05 5.5306E−05 4.9659E−05 4.4734E−05 4.0631E−05

1.6321E−05 2.9149E−05 3.6903E−05 4.0056E−05 4.0084E−05 3.8450E−05 3.6189E−05 3.3894E−05 3.1847E−05 3.0144E−05

53.75 54.72 56.22 58.12 60.33 62.78 65.43 68.25 71.19 74.19

Compared with the quadratic IGA results for K I and K II in Table 5, the quadratic FEA gives comparable results using a similar number of DOFs. As mentioned previously, the parabolic crack can be represented using quadratic Lagrange polynomial basis functions. Both the quadratic NURBS and the Lagrange polynomials could describe the exact geometry of crack-face and enable to compute the κc exactly. However, compared with the linear FEA results in Table 4, the quadratic FEA provides significantly improved results in Table 5, which implies that the effects of configuration design variation are very important in the shape sensitivity analysis of curved crack problems. 5.3. Circular arc crack Consider a circular arc crack in a plate under the remote tension σ ∞ = 1. The SIF analysis is performed for various angle β with the crack-tip position fixed. The plate has a width (2W ) of 20 and a length (2L) of 20. For the circular

488

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Table 5 DSA results of SIFs for various α. α

Quadratic FEA (a) FDM variation

(b) Analytical variation

(b)/(a) × 100%

Quadratic IGA (a) FDM variation

(b) Analytical variation

(b)/(a) × 100%

KI

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

5.1835E−05 2.7813E−05 2.0432E−06 −1.7625E−05 −2.9645E−05 −3.5704E−05 −3.8056E−05 −3.8368E−05 −3.7705E−05 −3.6652E−05

5.1838E−05 2.7820E−05 2.0589E−06 −1.7594E−05 −2.9596E−05 −3.5637E−05 −3.7974E−05 −3.8276E−05 −3.7609E−05 −3.6555E−05

100.01 100.02 100.77 99.83 99.83 99.81 99.78 99.76 99.74 99.73

5.3042E−05 2.8465E−05 2.1982E−06 −1.7774E−05 −2.9955E−05 −3.6121E−05 −3.8553E−05 −3.8959E−05 −3.8384E−05 −3.7404E−05

5.3045E−05 2.8472E−05 2.2138E−06 −1.7744E−05 −2.9906E−05 −3.6055E−05 −3.8473E−05 −3.8871E−05 −3.8292E−05 −3.7311E−05

100.00 100.02 100.71 99.83 99.84 99.82 99.79 99.77 99.76 99.75

K II

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

3.5205E−05 5.9264E−05 6.9263E−05 6.9025E−05 6.3632E−05 5.6623E−05 4.9811E−05 4.3924E−05 3.9110E−05 3.5276E−05

3.5203E−05 5.9246E−05 6.9213E−05 6.8929E−05 6.3488E−05 5.6435E−05 4.9588E−05 4.3677E−05 3.8849E−05 3.5010E−05

99.99 99.97 99.93 99.86 99.77 99.67 99.55 99.44 99.33 99.24

3.6000E−05 6.0532E−05 7.0635E−05 7.0297E−05 6.4757E−05 5.7621E−05 5.0715E−05 4.4741E−05 3.9843E−05 3.5929E−05

3.5997E−05 6.0514E−05 7.0583E−05 7.0199E−05 6.4610E−05 5.7430E−05 5.0490E−05 4.4494E−05 3.9584E−05 3.5666E−05

99.99 99.97 99.93 99.86 99.77 99.67 99.56 99.45 99.35 99.27

(b) β = 30◦ .

(a) Design perturbation.

(c) β = 60◦ .

Fig. 11. Design velocity field in the circular arc crack problem.

arc crack in an infinite plate (a/W → 0), the exact solutions of K I and K II are given as    2 2 (β/2) cos (β/2) 1 − sin cos (β/2) σ∞ KI = + cos (3β/2) (π R sin (β))1/2 2 1 + sin2 (β/2)

(90)

and σ∞ K II = (π R sin (β))1/2 2



 1 − sin2 (β/2) cos2 (β/2) sin (β/2) 1 + sin2 (β/2)

 + sin (3β/2) ,

where β and R are the angle and the radius of circular arc as shown in Fig. 11(a).

(91)

489

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

a

b

Fig. 12. Plot of mesh in IGA for the curved crack problem. (a) Control net. (b) Knot mesh mapped into physical space.

Table 6 Comparison of computed SIFs. β

Exact solution (a)

Linear FEA (b)

Linear FEA∗ (c)

Quadratic FEA (d)

Quadratic IGA (e)

KI

15◦ 30◦ 45◦ 60◦

1.1796 0.9750 0.6817 0.3528

1.1787 (99.93%) 1.0357 (106.2%) 0.8213 (120.5%) 0.5652 (160.2%)

1.1597 (98.31%) 0.9646 (98.94%) 0.6785 (99.53%) 0.3489 (98.91%)

1.1708 (99.26%) 0.9714 (99.64%) 0.6845 (100.4%) 0.3625 (102.8%)

1.1773 (99.80%) 0.9745 (99.95%) 0.6812 (99.93%) 0.3494 (99.05%)

K II

15◦ 30◦ 45◦ 60◦

0.3189 0.5856 0.7620 0.8303

0.2719 (85.27%) 0.5109 (87.25%) 0.6902 (90.58%) 0.7927 (95.47%)

0.3087 (96.81%) 0.5734 (97.92%) 0.7590 (99.60%) 0.8451 (101.8%)

0.3158 (99.03%) 0.5811 (99.23%) 0.7581 (99.49%) 0.8289 (99.83%)

0.3192 (100.1%) 0.5872 (100.3%) 0.7652 (100.4%) 0.8341 (100.5%)

The SIF analysis results by the IGA using quadratic NURBS basis functions are compared with those by the conventional FEA using quadratic and linear basis functions. Fig. 12 shows that the geometry of crack is exactly represented by the NURBS in the IGA, unlike the conventional FEA. The SIF analyses are conducted as increasing the angle β from 15◦ to 60◦ . Table 6 compares the performance of the IGA (1308 DOFs) and the conventional FEA (1360 DOFs) in the SIF analysis for the circular arc crack problem. Almost the same numbers of DOFs are used for both the IGA and the FEA. As the angle β increases, significant errors occur in the linear FEA (b). As mentioned in the numerical example of the parabolic crack problem, it turns out that the linearly approximated crack geometry results in an inaccurately defined local coordinate system at the crack-tip, which leads to inaccurate evaluation of auxiliary stress and displacement fields and consequently results in significant errors in the SIF analysis. Fig. 13 illustrates the possibility of an inaccurate tangential direction if a linearly approximated geometry is used in the FEA. At the crack-tip, the exact tangential line should be the red one by the regular extension of the crack line. However, the dotted blue line is used due to the mathematical modeling of the crack line in the FEA. The tangential line is very important to determine a local coordinate system at the crack-tip. This explanation could be validated through the SIF analysis results using the linear FEA provided that the local coordinates at the crack-tip is corrected (linear FEA∗ (c)) through using exact tangential angle. The results are significantly improved through the correction of the local coordinate system. The results of quadratic IGA (e) are better than any other results. The IGA is much more advantageous than the FEA due to its natural representation of exact geometry.

490

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Fig. 13. Tangential line in the finite element model.

Table 7 Comparison of design sensitivity of SIFs for various analysis methods. β

Exact variation (a)

Analytical variation (b)

FDM variation (c)

(b)/(a) × 100%

(b)/(c) × 100%

15◦

Linear FEA Quadratic FEA Quadratic IGA

4.3174E−04

5.5820E−04 4.2910E−04 4.3772E−04

4.6894E−04 4.2785E−04 4.3744E−04

129.29 99.39 101.39

119.03 100.29 100.07

30◦

Linear FEA Quadratic FEA Quadratic IGA

−7.2792E−05

4.2941E−04 −5.5545E−05 −7.4043E−05

9.4536E−05 −6.9596E−05 −7.4601E−05

−589.91 76.31 101.72

454.23 79.81 99.25

45◦

Linear FEA Quadratic FEA Quadratic IGA

−6.9241E−04

2.6963E−04 −6.2372E−04 −7.0847E−04

−4.0869E−04 −6.8168E−04 −7.0931E−04

−38.94 90.08 102.32

−65.97 91.50 99.88

60◦

Linear FEA Quadratic FEA Quadratic IGA

−1.2108E−03

1.3170E−04 −1.0508E−03 −1.2460E−03

−9.0693E−04 −1.1922E−03 −1.2471E−03

−10.88 86.78 102.91

−14.52 88.13 99.91

15◦

Linear FEM Quadratic FEM Quadratic IGA

4.5650E−04

2.0807E−04 4.5348E−04 4.6560E−04

3.9332E−04 4.5516E−04 4.6567E−04

45.58 99.34 101.99

52.90 99.63 99.99

30◦

Linear FEM Quadratic FEM Quadratic IGA

7.1809E−04

3.3659E−04 7.0709E−04 7.3334E−04

6.5798E−04 7.1758E−04 7.3357E−04

46.87 98.47 102.12

51.15 98.54 99.97

45◦

Linear FEM Quadratic FEM Quadratic IGA

6.6987E−04

3.3051E−04 6.4984E−04 6.8413E−04

7.0314E−04 6.7245E−04 6.8458E−04

49.34 97.01 102.13

47.01 96.64 99.93

60◦

Linear FEM Quadratic FEM Quadratic IGA

3.0494E−04

1.6935E−04 2.9116E−04 3.0921E−04

4.9495E−04 3.1464E−04 3.0944E−04

55.53 95.48 101.40

34.21 92.54 99.93

KI

K II

A shape DSA for the SIFs in the circular arc crack problem is performed. The angle β is selected as a design variable. Fig. 11(b) and (c) show the contours of a design velocity field, respectively, obtained by solving a displacement–loaded elasticity problem. The design perturbation is made such that the shape of the crack maintains circular shape as shown in Fig. 11(a). The amount of design perturbation δ takes 0.1% of initial angle (δ = 10−3 β) and is used as a prescribed displacement. Table 7 shows the DSA results of SIFs by various analysis methods. The exact variations (a) are obtained by differentiating Eqs. (89) and (91) with respect to the design variable β and multiplying the perturbation amount. As expected, the variation by linear FEA yields significant errors in both SIFs for all the cases of angle β. Compared with the linear FEA results, the quadratic FEA provides significantly improved results, which implies that the effects of configuration design variation are very important in the shape sensitivity analysis of the circular arc crack problems. Even though the quadratic FEA provides accurate results for small angle of β = 15◦ , the results become inaccurate as the angle β increases. As the angle β increases, the geometric approximation increases the errors in finite element discretization and consequently κc that is essential for the consideration of configuration design variation. However, the IGA could represent exact geometry even in coarse

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

(a) Model description.

491

(b) Cubic B-spline interpolation. Fig. 14. Sinusoidal wave crack.

Table 8 Tangential angle and curvature of crack-face curve at crack-tip. N Quadratic FEA Quadratic FEA Quadratic IGA Cubic IGA

17 21 15 11

Tangential angle (rad) Exact

1.0039

Analytical 1.0252 1.0180 1.0039 1.0039

Curvature Exact

0

Analytical 0.5211 0.4380 0.0433 0.0000

mesh. For all the angles, the quadratic IGA provides accurate results; less than 3% errors with exact solution and more than 99% agreement with the results of finite differencing. 5.4. Sinusoidal wave crack Consider a plate with a sinusoidal wave crack under a remote tension of σ ∞ = 1 as shown in Fig. 14(a). The plate has a width (2W ) of 20 and a length (2L) of 20. The sinusoidal wave crack is represented by y = 0.25 sin [2π (x + 0.5)] ,

|x| ≤ 0.5.

(92)

The shape of the sinusoidal wave crack is represented in two ways; piecewise Lagrange interpolation for the conventional FEM and B-spline interpolation of given data for the IGA. The B-spline basis function of order p is able to interpolate the data points of sinusoidal wave, satisfying the continuous end derivatives of order ( p − 1). See Piegl and Tiller [25] and Piegl and Tiller [32] for details of global interpolation satisfying the conditions of end derivatives. Fig. 14(b) shows the curve of sinusoidal wave using cubic B-splines, which interpolates 11 data points from Eq. (92) and satisfies the end conditions of the first and second derivatives. Table 8 compares the tangential angles and curvatures at the crack-tip by various interpolation schemes using N data points for Eq. (92). The FEA uses piecewise Lagrange interpolation using 17 and 21 data points and the IGA uses B-spline interpolation using 15 and 11 data points. The quadratic and cubic IGAs respectively satisfy the conditions of specified first and second order end derivatives. As the number of data points was increased, we can observe that the geometric errors are decreasing. When using the piecewise Lagrange polynomials, further refinements are still required to represent the precise geometry around the crack-tip. However, the tangential angle of the first order derivative is exactly satisfied by the quadratic B-splines and the curvature as well as the tangential angle by the cubic B-splines. Thus, the B-spline basis functions enable to

492

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Table 9 SIF and DSA results by various analysis methods. DOF/N

KI

Quadratic FEA Quadratic FEA Quadratic IGA Cubic IGA

1532/17 2420/21 1458/15 1458/11

K II

Quadratic FEA Quadratic FEA Quadratic IGA Cubic IGA

1532/17 2420/21 1458/15 1458/11

Reference SIF [31]

Analytical SIF

FDM variation (a)

DDM variation (b)

(b)/(a) × 100%

0.6943

0.5935 0.6336 0.6705 0.6797

2.4588E−06 −1.3348E−06 −1.4486E−06 −1.6586E−06

−4.6845E−06 −4.9225E−06 −1.8356E−06 −1.6589E−06

−190.52 368.77 126.72 100.01

0.6292

0.6762 0.6449 0.6389 0.6310

−1.4271E−05 5.8027E−06 5.6920E−06 5.2361E−06

−1.1627E−05 7.3230E−06 5.8818E−06 5.2362E−06

81.47 126.20 103.33 100.00

represent the exact geometry around the crack-tip and to consider the exact geometry in the isogeometric analysis. Consequently, more accurate SIF analysis and DSA results can be expected in the isogeometric approach. Table 9 shows the results of SIF and DSA by various analysis methods. Since the exact SIFs for the sinusoidal wave crack problem are unavailable in other publications, the IGA results are compared with those from other publication [31] using the singular integral equation. The DSA results are compared with finite difference ones. The design perturbation is made such that the shape of the crack maintains sinusoidal shape as shown in Fig. 14(a), and the amount of design perturbation δ takes 10−3 % of initial length a (δ = 10−5 a). Compared with the FEA results, the IGA gives better ones in spite of using less number of DOFs in Table 9. Table 9 demonstrates the importance of exact geometry around the crack-tip. For the accurate SIF analysis, the first order derivative condition around the crack-tip is required to determine the exact local coordinate system at the cracktip. Also, for the accurate DSA of SIFs, the second order derivative condition around the crack-tip is additionally required to consider the effect of configuration design variation. In the quadratic FEA, the accuracy of SIF results is generally improved through mesh refinement but the DSA shows very erroneous results and even shows a different sign (−1.1627E−05) in the quadratic FEA using 17 data points even though improved accuracy is observed though refinement. On the other hand, in the IGA, the accuracy of SIF results is significantly improved due to the consideration of exact geometry around the crack-tip. Accurate SIF results are observed in the case of the first order end derivative condition (Quadratic IGA) and accurate DSA results in the case of the second order end derivative condition (Cubic IGA). 6. Conclusions In this paper, an isogeometric DSA method is developed for the curved crack problems. The differences between the curved and the straight crack problems are (1) the shape dependence of auxiliary field and the q-function in the M-integral, (2) the exact determination of the local coordinate system at the crack-tip, and (3) the consideration of crack-face integral. We demonstrate the significance of the local coordinate system when using linearly approximated element geometry in the conventional FEA. Also, the exact representation of geometry using the NURBS basis functions allows the exact determination of the local coordinate system at the crack-tip and the improved evaluation of the auxiliary field. In the design sensitivity expressions of SIFs, we consider the configuration design variation due to the local coordinate system at the crack-tip and show the importance of curvature of the crack-face curve. Even though having C 0 continuity at the crack-tip due to the repeated knots, the higher order NURBS basis function could provide smooth stress, strain, and q-function fields. Through numerical examples, we demonstrate that the developed IGA method could represent the exact geometry of crack around the crack-tip even with coarse control nets, which yields better SIF results than the conventional FEA. Acknowledgments This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (No. 2010-0018282). The support is gratefully acknowledged. The authors would also like to thank Ms. Inyoung Cho at Korea University for editing assistance.

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

493

Appendix A Consider the slope of a curve f that represents a crack-face curve Γ , tan θ = f ,X 1

(A.1)

in the neighborhood of Xc ⊂ Γ . Introducing a parameter ξ for the parameterization X 1 = X 1 (ξ ) and θ = θ (ξ ), the differentiation of both sides of Eq. (A.1) with respect to ξ gives θ,ξ sec2 θ = f ,X 1 X 1 X 1,ξ

(A.2)

    2 θ,ξ sec2 θ = θ,ξ 1 + tan2 θ = θ,ξ 1 + f ,X . 1

(A.3)

and

Thus, the following holds;   2 θ,ξ 1 + f ,X = f ,X 1 X 1 X 1,ξ 1

(A.4)

and θ,ξ =

 f ,X 1 X 1 f ,X 1 X 1 2 X X = 1 + f ,X 1,ξ 1,ξ = κ X 1,ξ sec θ.   2 3/2 1 1 + f ,X 2 1 1 + f ,X 1

(A.5)

At the position of crack-tip Xc , θ,X 1 =

θ,ξ = κ c sec θ c X 1,ξ

(A.6)

where κ c ≡ κ(Xc ) is a curvature of the crack-face curve Γ at the crack-tip Xc . Then, substituting Eq. (A.6) into Eq. (69) gives   ′ θ˙ c = θ c = κ c V1 (X c ) sec θ c = κ c Vc  (A.7) where Vc ≡ V (Xc ) is the design velocity vector at the crack-tip Xc . Since the magnitude of the vector is invariant for a coordinate transformation, the following holds.   θ˙ c = κ c vc  , (A.8) where vc ≡ v(Xc ) is a design velocity vector at the crack-tip Xc in the local coordinate system at the crack-tip. Appendix B Detailed derivations for H¯ i , (i = 1–8) are as follows.    ′ ′  (2) (2)  (1) (2) (1) (2) z 1,11 x1′ + z 1,12 x2′ q,1 + z 1,1 q,1 , H¯ 1 = σ11 z 1,1 q,1 = σ11  ′   ′  (2)  (2) (1) (2) (2) (1) z 1,11 x1′ + z 1,12 x2′ q,2 + z 1,1 q,2 , H¯ 2 = σ12 z 1,1 q,2 = σ12   ′  ′  (2) (2) (2)  (1) (2) (1) z 2,11 x1′ + z 2,12 x2′ q,1 + z 2,1 q,1 , H¯ 3 = σ21 z 2,1 q,1 = σ21   ′  ′  (1) (2) (1) (2) (2) (2)  H¯ 4 = σ22 z 2,1 q,2 = σ22 z 2,11 x1′ + z 2,12 x2′ q,2 + z 2,1 q,2 ,  ′   ′  (1) (2) (1) (2) (2) (2)  H¯ 5 = z 1,1 σ12 q,2 = z 1,1 σ12,1 x1′ + σ12,2 x2′ q,2 + σ12 q,2 ,  ′   ′  (1) (2) (1) (2) (2) (2)  H¯ 6 = z 2,1 σ22 q,2 = z 2,1 σ22,1 x1′ + σ22,2 x2′ q,2 + σ22 q,2 ,

(B.1) (B.2) (B.3) (B.4) (B.5) (B.6)

494

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

 ′   ′  (2)  (2) (1) (2) (1) (2) H¯ 7 = z 1,2 σ12 q,1 = z 1,2 σ12,1 x1′ + σ12,2 x2′ q,1 + σ12 q,1 ,

(B.7)

 ′   ′  (2)  (2) (1) (2) (1) (2) H¯ 8 = z 2,2 σ22 q,1 = z 2,2 σ22,1 x1′ + σ22,2 x2′ q,1 + σ22 q,1 .

(B.8)

and

Detailed derivations for Hiθ , (i = 1–8) are as follows.  (1) (2)  E  ′ ′ ′ z I,J z 1,1 q,1 T1I T1J + T1I T1J + ν T2I′ T2J + T2I T2J 2 1−ν     E (1) (1) (2) = z + z 2,1 z 1,1 q,1 κ c vc  , 1 + ν 1,2  (1) (2)  ′ E ′ ′ z I,J z 1,1 q,2 = T T2J + T1I T2J + T2I′ T1J + T2I T1J 2(1 + ν) 1I     E (1) (1) (2) = z − z 1,1 z 1,1 q,2 κ c vc  , 1 + ν 2,2  ′  (1) (2) E ′ ′ = T1I T2J + T1I T2J + T2I′ T1J + T2I T1J z I,J z 2,1 q,1 2(1 + ν)     E (1) (1) (2) = z 2,2 − z 1,1 z 2,1 q,1 κ c vc  , 1+ν   (1) (2) E  ′ ′ ′ T2I T2J + T2I T2J + ν T1I′ T1J + T1I T1J z I,J z 2,1 q,2 = 2 1−ν     −E (1) (1) (2) z + z 1,2 z 2,1 q,2 κ c vc  , = 1 + ν 2,1      (1) (2)  (2) (1) (1) ′ = σ12 T1I′ T1J + T1I T1J z I,J q,2 = σ12 z 2,1 + z 1,2 q,2 κ c vc  ,      (1) (2)  (2) (1) (1) ′ = σ22 T2I′ T1J + T2I T1J z I,J q,2 = σ22 z 2,2 − z 1,1 q,2 κ c vc  ,      (1) (2)  (2) (1) (1) ′ = σ12 T1I′ T2J + T1I T2J z I,J q,1 = σ12 z 2,2 − z 1,1 q,1 κ c vc  ,

H1θ =

H2θ

H3θ

H4θ

H5θ H6θ H7θ

(B.9)

(B.10)

(B.11)

(B.12) (B.13) (B.14) (B.15)

and      (1) (2)  (2) (1) (1) ′ H8θ = σ22 T2I′ T2J + T2I T2J z I,J q,1 = −σ22 z 1,2 + z 2,1 q,1 κ c vc  .

(B.16)

Appendix C The material derivatives of the actual strains are as follows.         (1) (1) (1) (1) (1) (1) (1) · (1) (1) z 1,1 = T˙1I z I,1 + z 1,J T˙1J + z˙ 1,1 − z 1,k Vk,1 = z 2,1 + z 1,2 κ c vc  + z˙ 1,1 − z 1,k Vk,1 ,         (1) · (1) (1) (1) (1) (1) (1) (1) (1) z 1,2 = T˙1I z I,2 + z 1,J T˙2J + z˙ 1,2 − z 1,k Vk,2 = z 2,2 − z 1,1 κ c vc  + z˙ 1,2 − z 1,k Vk,2 ,         (1) · (1) (1) (1) (1) (1) (1) (1) (1) z 2,1 = T˙2I z I,1 + z 2,J T˙1J + z˙ 2,1 − z 2,k Vk,1 = −z 1,1 + z 2,2 κ c vc  + z˙ 2,1 − z 2,k Vk,1 ,         (1) (1) (1) (1) (1) · (1) (1) (1) (1) z 2,2 = T˙2I z I,2 + z 2,J T˙2J + z˙ 2,2 − z 2,k Vk,2 = −z 1,2 − z 2,1 κ c vc  + z˙ 2,2 − z 2,k Vk,2 . The material derivative of normal vector in crack-tip local coordinate system is represented as       n˙ 1 = T˙1I N I + n j v j,k n k n 1 − vk,1 n k = n 2 κ c vc  + n j v j,k n k n 1 − vk,1 n k ,       n˙ 2 = T˙2I N I + n j v j,k n k n 2 − vk,2 n k = −n 1 κ c vc  + n j v j,k n k n 2 − vk,2 n k .

(C.1) (C.2) (C.3) (C.4)

(C.5) (C.6)

495

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

Appendix D The NURBS basis function in two-dimensional space is defined as  n  m Wi j (ξ, η) = Ni (ξ )M j (η)wi j Nk (ξ )Ml (η)wkl

(D.1)

k=1 l=1

where Ni and M j are  the B-spline basis functions. wi j is the weight corresponding to each control point. The weight wi j is perturbed as wi j τ = wi j + τ δwi j = wi j + τ Viwj and the material derivative of NURBS basis function is obtained by    m n  d w w Nk (ξ )Ml (η)(wkl + τ Vkl ) W˙ i j (ξ, η) = Ni (ξ )M j (η)(wi j + τ Vi j ) dτ k=1 l=1 τ =0  2 n  m n  m     = Ni (ξ )M j (η) Nk (ξ )Ml (η) Viwj wkl − Vklw wi j Nk (ξ )Ml (η)wkl k=1 l=1

Ni =

NT



Viwj W − wi j Vw 

2 NT WM

k=1 l=1



MjM

.

(D.2)

The partial derivatives of R˙ i j with respect to the parametric coordinates are derived as        w W − w Vw M M T V w W − w Vw M M NT WM T + N NT 2N N N N V i i j i,ξ i i j j j ,ξ ,ξ ij ij ∂ ˙ − Wi j =  2  3 T T ∂ξ N WM N WM (D.3) and        2Ni NT Viwj W − wi j Vw M j M NT WM,η Ni NT Viwj W − wi j Vw M j,η M + M j M,η ∂ ˙ − . Wi j =  2  3 ∂η NT WM NT WM (D.4) Note that the operations of partial derivative with respect to the parametric coordinates and the material derivatives are commutative. References [1] K.H. Nam, I.H. Park, S.H. Ko, Patterning by controlled cracking, Nature 485 (2012) 221–224. [2] T.L. Anderson, Fracture Mechanics—Fundamentals and Applications, second ed., CRC press, Florida, 1995. [3] J.F. Yau, S.S. Wang, H.T. Corten, A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity, J. Appl. Mech. 47 (1980) 335–341. [4] G. Chen, S. Rahman, Y.H. Park, Shape sensitivity and reliability analyses in linear-elastic fracture mechanics, Int. J. Fract. 112 (3) (2001) 223–246. [5] G. Chen, S. Rahman, Y.H. Park, Shape sensitivity analysis in mixed-mode fracture mechanics, Comput. Mech. 27 (2001) 282–291. [6] B.N. Rao, S. Rahman, Continuum shape sensitivity analysis of a mixed-mode fracture in functionally graded materials, Comput. Methods Appl. Mech. Engrg. 194 (2005) 1913–1946. [7] R.A. Feijoo, C. Padra, R. Saliva, E. Taroco, M.J. Venere, Shape sensitivity analysis for energy release rate evaluation and its application to the study of three-dimensional cracked bodies, Comput. Methods Appl. Mech. Engrg. 188 (4) (2000) 649–664. [8] E. Taroco, Shape sensitivity analysis in linear elastic fracture mechanics, Comput. Methods Appl. Mech. Engrg. 188 (2000) 697–712. [9] B.N. Rao, S. Rahman, A continuum shape sensitivity method for fracture analysis of isotropic functionally graded materials, Comput. Mech. 38 (2006) 133–150. [10] M. Gosz, B. Moran, An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions, Eng. Fract. Mech. 69 (3) (2002) 299–319. [11] M.C. Walters, G.H. Paulino, R.H. Dodds, Interaction integral procedures for 3-D curved cracks including surface tractions, Eng. Fract. Mech. (2005) 1635–1663.

496

M.-J. Choi, S. Cho / Comput. Methods Appl. Mech. Engrg. 279 (2014) 469–496

[12] H. Yu, L. Wu, L. Guo, H. Wu, S. Du, An interaction integral method for 3D curved cracks in nonhomogeneous materials with complex interfaces, Int. J. Solids Struct. 47 (2010) 2178–2189. [13] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [14] D.J. Benson, Y. Bazilevs, M.C. Hsu, T.J.R. Hughes, Isogeometric shell analysis: the Reissner–Mindlin Shell, Comput. Methods Appl. Mech. Engrg. 199 (2010) 276–289. [15] S.H. Ha, K.K. Choi, S. Cho, Numerical method for shape optimization using T-spline based isogeometric method, Struct. Multidiscip. Optim. 42 (2010) 417–428. [16] S. Cho, S.H. Ha, Isogeometric shape design optimization: exact geometry and enhanced sensitivity, Struct. Multidiscip. Optim. 38 (1) (2009) 53–70. [17] W.A. Wall, M.A. Frenzel, C. Cyron, Isogeometric structural shape optimization, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2976–2988. [18] N.D. Manh, A. Evgrafov, A.R. Gersborg, J. Gravesen, Isogeometric shape optimization of vibrating membranes, Comput. Methods Appl. Mech. Engrg. 200 (2011) 1343–1353. [19] X. Qian, O. Sigmund, Isogeometric shape optimization of photonic crystals via Coons patches, Comput. Methods Appl. Mech. Engrg. 200 (2011) 2237–2255. [20] C.V. Verhoosel, M.A. Scott, R. de Borst, T.J.R. Hughes, An isogeometric approach to cohesive zone modeling, Internat. J. Numer. Methods Engrg. 87 (2011) 336–360. [21] C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, R. de Borst, An isogeometric analysis approach to gradient damage models, Internat. J. Numer. Methods Engrg. 86 (2011) 115–134. [22] E. De Luycker, D.J. Benson, T. Belytschko, Y. Bazilevs, M.C. Hsu, X-FEM in isogeometric analysis for linear fracture mechanics, Internat. J. Numer. Methods Engrg. 87 (2001) 541–565. [23] S.S. Ghorashi, N. Valizadeh, S. Mohammadi, Extended isogeometric analysis for simulation of stationary and propagating cracks, Internat. J. Numer. Methods Engrg. 89 (2012) 1069–1101. [24] D.F. Rogers, An Introduction to NURBS With Historical Perspective, Academic Press, San Diego, CA, 2001. [25] L. Piegl, W. Tiller, The NURBS Book, second ed., in: Monographs in Visual Communication, Springer, New York, 1997. [26] B. Koo, M. Yoon, S. Cho, Isogeometric shape design sensitivity analysis using transformed basis functions for Kronecker delta property, Comput. Methods Appl. Mech. Engrg. 253 (2013) 505–516. [27] E.E. Gdoutos, Fracture Mechanics Criteria and Applications, Kluwer Academic Publishers, 1990. [28] E.J. Haug, K.K. Choi, V. Komkov, Design Sensitivity Analysis of Structural Systems, Academic, New York, 1986. [29] A.P. Nagy, M.M. Abdalla, Z. G¨urdal, Isogeometric sizing and shape optimization of beam structures, Comput. Methods Appl. Mech. Engrg. 199 (2010) 1216–1230. [30] X. Qian, Full analytical sensitivities in NURBS based isogeometric shape optimization, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2059–2071. [31] Y.Z. Chen, Stress intensity factors for curved and kinked cracks in plane extension, Theor. Appl. Fract. Mech. 31 (1999) 223–232. [32] L. Piegl, W. Tiller, Curve interpolation with arbitrary end derivatives, Eng. Comput. 16 (2000) 73–79.