Semi-analytical isogeometric analysis for wrinkling instability of stiff films bonded to cylindrical modulus-graded compliant substrates

Journal Pre-proofs Semi-analytical isogeometric analysis for wrinkling instability of stiff films bonded to cylindrical modulus-graded compliant substrates Chunlei Li, Qiang Han, Zhan Wang PII: DOI: Reference:

S0263-8223(19)32377-3 https://doi.org/10.1016/j.compstruct.2019.111787 COST 111787

To appear in:

Composite Structures

Received Date: Accepted Date:

20 June 2019 9 December 2019

Please cite this article as: Li, C., Han, Q., Wang, Z., Semi-analytical isogeometric analysis for wrinkling instability of stiff films bonded to cylindrical modulus-graded compliant substrates, Composite Structures (2019), doi: https:// doi.org/10.1016/j.compstruct.2019.111787

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

Semi-analytical isogeometric analysis for wrinkling instability of stiff films bonded to cylindrical modulus-graded compliant substrates Chunlei Li, Qiang Han∗, Zhan Wang School of Civil Engineering and Transportation, South China University of Technology, Guangzhou, Guangdong Province 510640, P.R.China

Abstract Surface wrinkles in film-substrate structures have received considerable attention in science and engineering. Understanding the wrinkling instability is the key to envisioned applications of these structures. In this paper, based on Floquet’s principle, an explicit semi-analytical isogeometric analysis method (SIGA) is proposed for predicting wrinkling instability of a stiff film on cylindrical graded substrate. In total Lagrangian framework, the wrinkle configuration is considered as the stressed state relative to the stress-free state. The linearized variational formulation with respect to stress (wrinkle) is derived from the principle of minimum potential energy and parameterized by non-uniform rational basis splines (NURBS). By using the SIGA method, only the radial direction of layered structure is parameterized and the dimension of model is reduced to the one-dimensional. Meanwhile, the analysis of wrinkling instability can be performed by a generalized eigenvalue problem, which explicitly determines critical wrinkling strain and wavelength. Compared to the finite element method, the proposed approach demonstrates distinct superior in computational efficiency and accuracy. Besides, two typical functions are employed to describe material distribution of inhomogeneous substrate. The critical conditions for wrinkling instability of the film-substrate system are computed and the influences of various geometric and material parameters are discussed in detail. Keywords: Semi-analytical isogeometric analysis; wrinkle; film-substrate; graded material; cylindrical structures

1. Introduction Film-substrate systems are ubiquitious in nature and engineering applications [1–3], by reasons of the factors such as the fabrication and use of thin film devices , the growth of biological tissues and the physical stimulation, a compressing stress is produced to be beyond the critical value in the thin film, which may leads to the onset of wrinkling instability. As a typical case of instability, wrinkle is not only a failure mode avoided in most engineering designs [4] but also a useful mechanical ∗ Corresponding

author Email address: [email protected] (Qiang Han)

Preprint submitted to Composite Structures

December 10, 2019

feature for designed performances [5]. Wrinkles with strain, wavelength and amplitude can be efficiently tuned in some layered systems under compression, which exhibits important applications in stretchable electronics, gels, tissues and other fields [6–8]. Much interest has been focused on understanding the intrinsic mechanism of surface wrinkling phenomenon of film-substrate layered systems with extensive theoretical and numerical efforts. For the surface wrinkling instability of cylindrical structures, there are many previous works about cylindrical shells with or without homogenous cores. In stability, the buckling analysis of thin cylindrical shells has been concerned about since 1940s, as a representive study Koiter [9] made a remarkable contribution to buckling theory and a classical solution of a cylindrical shell under axial compression is established. By using Donnell’s equations, Seide [10] considered theoretically the axially and laterally compressed cylindrical shell supported by a soft core. Yao [11] used a theoretical solution to determine the buckling characteristics of a long cylinder shell with a solid or elastic soft core under axial compression with the sinusoidal buckling assumption. Herrmann and Forrestal [12] investigated the stability of a long thin cylindrical shell bonded to an elastic core, which is under a uniform external hydrostatic pressure. The linearized theory of elasticity is used to derive an expression for the buckling pressure. Later, Weingarten and Wang [13] utilized the finite element method to explore the buckling of a composite core-shell structures subjected to axisymmetric loads. Karam and Gibson [14] considered the thin cylindrical shell structure with a cellular core and gave out the potential of biomimicking of natural cylindrical shell structures on improving the mechanical efficiency of engineering structures. Hutchinson and He [15] studied the buckling of cylindrical sandwich shells with metal foam cores for the optimal design for structural weight and carrying capacity. By means of the energy method, Arani et al. [16] analyzed the elastic axisymmetric buckling of a thin simply supported and axially compressed cylindrical shell with an elastic core. In order to further investigate the effect of partially filled foam core on the buckling behavior, Lu et al. [17] proposed a theoretical formula to predict the critical buckling stress in a thin cylindrical shell with various percentage of infill. Liu et al. [18] studied the compressed buckling behaviors of a single-walled cylindrical shell filled with fluidlike materials via the energy-based analytical and numerical methods. Cao et al. [19] invoked the nonlinear Donnell-Mushtari-Vlassov shell equations to consider theoretically a core-shell structure on surface wrinkling along the axial and circumferential directions. The same theoretical model and the nonlinear shell elements were combined by Xu and Potier-Ferry [20] to discuss the wrinkling modes from a quantitative point of view. Besides, based on Euler-Bernoulli beam theory and Carrera’s Unified formulation, Hu et al. [21] presented the one-dimensional finite elements for the analysis of wrinkling in a stiff thin film resting on a thick elastic substrate. Then, stemming from the multi-scale method [22], they utilized the Fourier series with slowly variable coefficients to carry out a series of double scale analyses for effectively predicting wrinkling phenomena in film-substrate systems [23, 24], membrane

2

[25] and sandwich structures [26], which show promising potential in periodic instability problems. These analyses for layered structures often treat the elastic core substrate as a homogeneous one. However, in practice, biological organs such as plant stems, human blood vessels, earthworm bodies, etc. are common inhomogeneous substrates. Besides, the deposition of a stiff film or the attachment of a transient system device on a cylindrical soft core may cause the variation of material properties of the substrate in the radial direction. These materials of substrates are regarded to be gradient. To date, there are few studies concerning the wrinkling instability of a film bonded to cylindrical nonhomogeneous substrate. Cao et al. [27] explored the surface wrinkling of an elastic cylinder with modulus arbitrarily varying along radial direction, using a semi-analytical finite element method based on Koiter’s elasticity stability theory and the work of Lee et al. [28]. Wu et al. [29] used the state space method to determine the critical condition of surface instability of a graded elastic cylinder, which is divided into a sufficient amount of homogeneous cylindrical sublayers. Similarly, there are also several theoretical studies about wrinkling phenomenon of elastic in-plane layers with various material property along the depth [28, 30–32]. As abovementioned, wrinkling instability of a cylindrical layered structure was conducted by the classical elasticity theory or the classical finite element method. The former is mainly based on model simplification or theoretical assumption, which usually seems uneasy for solving the inhomogeneous media. When it comes to complex geometries, the latter has disadvantages of the geometric discrete errors and the high consumption of resources to improve accuracy. Therefore, in the present work a novel numerical method, namely isogeometric analysis (IGA) [33, 34] is employed to consider the wrinkling instability analysis. As a prevailing standard in computer-aided design (CAD), the nonuniform rational B-splines (NURBS) are used as the basis functions for exactly describing the geometry and physical fields. The high smoothness and continuity of the NURBS can efficiently avoid the discretization error and reduce costs [35]. In virtue of these advantages, IGA has been implemented to carry out the explicit finite deformation analysis of thin shell/membrane structures undergoing arbitrary deformations, which presents exciting results [36]. Meanwhile, morphological instabilities of living systems were studied quantitatively by IGA enhanced with a concurrent eigenvalue analysis [37]. By combining the IGA and the Fourier series, an explicit semi-analytical isogeometric analysis (SIGA) approach is developed, via the Fourier integral transformation in time and axial direction, specially for periodic wrinkling instability problems of cylindrical layered structures with modulus-gradient compliant substrate. The paper is organized as follows: in Section 2, the basic equations are briefly introduced by the linearized incremental theory. In Section 3, based on Floquet’s principle [38], the governing formulation with respect to critical wrinkling of cylindrical layered structures is obtained and parameterized by the NURBS basis functions, corresponding to the SIGA method. Section 4 represents the critical conditions for wrinkling instability of different film-substrate systems. The applications of the SIGA

3

method with two typical continuously graded modulus functions are demonstrated and the changing effects of structure parameters and material properties on critical wrinkling strain and wavelength are demonstrated in detail. Section 5 concludes the present study with a short summary.

2. Basic equations in incremental form

",

/

",

v

z

/

/ /

/ / /

/ / /

/ / / /

/

thinfilm

substγate

(a)

(b)

(c) Figure 1: Schematic diagram of cylindrical thin film-substrate system with modulus graded substrate in cylindrical coordinate system: (a) the initial (reference) configuration; (b) the wrinkle configuration; (c) the cross section of filmsubstrate structure; (c) observation points along the radial direction

2.1. Measures of Lagrangian deformation In the winkling problem of cylindrical layered systems as shown in Fig. 1(a), the total Lagrangian incremental formulation is typically employed for describing the loading and motion of bodies [39, 40]. According to the incremental theory, the deformation of the structure can be regarded as two configurations, namely the initial C0 and the current Cw increments, which represent undeformed (stress-free) and deformed (wrinkled) states, respectively. V0 and Vw denote the spatial domains with the boundaries ∂V0 and ∂Vw . The layered cylindrical system is under the static case without the inertia load, and two configurations can be expressed as x ≡ X0

in V0

x = φ(X, t) = X + u(X, t)

in Vw

4

(1)

in which, X and x are the material and spatial position vectors, respectively, u is the displacement of a material point corresponding to the wrinkling mode and φ is a mapping function from the initial state to the current. In the material space any position vector can be described as X(r, θ, z) = r cos θex + r sin θey + zez

(2)

In the cylindrical system, the natural bases at a point X is expressed by gi = ∂X/∂Xi , that is gr = cos θex + sin θey , gθ = −r sin θex + r cos θey , gz = ez

(3)

where ei (i = x, y or z) are the basis vectors in cartesian space, and the covariant bases in curvilinear coordinates can be defined by the vectors in Eq. (3), that is

g1 = cos θex + sin θey , g2 = − sin θex + cos θey , g3 = ez

(4)

The related covariant metric tensor Gij = gi gj are written as 



1

0

0

  G= 0  0

r2

  0   1

0

(5)

Accordingly, the contravariant bases gi satisfying gi gi = δij are given as 1 g , g3 = g3 r 2

g1 = g1 , g2 =

(6)

From Eqs. (4) and (6), by the definition gi = βij g(j) , the transformation tensor is expressed by  1   β = 0  0

0 r−1 0

 0   0  1

(7)

The deformation gradient F in material space is defined by F=

∂x = I + ∇0 u ∂X

(8)

where I is the unit tensor, ∇0 (= ∂/∂X) menas the material gradient and the Green strain tensor in the material space is E=

1 T (F F − I) = EL + EN 2

5

(9)

with EL =

1 [(∇0 u)T + (∇0 u)], 2

EN =

1 [(∇0 u)T (∇0 u)] 2

in which, T denotes the transpose of the matrix, EL and EN denote the linear and nonlinear terms of the strain tensor, respectively. ∇0 u = ∇uF = Uij gi gj (∇ = ∂/∂x means the spatial gradient) is the displacement gradient [41, 42], 

∂r ur

∂r uθ

  U = ∂θ ur /r − uθ /r  ∂z ur

∂r uz

∂θ uθ /r + ur /r ∂z uθ



  ∂θ uz /r  ∂z uz

(10)

represents the components of the gradient displacements. 2.2. Linearized variational governing equation of stressed structure On account of the principle of minimum potential energy, the potential energy of the film-substrate layered cylinder can be obtained from the difference between the stored strain energy by integrating the strain energy density over the undeformed configuration and the work caused by external loads. It is assumed that the applied loads are conservative and can be transformed to the undeformed state. The displacement space u lives in the space given by u(X) ∈ U, U = {u|u ∈ C0 (X), u = u on Γu }

(11)

where the displacements in U is often also called kinematically admissible displacements or compatible displacements. Meanwhile, they satisfy the displacement boundary conditions Γu and the continuity condition C 0 . The variation of functional can be written as δΠ = δPint − δPext Z Z = S : δEdV − V0

(12) tδudS

∂V

in which, δPint and δPext are the variations of the internal and external virtual energies. S means the second Piola-Kirchhoff stress with respect to the Green strain, t represents the external load distributed on the surface of body, δE and δu are the variations of the strain tensor and displacement. According to the Taylor series expansion, the deformation gradient F can be linearized in the direction of small displacement ∆u as DF(φ)[u] =

d F (φ + ξu) dξ ξ=0

= ∇0 ∆u

6

(13)

where D denotes the linearized operator, the term extracted by D proposes the first order of Taylor’s series, and ξ is a small parameter. With the directional derivative, the internal virtual power is linearized in the incremental form as Z DδPint [u] =

(D(S : δE)[u])dV V0

(14)

Z =

(δE : C : DE[u] + S : DδE[u])dV V0

where C is the 4th-order elasticity tensor. The variation of Green strain tensor is given as 1 T δF F + FT δF 2 1 = [(∇0 δu)T F + FT ∇0 δu] 2

δE =

(15)

= sym[(∇0 δu)T F] Similar to Eqs. (14) and (15), the linearized increment of strain tensor is obtained as DE[u] = sym[(∇0 ∆u)F]

(16)

where sym[·] denotes the symmetric part of a tensor, ∆u is the increment of displacement vector. The increment of the variation of the strain tensor is expressed as DδE[u] =

1 [(∇0 δu)T ∇0 ∆u + (∇0 ∆u)T ∇0 δu] 2

(17)

= sym[(∇0 δu)T ∇0 ∆u] and Z DδPext =

∆tδudS

(18)

∂V0

in which, ∆t is the increment of the mechanical forces. Here, it is assumed that the mechanical forces are constant and conservative, so that the linearized increment of external virtual work vanishes. Moreover, from the principle of minimum potential energy the variation in Eq. (12) must vanish if the system is in equilibrium. Therefore, we can get DδΠ = DδPint = 0

(19)

By combining Eqs. (14), (15), (18) and (19), the linearization of variation of the energy in Eq. (19) is obtained as Z

{sym[(∇0 δu)T F] : C : sym[(∇0 ∆u)F] + S : sym[(∇0 δu)T ∇0 ∆u]}dV = 0

V0

7

(20)

The above formulation is the governing equation for wrinkling instability of cylindrical layered structures, but it just provides the implicit expression. The critical wrinkling wavenumber and wavelength can not be obtained directly when the structure is about to wrinkle. Here is the explicit approach as the following.

3. SIGA formulation for wrinkling in cylindrical layered structures 3.1. Basic concepts of NURBS As the most prevalent approach to geometric representation, the NURBS are composed of linearly combination of piecewise polynomial B-spline basis functions and control points that can be visualized as points in coordinate space. Unlike the Lagrangian polynomial functions describing the geometry approximately, the NURBS basis functions furnish valuable properties of highly smoothness and continuity, which efficiently eliminate the Gibbs phenomena and provide a considerable reduction of computational time. The basic concept of NURBS for the exact representation of geometric model is briefly summarized and more details about the isogeometric can be found in the works [43]. In isogeometric analysis, B-spline functions are used to represent the geometric structures and parameterize the physical fields. A B-spline curve can be created by a knot vector defined by a nondecreasing real sequence Ξ = [ξ1 , ξ2 , · · · , ξn+p+1 ] (n and p denote the number of basis functions and the order) in the parametric space with 0 ≤ ξ ≤ 1. The knot vector is open and uniform. The first and last knots are repeated p + 1 times to satisfy the interpolating properties at boundaries. The knots partition the parameter space with elements. Given a knot vector, the B-spline basis functions can be constructed via the following recursive formula [44].   1, if ξ 6 ξ < ξ i i=1 Ni,0 (ξ) =  0, else Ni,p (ξ) =

(21)

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

Using the NURBS basis functions, the physical field can be described as

a=

Cpts X

Rpi (r(ξ))ai

(22)

i=1

with the one-dimensional basis functions Nip (ξ)ωi Rpi (ξ) = P n Njp (ξ)ωj j=1

where Cpts is the number of control points in isogeometric parametric mesh and ωi is the weight corresponding to the ith basis function, ξ represents the knot, i.e. the parametric coordinate and r is 8

the physical cartesian coordinate, a and ai denote the interpolated physical vector and the physical vector of the ith control point. 3.2. Explicit semi-analytical formulation for wrinkling instability For the axially compressed cylindrical layered structure, the physical fields of the wrinkle patterns after deformation are similar to that of elastic wave modes in waveguides. The linear wrinkle mode along the axial direction is considered and the pattern is periodic, as depicted in Fig. 1(b). By the virtue of Floquet’s principle [36] and the total Lagrangian formulation, the wrinkle phenomenon can be regarded as the static situation. The displacements and strains in polar coordinate system (r, θ, z) are functions of three spatial coordinates depending on eImθ (the integer m is the order in the circumferential direction, I is the imaginary unit) and eIkz (k is the wavenumber of linearly longitudinal wrinkle mode) along the longitudinal direction. Since this work focuses on the critical conditions of wrinkling instability, only the circumferential order m = 0 is considered, which means that the linear wrinkle patterns perform axisymmetric buckling modes. Thus, the displacement and strain fields are θ-independency, the three-dimensional axisymmetric problem can be reduced to the two-dimensional problem in spatial polar coordinates (r, z). The displacement field is expressed as a Fourier series along the longitudinal direction, that is

˜ (r)eIkz u(r, z) = u

(23)

with u(r, z) = [ur , uz ]T ,

˜ (r) = [˜ u ur , u ˜ z ]T

˜ represent the displacement vector and the amplitude vector in the deformed where, u(r, z) and u configuration Cw , I is the imaginary unit, z is the coordinate. The above formulation represents the generality of the wrinkle pattern avoiding the hypothesis of distribution of displacements for the wrinkle (buckle) mode. The related variational and incremental forms of displacements are expressed as ˜ (r)eIkz δu(r, z) = δ u (24) ∆u(r, z) = ∆˜ u(r)eIkz Since the physical field is θ-independency, it is sufficient to parameterize the cross-section of the axisymmetric structure only in the radial direction. In virtue of the use of a Fourier series in the longitudinal direction, the geometry of the cylinder can be described by one-dimensional NURBS basis function. The displacement vector at a generic point in geometric volume and boundary after the discretization can be described as u(ξ, z) = R(ξ)U(z),

χ ∈ V0

un (ξ, z) = Rn (ξ)Un (z),

χ ∈ ∂V0

9

(25)

where R and R represent the matrices containing the NURBS basis functions, which are taken in parametric coordinates ξ and η on V0 and ∂V0 , U and Un are the displacement vector of control points. In a parametric element space, the radial coordinate ri is represented by r(ξ) = Rpi ri

(26)

in which Rpi means the matrix of the NURBS basis function with the order p and i denotes the control point number. Each element can be described by the mapping from the parent element to the parametric space to the physical space and the global mapping determinant J = J1 ∗ J2 . The knot quadrature points in the parent element ξ ∈ [−1, 1] and the knot vectors are employed to calculate the parametric coordinates, the Jacobian determinant J2 =| (ξi+1 − ξi )/2 | for the ith element. Hence, the parametric coordinates are used to determine the physical coordinates and the derivatives of basis functions. The transformation between two kinds of coordinates is given by 

   ∂ ∂ = J1 ∂ξ ∂r

(27)

with the Jacobian determinant J1 =|

∂r |=| Rpi,ξ ri | ∂ξ

From Eqs. (9) and (10), the linear strain-displacement part can be rearranged as ε = Lr ue + Lz ue,z = (L0 /r + L1 ∂r )ue + Lz ue,z

(28)

with Eqs. (25) and (28), the linear strain vector that realizes the explicit wavenumber of the linear wrinkle mode is expressed as ε = (Lr + IkLz )Ru

(29)

which indicates that the linear part of strain can be contracted from the longitudinal coordinate z to the wavenumber k of the linear wrinkle mode using the spatial Fourier transform, which makes the abovementioned two-dimensional model reduced into the one-dimensional problem, that is the development of the present Fourier-based model, e means the parametric element, ue represents the displacements of knots in parametric element, u,z is the derivative of the displacement with respect to

10

the coordinate z, Lr = L0 /r + L1 ∂r and the compatibility operators  L0 =   L1 =   Lz = 

0

1

0

0

0

0

0

0

0

0

0

0

1

0

0

0

0

0

0

0

0

1

0

0

0

0

0

1

0

0

0

0

0

0

0

0

T  T (30)

 T 

For clarity, in the absence of body force, Eq. (20) can be rewritten in the spatial description as Z Z 2π r

{δ : C : d + σ : sym[(∇δu)T ∇0 ∆u]}drdz = 0

(31)

z

in which, the second-order tensors δ and d represent the first variation and the linearization of the Green strain tensor, the influence of the axial compression is included in the second term, which determines the critical condition of wrinkling instability. Besides, from Eqs. (25) to (29), the threedimensional problem can be turned into the one-dimensional in polar coordinates. By using Eq. (29), the above equation is converted to the matrix form of the wrinkle governing equation as follows ˜ δU

T

h   i ˜ =0 K1 + Ik K2 − KT2 + k 2 K3 ∆U

(32)

with Ke1 = Ke2 = Ke3 =

Z ZV

e

ZV

e

Ve

  RT LTr CLr + GTr Σσ Gr Rrd ξ

(33a)

  RT LTr CLz + GTr Σσ Gz Rrdξ

(33b)

  RT LTz CLz + GTz Σσ Gz Rrdξ

(33c)

with the stress matrix  σ  0  Σσ =  0  0

0 σ0 0

0



  0  σ0



σ0rr

  σ0 = σ0θr  σ0zr

σ0rθ σ0θθ σ0zθ

σ0rz



  σ0θz   σ0zz

(34)

where Σσ is the stress matrix, σ0 means the stress state at a point, V e means the element domain in parametric space, the stress state caused by axial compression can be obtained from the input strain ε and the differentiation operator matrices

11

 Gr =   Gz = 

∂r

0

0

0

r−1

0

0

0

0

0

0

0

0

0

0

∂r

0

0

0

0

1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

1

T  (35)

T 

Now, the critical wrinkling problem for the film-substrate layered cylinder under axial compression is transferred to a second-order polynomial eigenvalue case with respect to wavenumber k and axial strain ε. By the classic technique of matrix transformation to recast Eq. (32) to the first-order generalized explicit eigenvalue equation as ˜k = 0 [Ak − kBk ]Q with

 Ak = 



0

K1

K1

K2 − KT2



 Bk = 

(36)

K1

0

0

−K3





˜k =  Q

˜k U



˜k kU



Considering the relation between wavelength and wavenumber of linear wrinkle mode λ = 2π/k

(37)

where λ is the wavelength of the wrinkling mode. Eq. (36) can be also rewritten as ˜λ = 0 [Aλ − λBλ ]Q with

 Aλ = 

0

4π 2 K3

4π 2 K3

2IπK2 − KT2





 Bλ = 

(38)

4π 2 K3

0

0

−K1





˜λ =  Q

˜λ U



˜λ λU



which indicates that Ai and Bi (i represents k or λ, corresponding to two solving ways) are symmetric matrices. At each value of variable ε, there exist 2n eigenvalues and eigenvectors(n is the order of matrix). Eqs. (36) and (38) present the explicit semi-analytical isogeometric analysis formulation for wrinkling instability, which directly provide wavenumber or wavelength by an eigenvalue analysis. When the value comes to a certain critical strain, the above eigenvalue equation has the purely real solution, namely the critical wavenumber kc or wavelength λc . The corresponding eigenvectors represent the wrinkle profiles of displacements along the radial direction.

12

4. Numerical examples and discussion 4.1. Verification and convergence First of all, in order to validate the efficiency and accuracy of this proposed method, an example from Cao et al. [27] concerning a cylindrical film-substrate composite structure under axial compression is considered in cylindrical coordinate system (O − rθz), as illustrated in Fig. 1(a). The film and substrate are both homogeneous. Poisson’s ratio of the materials are assumed to be constant. hf and hs denote the thickness of stiff thin film and substrate, R is the outer radius. It is defined that the material properties are invariant along the axial direction, the compressive force is approximately employed by the axial strain uniformly distributed in the cylindrical layered structure. In this example, hf = 1 mm and hs = 100 mm, Poisson’s ratio is set as 0.4. The elastic modulus of substrate equals to 1 MPa while that of thin film varies from 40 to 10000 MPa by altering the modulus ratio κ = Ef /Es of stiff film to substrate. Fig. 2(a) shows the effects of various modulus ratio on the critical wrinkling strain. The relevant critical wrinkling wavelength is obtained directly from Eq. (38) and its dimensionless form can be taken as λc = λ/hf , as depicted in Fig. 2(b). 6

P re s e n t m e th o d C a o 's m e t h o d S h e ll th e o ry b y K o ite r

3 0

C r i t i c a l w a v e l e n g t h λc / h

5

C ritic a l s tra in (% )

4

3

2

2 0

P re s e n t m e th o d C a o 's m e t h o d S h e ll th e o ry b y K o ite r

1

1 0 0 1 0 0

1 0 0 0

E f/ E

1 0 0

1 0 0 0 0

1 0 0 0

E f/ E s

(a)

1 0 0 0 0

s

(b)

Figure 2: The variation of the critical condition for onset of wrinkling with the modulus ratio of film to substrate increasing: (a) the critical wrinkling strain εc vs the modulus ratio κ; (b) the critical wrinkling wavelength λc vs the modulus ratio κ.

It is obvious that the results of the present method for the critical condition of surface wrinkling of cylindrical film-substrate are consistent with that of Cao’s method. When the ratio κ increases, namely the film becomes gradually stiffer than the substrate, the bi-layer cylindrical structure degrades approximately to a single-layer homogeneous hollow cylindrical shell. The critical wrinkling condition converges to the classic solution of the thin hollow cylindrical shell given by Koiter [9]. Furthermore, for the sake of showing the advantage of the isogeometric analysis(IGA) over the finite element analysis(FEA), a comparative study of IGA versus FEA on solving the onset of surface wrinkling phenomenon in cylindrical film-substrate system is carried out using the proposed semianalytical approach based on IGA and FEM, respectively. The deduction process of the FEM-based

13

E f/E s= 1 0

1 e le m e n t 2 e le m e n ts

1 0

-2

1 0

-4

C ritic a l s tra in (re la tiv e e rro r % )

C ritic a l s tra in (re la tiv e e rro r % )

1 0

0

0

1 0

IG A F E A

4 e le m e n ts 5 e le m e n ts 6 e le m e n ts 8 e le m e n ts

1 0

-6

1 0

-8

1 0

2

1 0 e le m e n ts

-1 0

2 0 0

3 0 0

4 0 0

1 e le m e n t 2 e le m e n ts

1 0

-2

1 0

-4

1 0

-6

1 0

-8

1 0

1 0 0

E f/E s= 1 0 0 0

1 0

IG A F E A

4 e le m e n ts 5 e le m e n ts 6 e le m e n ts 8 e le m e n ts

1 0 e le m e n ts

-1 0

1 0 0

2 0 0

3 0 0

4 0 0

D o fs

D o fs

(a)

(b)

Figure 3: The relative error of the critical condition for onset of wrinkling instability when Ef /Es = 10 and Ef /Es = 1000, demonstrating comparison of accuracy with h-refinement between IGA and FEA: (a) the relative error of the critical strain; (b) the relative error of the critical wavelength.

approach is same as that of the IGA-based approach. The film-substrate system is regarded to be constructed by two patches(film and substrate) with a knot vector Ξ = {0, 0, 0, 1, 1, 1} and the order p = 2 of NURBS. An h-refinement analysis is performed by varying the mesh size parameter h (that is the length of the element along the radial direction) between hmin /hf = 0.1 and hmax /hf = 1 to evaluate the convergence of both theories. Here, the IGA wrinkling instability solution with respect to the finest mesh discretization (hmin /hf = 0.1) is chosen as the reference value for the analysis of verification and convergence of IGA- and FEA-based semi-analytical methods. Fig. 3 illustrates the relative error(percent error) of the critical condition for the onset of wrinkling of the thin-substrate layered cylinder. Two cases of the modulus ratio Ef /Es = 10 and Ef /Es = 1000 are considered, respectively. The convergence characteristics of semi-analytical formats based on both theories are shown from the coarsest to the finest mesh discretization. It is clearly noted that IGA furnishes more accurate results compared to FEA with a smaller number of degrees of freedom. The variation trend of relative error of critical wrinkling strain reveals the importance of the level of refinement across the film. Meanwhile, the overall trends for both cases are similar to each other, which states the convergences are not sensitive to variation of the modulus ratio. From the above observations in the analyses of verification and convergence, it is proved that the proposed IGA-based approach is reasonable and accurate. Subsequently, the semi-analytical isogeometric analysis approach is applied to wrinkling instability analysis of cylindrical composite structures with two typical graded materials, varying continuously along the radial direction of the substrate. The substrate modulus follows the exponential decaying function and the power-law distribution, as shown in Fig. 4, respectively. These two gradient functions

14

are expressed as Es = (Et − Ec )e(−α1 (R−r)) + Ec r Es = (Et − Ec )( )α2 + Ec R

(39)

where Et (= Ef /κi , i = 1, 2) and Ec (= βj Et , j = 1, 2) are the Young’s moduli at the surface and the center of the substrate, respectively, κ and β mean the modulus ratio of the film to the substrate surface and the modulus ratio of the surface to the center of cylindrical substrate, α1 and α2 denote the decaying rate and the power-law index, R is the radius of the cylindrical substrate and r represents the radial coordinate. The changes of substrate modulus along the radial direction, caused by the two gradient functions with various αi (i = 1, 2), are shown in Fig. 4. = 1

 1

= 1 0 0

 1

= 1 0 0 0

S u b s tra te m o d u lu s E

S u b s tra te m o d u lu s E

8 0

 2

= 1

 2

= 5

 2

= 1 0

s

8 0

1

(M P a )

1 0 0 

s

(M P a )

1 0 0

6 0

4 0

2 0

6 0

4 0

2 0

0

0 0

5

R a d ia l c o o rd in a te r /h

1 0

0

2 0

4 0

6 0

R a d ia l c o o rd in a te r /h f

(a)

8 0

1 0 0

f

(b)

Figure 4: Changing schematics of substrate modulus along the radial direction: (a) modulus variation following the exponential decaying dunction; (b) modulus variation following the power-law distribution.

4.2. The effects of significant parameters on wrinkling instability For the nonhomogeneous media, the semi-analytical IGA-based method is an effective approach to solve such problems. Firstly, an film-substrate system with the exponential decaying substrate is considered. According to the previous studies, the critical wrinkling strain εc and the non-dimensional wavelength λc /hf depend on the moduli of thin film and substrate, the thicknesses of both two, Poisson’s ratio and other coefficients. The influences of these parameters on wrinkling phenomenon should be systematically studied here. As shown in Fig. 5, the variation of thickness ratio hs /hf between substrate and film is investigated for the onset of wrinkling of the cylindrical film-substrate structure. The decaying ratio α1 is set as 1000, the modulus ratios β1 and κ1 are taken as 0.01 and 0.001. In Fig. 5(a) it is noted that the critical wrinkling strain decreases as the thickness ratio hs /hf of film to substrate increases while the wavelength is bigger, as shown in Fig. 5(b). For a small ratio hs /hf below 50, the critical wrinkling strain or the critical wrinkling wavelength is not sensitive to the variation of decaying ratio of substrate modulus. While for a large hs /hf , raising the decaying ratio is beneficial to the decrease of the wavelength. It demonstrates that the faster decaying ratio results in the smaller critical wrinkling strain and the longer dimensionless wrinkling wavelength. 15

= 1 1

= 1 0 0

 1

= 1 0 0 0

/h

1



C ritic a l W a v e le n g th

C ritic a l S tra in (% )

4





5

f

6

3

2

2

1 0

1

1 0 1 0

1

1 0

2

1 0

h s/h

3

1 0

4

1

1

1 0

1 0

2

1 0

h s/h f

(a)

3

 1

= 1

 1

= 1 0 0

 1

= 1 0 0 0 1 0

4

f

(b)

Figure 5: The effects of thickness variation on the critical condition for wrinkling of cylindrical film-substrate structure: (a) the critical wrinkling strain εc ; (b) the dimensionless wavelength λc . 3 6

C ritic a l W a v e le n g th

C ritic a l S tra in (% )



/h

f

0 .9

0 .8

0 .7

0 .6 1 0

0

1 0

1

1 0

2

D e c a y in g R a te

1 0 

3

1 0

3 4

3 2

3 0 4

1 0

0

1 0

1

1 0

2

D e c a y in g R a te 1

(a)

1 0 

3

1 0

4

1

(b)

Figure 6: The critical wrinkling condition of cylindrical film-substrate for various decaying ratio: (a) the critical wrinkling strain εc ; (b) the critical wrinkling wavelength λc .

Besides, Fig. 6 shows the critical wrinkling strain c and the dimensionless critical wrinkling wavelength λc /hf versus the decaying ratio α1 of graded modulus for a certain thickness ratio hs /hf = 100. For α1 < 50, the critical strain remains nearly unchanged at about 0.85% and the critical wrinkling wavelength is near 31. With the further increase of the decaying ratio, the strain gradually decreases while the wavelength increases. However, the overall changes of the critical wrinkling strain and the critical wrinkling wavelength, affected by the decaying ratio, are in small range. This situation at hs /hf = 100 can be also observed in Fig. 5. The effect of the modulus ratio κ of thin film to the surface of substrate on the critical condition for the onset of wrinkling is illustrated in Fig. 7, where the thickness ratio hs /hf = 100 and the decaying ratio α1 = 1000. The modulus ratio β1 of surface to center of cylindrical substrate is taken as 0.001, 0.01 and 0.1, respectively. The stiffness of the cylindrical substrate is dependent on the modulus ratio κ1 , which dominates wrinkling instability of cylindrical layered structures. From Fig.7(a), for the smaller modulus ratio κ1 the substrate modulus is comparable to that of the thin film, which represents 16

0 .2

= 0 .0 0 1 1

= 0 .0 1 1

= 0 .1

= 0 .0 0 1 1



1



1

/h 





= 0 .0 1 = 0 .1

C ritic a l W a v e le n g th

C ritic a l S tra in (% )



2

1 0

f

1



0 .1

1 0

1

0 .0 1 0

0

1 0

1

1 0

2

1 0

M o d u lu s ra tio



3

1 0

4

0

1 0

1

1 0

2

1 0

1

(a)

3

1 0

M o d u lu s ra tio



4

1 0

1

(b)

Figure 7: The critical wrinkling condition versus the modulus ratio κ1 and the influences of various decaying ratio β1 : (a) the critical wrinkling strain c ; (b) the critical wrinkling wavelength λc . 3 4

= 1 1

= 1 0 0 1

= 1 0 0 0

 

3 3



/h

1



0 .9 0

C ritic a l W a v e le n g th

C ritic a l S tra in (% )

f

0 .9 5

0 .8 5

3 2

 1

= 1

 1

= 1 0 0

 1

= 1 0 0 0

3 1

0 .8 0 3 0

1 0

-4

1 0

-3

1 0

-2

M o d u lu s ra tio

1 0 

-1

1 0

1 0 0

-4

1 0

-3

1 0

-2

M o d u lu s ra tio

1

(a)

1 0 

-1

1 0

0

1

(b)

Figure 8: The critical condition for the onset of wrinkling instability versus the modulus ratio β1 and the influences of decaying ratio α1 : (a) the critical wrinkling strain εc ; (b) the critical wrinkling wavelength λc .

the high resistance capacity to wrinkle. When the modulus ratio κ1 of thin film to substrate surface is below 10, the critical wrinkling strain c can be influenced by the modulus ratio β1 ratio of the top surface to the center of cylindrical substrate. Actually, the change of modulus ratio β1 leads to the variation of the elastic graded modulus of substrate. The critical strain increases as the modulus ratio raises while the wavelength reduces below the value 10. The curves for the critical wavelength λc /hf versus modulus ratio κ1 are not smooth and vary rapidly near the value 10. This abnormal phenomenon may be because the occurrence of mode transformation between global wrinkling mode and local surface wrinkling mode. The former means the wrinkling instability of the whole bilayer system with a smaller modulus ratio. For this case, the substrate modulus is comparable to that of stiff film, which makes two layers deform synchronously similar to a homogeneous solid cylinder. However, as the modulus ratio increases, the substrate becomes softer and softer. The surface stiff thin film firstly wrinkles before the deformation of the substrate reaches the threshold. This dramatic pattern transformation may holds promise for application in surface wrinkle pattern switching. Meanwhile,

17

the corresponding wrinkling wavelength decreases along with the modulus ratio β1 increasing. The minimal value of the curves gradually shifts to the negative direction. After the extreme, the curves of the critical wrinkling strain or wavelength are fitted into one curve, which has obviously positive correlation with the modulus ratio κ1 . To the end, the curves converges to one value, which is similar to the case of a stiff hollow cylinder. Moreover, when the thickness ratio hs /hf = 100, the modulus ratio κ1 = 1000 of film to substrate surface are invariant, the effects of various modulus ratio between the surface to the center on the onset of wrinkling instability for dfifferent decaying ratio α1 are investigated as shown in Fig. 8. The modulus ratio β1 determines the range of material graded modulus of substrate along the radial direction. The decaying ratio α1 indicates the distribution of material in substrate. When α1 = 1, the critical wrinkling strain remains almost constant but the critical wavelength decreases lightly. For the larger decaying ratio, the inhomogeneous is more evident, which can reduces the resistance capacity to wrinkling of composite structure, that is, the critical strain becomes smaller and the wavelength longer. Besides, it is observed that the curves for the critical wrinkling strain or wrinkling wavelength gather into a point. This is because the modulus ratio β1 = 1 means that the substrate is taken as homogeneous with the modulus Es = Et , which is independent of the decaying ratio α1 . 0 .3

= 1



/h

2

= 5

 2

= 1 0

C ritic a l W a v e le n g th

1 0 0 2

 2

= 1

 2

= 5

 2

= 1 0

f

0 .2



C ritic a l S tra in (% )



0 .1

1 0

0 .0 1 0

0

1 0

1

1 0

2

M o d u lu s ra tio

1 0 

3

1 0

4

1 0

0

1 0

1

1 0

2

M o d u lu s ra tio

2

(a)

1 0 

3

1 0

4

2

(b)

Figure 9: The critical condition for the onset of wrinkling instability versus the modulus ratio κ2 and the influences of power-law index α2 : (a) the critical wrinkling strain εc ; (b) the critical wrinkling wavelength λc .

Furthermore, the power-law distribution of substrate modulus is also considered. The modulus ratio between the stiff film and the substrate surface κ2 = 1000, the thickness ratio hs /hf = 100, the power-law index α2 = 1, 5 and 10 depicting the modulus variation in Fig. 4(b). Here, the effects of the modulus ratio κ2 and the power-law index α2 are studied and shown in Fig. 9. For the smaller value of modulus ratio κ2 , the critical wrinkling strain decreases as α2 increases. The material of substrate becomes softer and softer and it is more beneficial to wrinkling of composite structures. When α2 = 1, the material of substrate is a linear elastic graded material. From the figure, the power-law index can modulate the critical wrinkling strain and the critical wrinkling wavelength. Similarly, it is clearly 18

observed that there is an extremum on the curve of the critical wavelength versus the modulus ratio and the extremum shifts to the positive direction of the abscissa axis. This phenomenon indicates the mode transformation at a certain value similar to that in Fig. 7(b). With the increasing of the power-law index and the modulus ratio, the critical wrinkling strain has no apparent variations while the critical wavelength gets longer lightly. To the end, the curves for the critical wrinkling strain or critical wavelength get together to one curve. For the lager power-law index α2 , the film-substrate structure will gradually evolve to the shell without core. The critical wrinkling strain c and the p critical wrinkling wavelength λc will both converge into the solutions εc = hf /( 3(1 − v 2 )R0 ) and q p λc = πhf 2R0 / 3(1 − v 2 )hf obtained by Koiter [9]. r* = r* = r* = r* = r* =

0 .1 5 0 .1 5

0 .1 0

2 0 8 0 9 0 1 0 0 1 0 1

0 .0 5

0 .0 0

U

r

U

r

0 .1 0

-0 .0 5

0 .0 5

-0 .1 0

-0 .1 5

0 .0 0 0

2 0

4 0

6 0

R a d ia l c o o rd in a te r /h

8 0

-3 0

1 0 0

-2 0

-1 0

0

1 0

A x ia l c o o rd in a te z /h f

(a)

2 0

3 0

f

(b)

0 .1 r* = r* = r* = r* = r* =

0 .1 5

0 .1 0

2 0 8 0 9 0 1 0 0 1 0 1

0 .0

U

U

z

z

0 .0 5

-0 .1

0 .0 0

-0 .0 5

-0 .1 0 -0 .2

-0 .1 5 0

2 0

4 0

6 0

R a d ia l c o o rd in a te r /h

8 0

1 0 0

-3 0

-2 0

-1 0

0

A x ia l c o o rd in a te z /h f

(c)

1 0

2 0

3 0

f

(d)

Figure 10: Schematics of cylindrical film-substrate substrate following the exponential decaying function: (a) the radial displacement eigenvectors Ur ; (b) wrinkle profiles of displacement Ur along the axial direction; (c) the axial displacement eigenvectors Uz ; (d) wrinkle profiles of displacement Uz along the axial direction.

4.3. The displacement profiles of wrinkling instability In order to clearly represent the profiles of wrinkling, the displacement profiles for wrinkling instability of cylindrical film-substrate layered structure is shown in Fig. 10 and 11, according to Eqs. (22)

19

r* = r* = r* = r* = r* =

0 .1 5

0 .1 5

0 .1 0

2 0 8 0 9 0 1 0 0 1 0 1

0 .0 5

0 .0 0

U

r

U

r

0 .1 0

-0 .0 5

0 .0 5

-0 .1 0

-0 .1 5

0 .0 0 0

2 0

4 0

6 0

R a d ia l c o o rd in a te r /h

8 0

-0 .0 0 3

1 0 0

-0 .0 0 2

-0 .0 0 1

0 .0 0 0

0 .0 0 1

A x ia l c o o rd in a te z /h f

(a)

0 .0 0 2

0 .0 0 3

2 0

3 0

f

(b) 0 .0 8 r* = r* = r* = r* = r* =

0 .0 0

0 .0 6

2 0 8 0 9 0 1 0 0 1 0 1

0 .0 4 -0 .0 2

U

U

z

z

0 .0 2

-0 .0 4

0 .0 0

-0 .0 2 -0 .0 4

-0 .0 6

-0 .0 6 -0 .0 8

-0 .0 8 0

2 0

4 0

6 0

R a d ia l c o o rd in a te r /h

8 0

-3 0

1 0 0

-2 0

-1 0

0

A x ia l c o o rd in a te z /h f

(c)

1 0 f

(d)

Figure 11: Schematics of cylindrical film-substrate substrate following the power-law distribution: (a) the radial displacement eigenvectors Ur ; (b) wrinkle profiles of displacement Ur along the axial direction; (c) the axial displacement eigenvectors Uz ; (d) wrinkle profiles of displacement Uz along the axial direction.

and (35), which are for the exponential decaying function and the power-law distribution. The components of displacements are described by the radial normalized eigenvector Ur and axial normalized eigenvector Uz . From Fig. 1(c), these observation points (r∗ = 20, 80, 90, 100 and 101) are chosen to detect the wrinkling profiles of perturbed displacement components at different depths. It is obviously observed that the displacements mainly concentrate in and around the film while the distribution is not almost apparent away from the surface, as illustrated in Fig. 10 and 11.

5. Conclusions In this paper, we have presented an explicit semi-analytical isogeometric analysis method for the onset of wrinkling instability of thin films rested on cylindrical substrates with radially graded modulus. The wrinkling instability of film-substrate system can be described as an incremental problem within the total Lagrangian framework. The governing equation of stressed cylindrical layered structure is derived and parameterized by the NURBS. Then, the efficiency and accuracy of this method are

20

verified by comparing to the results in the reference. It is demonstrated that the isogeometric analysis is superior to the finite element analysis on solving the critical condition of wrinkling. Besides, the exponential decaying function and power-law distribution are employed to describe elastic graded material of substrates and wrinkling of thin films is predicted. The influences of various geometric and material parameters are considered on the critical wrinkling conditions. The displacement profiles for the onset of wrinkling are also shown.

Acknowledgements The authors wish to acknowledge the support from the Natural Science Foundation of China (Grant No. 11772130, 11972160 and 11902117) and the Fundamental Research Funds for the Central Universities, SCUT(Grant No. D2192630).

References [1] B. Li, Y. P. Cao, X. Q. Feng, H. J. Gao, Mechanics of morphological instabilities and surface wrinkling in soft materials: a review, Soft Matter 8 (2012) 5728–5745. [2] N. Bowden, S. Brittain, A. G. Evans, J. W. Hutchinson, G. M. Whitesides, Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer, Nature 393 (6681) (1998) 146. [3] M. B. Amar, A. Goriely, Growth and instability in elastic tissues, J Mech Phys Solids 53 (10) (2005) 2284–2319. [4] H. G. Allen, Analysis and design of structural sandwich panels: the commonwealth and international library: structures and solid body mechanics division, Elsevier, 2013. [5] H. Q. Jiang, D. Y. Khang, H. Y. Fei, H. Kim, Y. G. Huang, J. L. Xiao, J. A. Rogers, Finite width effect of thin-films buckling on compliant substrate: experimental and theoretical studies, J Mech Phys Solids 56 (8) (2008) 2585–2598. [6] J. A. Rogers, T. Someya, Y. G. Huang, Materials and mechanics for stretchable electronics, science 327 (5973) (2010) 1603–1607. [7] M. Guvendiren, S. Yang, J. A. Burdick, Swelling-induced surface patterns in hydrogels with gradient crosslinking density, Adv Funct Mater 19 (19) (2009) 3038–3045. [8] S. Budday, P. Steinmann, E. Kuhl, The role of mechanics during brain development, J Mech Phys Solids 72 (2014) 75–92. [9] W. T. Koiter, On the stability of elastic equilibrium. delft, Ph.D. thesis, Doctoral Thesis (1945). 21

[10] P. Seide, The stability under axial compression and lateral pressure of circular-cylindrical shells with a soft elastic core, J Aerosp Sci 29 (7) (1962) 851–862. [11] J. C. Yao, Buckling of axially compressed long cylindrical shell with elastic core, J App Mech 29 (2) (1962) 329–334. [12] M. J. FORRESTAL, G. HERRMANN, Buckling of a long cylindrical shell containing an elastic core, AIAA J 3 (9) (1965) 1710–1715. [13] V. I. Weingarten, Y. S. Wang, Stability of shells attached to an elastic core, J Eng Mech Div 102 (5) (1976) 839–849. [14] G. N. Karam, L. J. Gibson, Elastic buckling of cylindrical shells with elastic coresi. analysis, Int J Solids Struct 32 (8-9) (1995) 1259–1283. [15] J. W. Hutchinson, M. Y. He, Buckling of cylindrical sandwich shells with metal foam cores, Int J Solids Struct 37 (46-47) (2000) 6777–6794. [16] A. G. Arani, S. Golabi, A. Loghman, H. Daneshi, Investigating elastic stability of cylindrical shell with an elastic core under axial compression by energy method, J Mech Sci Technol 21 (7) (2007) 983–996. [17] L. Ye, G. Lu, L. S. Ong, Buckling of a thin-walled cylindrical shell with foam core under axial compression, Thin Wall Struct 49 (1) (2011) 106–111. [18] J. Wu, Q. H. Cheng, B. Liu, Y. W. Zhang, W. B. Lu, K. C. Hwang, Study on the axial compression buckling behaviors of concentric multi-walled cylindrical shells filled with soft materials, J Mech Phys Solids 60 (5) (2012) 803–826. [19] Y. Zhao, Y. P. Cao, X. Q. Feng, K. Ma, Axial compression-induced wrinkles on a core–shell soft cylinder: Theoretical analysis, simulations and experiments, J Mech Phys Solids 73 (2014) 212–227. [20] F. Xu, M. Potier-Ferry, On axisymmetric/diamond-like mode transitions in axially compressed core–shell cylinders, J Mech Phys Solids 94 (2016) 68–87. [21] J. Yang, Q. Huang, H. Hu, G. Giunta, S. Belouettar, M. Potier-Ferry, A new family of finite elements for wrinkling analysis of thin films on compliant substrates, Composite Structures 119 (2015) 568–577. [22] N. Damil, M. Potier-Ferry, H. Hu, Membrane wrinkling revisited from a multi-scale point of view, Advanced Modeling and Simulation in Engineering Sciences 1 (1) (2014) 6.

22

[23] Q. Huang, J. Yang, W. Huang, Y. Liu, H. Hu, G. Giunta, S. Belouettar, A new fourier-related double scale analysis for wrinkling analysis of thin films on compliant substrates, Composite Structures 160 (2017) 613–624. [24] Q. Huang, X. Rui, Y. Liu, H. Hu, G. Giunta, S. Belouettar, M. Potier-ferry, A two-dimensional fourier-series finite element for wrinkling analysis of thin films on compliant substrates, ThinWalled Structures 114 (2017) 144–153. [25] W. Huang, Q. Huang, Y. Liu, J. Yang, H. Hu, F. Trochu, P. Causse, A fourier based reduced model for wrinkling analysis of circular membranes, Computer Methods in Applied Mechanics and Engineering 345 (2019) 1114–1137. [26] Q. Huang, Y. Liu, H. Hu, S. Qian, K. Yu, G. Giunta, S. Belouettar, M. Potier-Ferry, A fourierrelated double scale analysis on the instability phenomena of sandwich plates, Computer Methods in Applied Mechanics and Engineering 318 (2017) 270–295. [27] F. Jia, Y. P. Cao, Y. Zhao, X. Q. Feng, Buckling and surface wrinkling of an elastic graded cylinder with elastic modulus arbitrarily varying along radial direction, Int J App Mech 6 (01) (2014) 1450003. [28] D. Lee, N. Triantafyllidis, J. R. Barber, M. D. Thouless, Surface instability of an elastic half space with material properties varying with depth, J Mech Phys Solids 56 (3) (2008) 858–868. [29] L. D. Han, C. X. Zhan, Y. H. Liu, Z. G. Wu, A state space solution for onset of surface instability of elastic cylinders with radially graded young’s modulus, Int J Solids Struct 126 (2017) 8–16. [30] Y. P. Cao, F. Jia, Y. Zhao, X. Q. Feng, S. W. Yu, Buckling and post-buckling of a stiff film resting on an elastic graded substrate, Int J Solids Struct 49 (13) (2012) 1656–1664. [31] Z. Chen, W. Q. Chen, J. ZSong, Buckling of a stiff thin film on an elastic graded compliant substrate, P Roy Soc A-Math Phy 473 (2208) (2017) 20170410. [32] J. J. Sui, J. B. Chen, X. X. Zhang, G. H. Nie, T. Zhang, Symplectic analysis of wrinkles in elastic layers with graded stiffnesses, J App Mech 86 (1) (2019) 011008. [33] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement, Comput Method App M 194 (39-41) (2005) 4135–4195. [34] C. L. Li, Q. Han, Y. J. Liu, Z. Wang, Wave isogeometric analysis method for calculating dispersive properties of guided waves in rotating damped cylinders, Meccanica 54 (1-2) (2019) 169–182. [35] G. Xu, B. Mourrain, R. Duvigneau, A. Galligo, Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications, Comput Aided Design 45 (2) (2013) 395–404. 23

[36] L. Chen, N. Nguyen-Thanh, H. Nguyen-Xuan, T. Rabczuk, S. P. A. Bordas, G. Limbert, Explicit finite deformation analysis of isogeometric membranes, Computer Methods in Applied Mechanics and Engineering 277 (2014) 104–130. [37] B. Dortdivanlioglu, A. Javili, C. Linder, Computational aspects of morphological instabilities using isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 316 (2017) 261–279. [38] M. Collet, M. Ouisse, M. Ruzzene, M. N. Ichchou, Floquet–bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems, Int J Solids Struct 48 (20) (2011) 2837–2848. [39] K. J. Bathe, Finite element procedures, Englewoodcliffs,NJ: Prentice-hall, 2006. [40] N. H. Kim, Introduction to nonlinear finite element analysis, Springer Science & Business Media, 2014. [41] W. Fl¨ ugge, Tensor analysis and continuum mechanics, Vol. 4, Springer, 1972. [42] Z. Huang, Fundamentals of continuum mechanics. [43] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA, John Wiley & Sons, 2009. [44] L. A. Piegl, W. Tiller, The Nurbs Book, Springer-Verlag, 1995.

24