An improved interpolating complex variable element free Galerkin method for the pattern transformation of hydrogel

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

An improved interpolating complex variable element free Galerkin method for the pattern transformation of hydrogel Yajie Deng a, Xiaoqiao He b,∗, Ligang Sun c, Shenghui Yi d, Ying Dai a,∗ a

School of Aerospace Engineering and Applied Mechanics, Tongji University, Shanghai 200092, China Department of Architecture and Civil Engineering, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong SAR, Hong Kong c College of Science, Harbin Institute of Technology, Shenzhen 518055, China d Centre for Advanced Structural Materials, City University of Hong Kong Shenzhen Research Institute, 8 Yuexing 1st Road, Shenzhen Hi-Tech Industrial Park, Nanshan District, Shenzhen, China b

a r t i c l e

i n f o

Keywords: Interpolating complex variable meshless method Hydrogel Pattern transformation Actuating behavior Area reduction

a b s t r a c t In this paper, an effective interpolating complex variable element free Galerkin method is proposed for the steady large deformation analysis of hydrogels. In this meshless method, the essential boundary conditions can be applied directly without using other special methods, which leads to less unknown coefficients. Through computing the swelling large deformation of a square hydrogel with three constrained boundaries, the presented meshless method is validated with higher accuracy and efficiency than other non-interpolating meshless methods. Based on the developed meshless method, different samples with elliptical, circular and square holes in a square hydrogel are simulated for analyzing the pattern transformation. The numerical results indicate that the shapes and initial sizes of holes in hydrogels lead to different pattern transformation rules, actuating behaviors and area reductions of holes. Finally, according to the new findings, the potential applications of a square hydrogel are discussed by designing different shapes and initial sizes of holes.

1. Introduction Hydrogel is a three-dimensional network of long cross-linked polymer chains and can absorb a mass of solvent molecules [1]. Because the three-dimensional network has the ability to support large and restorable deformation, the hydrogel can swell and shrink to different deformations by sensing and responding to the environmental stimulus, such as solvent, load, pH, light, temperature, electric field and so on [2–6]. Therefore, the hydrogel has attracted much attention from researchers and has the potential to be applied in engineering machines and biomedicine, such as the flow control actuators, soft machines, tissue engineering and drug delivery [7–10]. As one of the most common deformation behaviors, the mechanical instabilities can cause different and complicated pattern transformations in the hydrogel structures with periodic distributed holes [11–14]. Under the external stimulus, changing the form of the hydrogel structures will produce noticeable deformation which can be used in various useful aspects. Researchers had done a lot of work both in experimental and numerical simulation aspects. In experimental aspect, Singamaneni et al. described the details of the triggering and mechanistic aspects about the instabilities for pattern transformation in periodic porous elasto-plastic solids [11]. Zhu et al.

∗

fabricated the PHEMA based on hydrogel membranes in a square lattice with cylindrical holes and reported the capillary force induced instability, then observed the change of transparent of the membranes which can be used as a stimuli response window possibly [15]. Yang et al. designed a soft gripper and soft robots using the buckling deformation to achieve useful rotary motion, and also described the vacuum-actuated muscle-inspired pneumatic structures (VAMPs) to actuate a polymer duplicate as a human arm [16,17]. In numerical simulation aspect, by using the finite element analysis, Mullin et al. found the critical value of applied load for the uniaxial compression and novel transformations of the patterned structures [18]; Bertoldi et al. investigated the deformation behavior of two-dimensional elastomeric sheets with periodic patterns [19]; Okumura et al. studied the effect of geometrical imperfections and prestrains on swellinginduced buckling pattern in gel films with a square lattice of holes using the inhomogeneous field theory [20,21]. In addition, using the meshless method, Li et al. and Liu et al. simulated the swelling induced pattern transformation of lattice in gels [22,23]. In the investigations of hydrogels, most of researchers use the finite element method to simulate the pattern transformation. However, in the finite element simulation, a mesh-refinement scheme must be studied to guarantee the computing accuracy and convergence which leads

Corresponding authors. E-mail addresses: [email protected] (X. He), [email protected] (Y. Dai).

https://doi.org/10.1016/j.enganabound.2019.12.004 Received 14 April 2019; Received in revised form 20 November 2019; Accepted 10 December 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

to a low computing efficiency [18]. To avoid the mesh distortion and improve the computing efficiency and accuracy, the meshless method was developed based on a series of discrete nodes [24]. Compared with the finite element method, the meshless method does not produce mesh distortion during the computing process and has higher computing efficiency and accuracy [25–27]. As we all known, the most common used meshless method is the element free Galerkin (EFG) method built on the moving least squares (MLS) approximation [28]. Then, the complex variable moving least squares (CVMLS) and the improved complex variable moving least squares (ICVMLS) approximations were developed according to the MLS approximation with high computing efficiency [29,30]. On the foundation of the CVMLS approximations, the complex variable element free Galerkin (CVEFG) and the improved complex variable element free Galerkin (ICVEFG) methods were built for different problems, such as the elasticity and viscoelasticity problems [31–33]. However, the trial functions in the MLS and CVMLS approximations cannot satisfy the property of Kronecker 𝛿 function, and the corresponding EFG and CVEFG methods must rely on special techniques to dispose the essential boundary conditions with increased unknown coefficients [34,35]. When using the penalty function method to handle the essential boundary conditions, the choice of penalty factor 𝛼 is dependent on experience. Therefore, for some complicated and new problems, choose an appropriate penalty factor is time consuming. To overcome this disadvantage in the meshless method, the interpolating moving least squares (IMLS) method was proposed by Lancaster et al. [36], in which the shape function satisfies the property of Kronecker 𝛿 function [37]. To improve the computing efficiency, Ren and Cheng improved the expression form of the shape function and used the corresponding interpolating meshless method to solve the potential problem and elasticity problem [38,39]. Then Deng et al. applied the interpolating idea to the complex variable meshless method and presented the interpolating complex variable moving least squares method with high computing accuracy and efficiency [40]. In this paper, an improved interpolating complex variable moving least squares (IICVMLS) method is applied, in which a complete basis function is used to guarantee the computing accuracy [41]. Then combining the IICVMLS method with the fundamental deformation field theory of hydrogel, and considering the multiplicative decomposition of deformation gradient [22,42,43], an improved interpolating complex variable element free Galerkin (IICVEFG) method for the large deformation analysis of hydrogels is proposed in present work. The accuracy and efficiency of the IICVEFG method are validated through simulating a simply constrained swelling deformation of hydrogel. Then using the proposed IICVEFG method, mechanical behaviors of a square hydrogel with various shapes and initial sizes of holes are simulated systematically, and the effects of shape and size of holes in hydrogels are discussed in detail for the actuation and area control of hydrogels.

In this work, the following well-known Flory–Rehner free energy function is adopted as [46,47] 𝑊 (𝑭 , 𝐶 ) =

[ ( ) 𝜒 ] 1 kT 1 vC ln 1 + + , NkT(𝐼 − 3 − 2 ln 𝐽 ) − 2 𝑣 vC 1 + vC

where N is the number of polymeric chains per reference volume, kT is the absolute temperature in the unit of energy, v is the volume of per solvent molecule, 𝜒 is a dimensionless measure of the enthalpy of mixing, I = Fij Fij and 𝐽 = 𝑑𝑒𝑡(𝑭 ) are the invariants of deformation gradient F. In this material model, the incompressibility of molecules 1 + vC = J is considered, and through the Legendre transformation the equilibrium Eq. (1) can be written as [42] ∫Ω

𝛿 𝑊̂ dΩ =

∫Ω

𝑏̄ 𝑖 𝛿𝑥𝑖 dΩ +

∫Γ𝑡

𝑡̄𝑖 𝛿𝑥𝑖 dΓ,

𝑭 = 𝑭 𝜇 𝑭 ′,

(4)

where F𝜇 is the free swelling deformation gradient from state C0 to C𝜇 , and F′ is the mechanical deformation gradient from state C𝜇 to C. In the reference state C𝜇 , the corresponding free energy function 𝑊̂ ′ (𝑭 ′ , 𝜇) can be written as ( )) 𝑣 ̂ ′( ′ ) Nv ( 𝑊 𝑭 ,𝜇 = 𝐼 − 3 − 2 ln 𝐽𝜇 𝐽 ′ 2𝐽𝜇 kT ( ) ) 𝐽𝜇 𝐽 ′ 𝜒 𝜇 ( ′ − − − 𝐽 ′ − 𝐽𝜇−1 ln 𝐽 − 𝐽𝜇−1 , (5) 𝐽𝜇 𝐽 ′ − 1 𝐽𝜇2 𝐽 ′ kT where 𝐽 ′ = 𝑑𝑒𝑡(𝑭 ′ ) and 𝐽𝜇 = 𝑑𝑒𝑡(𝑭 𝜇 ). 2.2. Mechanical parameters According to the above free energy function, the nominal second Piola-–Kirchhoff stress of hydrogel is defined as ( ′ ) ( ′ ) ′ ′ 𝑣 ̂′ 𝑣 𝜕 𝑊̂ 𝑭 , 𝜇 𝑣 𝜕 𝑊̂ 𝑭 , 𝜇 𝑆ij = = 2 𝜕 𝐸 ′ ij 𝜕 𝐺′ ij kT kT kT ( ) 2 Nv𝜆𝜇 𝐽𝜇 𝐽 ′ 𝜒 𝜇 ′ ′ −1 1 − Nv = 𝛿ij + −𝐽 ′ ln + − 𝐽 𝐺 ji , 𝐽𝜇 𝐽𝜇 𝐽𝜇 𝐽 ′ − 1 𝐽𝜇2 𝐽 ′ kT

2.1. Free energy function

where E′ is the Green–Lagrange strain tensor

When a dry network of hydrogel, with geometric constraints and external load, is soaked in a solvent of chemical potential 𝜇, the hydrogel will swell and produce large deformation until the hydrogel reaches the equilibrium state. According to the equilibrium theory presented by Hong et al. [42,44,45], the equilibrium equation is

𝑬′ =

𝛿𝑊 (𝑭 , 𝐶 )dΩ =

∫Ω

𝑏𝑖 𝛿𝑥𝑖 dΩ +

∫Γ𝑡

𝑡𝑖 𝛿𝑥𝑖 dΓ + 𝜇

∫Ω

𝛿𝐶dΩ,

(3)

where 𝑊̂ (𝑭 , 𝜇) = 𝑊 (𝑭 , 𝐶) − 𝜇𝐶 is a new free energy function. When the hydrogel swells without mechanical load, the above Eq. (3) can be expressed as ∫Ω 𝛿 𝑊̂ dΩ = 0. To solve this equation, Li et al. built the swelling deformation model of hydrogel under geometric constraints by using the multiplicative decomposition method [22]. According to the multiplicative decomposition method, the total deformation gradient F is yielded during the deformation from the original state C0 to the current state C. Supposing a reference state C𝜇 with the chemical potential 𝜇 and free stress, the whole deformation gradient F can be divided into two parts

2. Equilibrium field theory of hydrogel

∫Ω

(2)

) 1( ′ ) 1 ( ′T ′ 𝑭 𝑭 −𝑰 = 𝑮 −𝑰 , 2 2

(6)

(7)

and I is the unit tensor, G′ is the right Cauchy–Green tensor and the expression form is G′ = F′T F′. In the free swelling state, the nominal stress 𝑣 ̂′ 𝑆 = 0 and Eq. (6) gives the relationship between the free swelling 𝑘𝑇 𝑖𝑗 stretch 𝜆′ and chemical potential 𝜇 [22]. And the nominal tangent constitutive matrix is

(1)

where Ω is the domain of hydrogel, Γ = Γu ∪ Γt is the boundary of Ω, Γu is the specified displacement boundary, Γt is the given force boundary, 𝑏̄ 𝑖 and 𝑡̄𝑖 are the forces applied on the volume and area, respectively. W(F, C) is the free energy function which relates to the deformation gradient F and solvent concentration C.

𝐷𝑖𝑗𝑘𝑙 =

𝑣 ̂′ 𝑆𝑖𝑗 𝜕 𝑘𝑇

𝜕𝐸 ′ (

=2

100

𝑘𝑙

=2

𝑣 ̂′ 𝑆𝑖𝑗 𝜕 𝑘𝑇

𝜕 𝐺′ 𝑘𝑙

) 𝐽𝜇 𝐽 ′ 𝜒 𝜇 ′ −1 1 − 𝑁𝑣 −1 − 𝐽 ′ ln + − 𝐽 𝑄𝑖𝑗𝑘𝑙 𝐽𝜇 𝐽 ′ − 1 𝐽 2 𝐽 ′ 3 𝑘𝑇 𝐽𝜇 𝐽 ′ 2 𝜇

Y. Deng, X. He and L. Sun et al.

( +

Engineering Analysis with Boundary Elements 113 (2020) 99–109

𝐽𝜇 𝐽 3𝜒 𝜇 ′ 2𝑁𝑣 − 2 𝐽′ + + 𝐽 ′ ln + 𝐽 − 𝐽𝜇 𝐽𝜇 𝐽 ′ − 1 𝐽𝜇 𝐽 ′ − 1 𝑘𝑇 𝐽𝜇2 𝐽 ′

× (𝐺′ 𝑗𝑖

′

−1

𝐺′ 𝑙𝑘

−1

),

)

The local function at point z can be expressed as 𝑢ℎ (𝑧, 𝑧̂ ) =

𝑖=1

(8)

𝑸 = 𝑄ijkl

)

⎡ 0 ⎢ = ⎢ 𝐺′ 33 ⎢ ⎣ 0

𝐺′ 33

0

0

0

0

−𝐺

′

33

⎤ ⎥ ⎥. ⎥ ⎦

(9)

2.3. Equilibrium equation of hydrogel

(10)

= 𝑝 𝑖 (𝑧 ) −

The displacement boundary conditions in Γu is Δu′ = (1 − 𝜆𝜇 ) × Xi , where 𝜆𝜇 is the free swelling stretch in the chemical potential 𝜇 and Xi is the coordinate in the dry network state. Substituting Eq. (8) into Eq. (10), it is obtained that T

∫Ω

𝛿𝜺′ ⋅ 𝑫 𝜺′ dΩ =

𝑣 𝑣 𝑏 𝛿𝑥′ dΩ + 𝑡 𝛿𝑥′ dΓ, ∫Ω kT 𝑖 𝑖 ∫Γ𝑡 kT 𝑖 𝑖

𝐼=1

𝒗(𝑧) = (𝑣(𝑧 − 𝑧1 ), 𝑣(𝑧 − 𝑧2 ), … , 𝑣(𝑧 − 𝑧𝑛 ))𝑇 .

(11)

(23) (24)

The normalized weight function is derived through letting the term T number m of the conjugate basis function 𝒑 (𝑧) equal 1 in the ICVMLS approximation [29]. Therefore in the case m = 1 the trial function is ∑𝑛 𝑢(𝑧𝐼 )𝑤(𝑧 − 𝑧𝐼 ) 𝑢ℎ (𝑧) = 𝑢𝑠 (𝑧) = 𝐼=1 = (𝑢, 𝛽𝑧(1) (𝑧))𝑧 𝛽𝑧(1) (𝑧), (25) ∑𝑛 𝐼=1 𝑤(𝑧 − 𝑧𝐼 ) which is a weighted average of the function values at nodes zI in the influence domain of point z. (𝑚) The 𝛽𝑧(1) (𝑧), 𝑏(2) 𝑧 (𝑧), … , 𝑏𝑧 (𝑧) constitute the new basis function of the IICVMLS method. Hence, the trial function in Eq. (17) is expressed as the following interpolating function

(12)

where Δ𝜺′ is the increment of strain and Δ𝜺′ can be divided into the linear term Δe and quadratic term Δ𝜼 which are related to the displacement increment Δu′i , Δ𝜺′ = Δ𝒆 + Δ𝜼,

𝑛 ∑ ( ) ( ) 𝑝𝑖 𝑧𝐼 𝑣 𝑧 − 𝑧𝐼 = 𝑝𝑖 (𝑧) − 𝒗T (𝑧)𝒑, (𝑖 = 2, 3, … , 𝑚), (22)

where v(z − zi ) is the normalized weight function, 𝑤(𝑧 − 𝑧𝑖 ) 𝑣 ( 𝑧 − 𝑧 𝑖 ) = ∑𝑛 , (𝑖 = 1, 2, … , 𝑛), 𝐼=1 𝑤(𝑧 − 𝑧𝐼 )

where 𝜺′ is the vector of strain. For the swelling large deformation of hydrogel, the increment form is introduced to solve the nonlinear Eq. (11). Suppose that the strain 𝜺′l is known, the unknown strain 𝜺′l + 1 is solved according to the following equation, 𝜺′𝑙+1 = 𝜺′𝑙 + Δ𝜺′ ,

(20)

where w(z − zI ) is a singular weight function with compact support, and zI are the nodes whose influence domains cover point z. Let pi (z) (i = 2, 3, ..., m) be orthogonal to 𝛽𝑧(1) (𝑧), it is obtained that ( ) 𝑏(𝑧𝑖) (𝑧) = 𝑝𝑖 (𝑧) − 𝑝𝑖 (𝑧), 𝛽𝑧(1) (𝑧) 𝑧 𝛽𝑧(1) (𝑧)

On the basis of the above equilibrium theory, in the reference state the new equilibrium equation is built as 𝑣 ̂ ′ 𝑣 ̄ ′ 𝑣 ̄ ′ 𝑏 𝛿𝑥 dΩ + 𝑡 𝛿𝑥 dΓ. 𝛿 𝑊 dΩ = ∫Ω 𝑘𝑇 ∫Ω 𝑘𝑇 𝑖 𝑖 ∫Γ𝑡 𝑘𝑇 𝑖 𝑖

𝑝𝑖 (𝑧̂ )𝑎𝑖 (𝑧) = 𝒑T (𝑧̂ )𝒂(𝑧),

where 𝑧̂ is the node whose influence domain covers the point z. In order to obtain an interpolating shape function, the inner product and corresponding norm are used [29], and the first term p1 (z) ≡ 1 is normalized on the space 𝑠𝑝𝑎𝑛 (𝑝1 , 𝑝2 , … , 𝑝𝑚 ) as [39] 𝑝1 1 𝛽𝑧(1) (𝑧) = = , (21) ‖𝑝 ‖ [∑𝑛 ]1 ‖ 1 ‖𝑧 𝑤 ( 𝑧 − 𝑧𝐼 ) 2 𝐼=1

where (

𝑚 ∑

𝑚 ∑ ( ) 𝑢ℎ (𝑧) = 𝑢, 𝛽𝑧(1) (𝑧) 𝑧 𝛽𝑧(1) (𝑧) + 𝑎𝑖 (𝑧)𝑏(𝑧𝑖) (𝑧) = 𝒗T (𝑧)𝒖 + 𝒃T (𝑧)𝒂(𝑧),

(13)

(26)

𝑖=2

Δ𝑒𝑖𝑗 =

1 (Δ𝑢′𝑖,𝑗 + Δ𝑢′𝑗,𝑖 + 𝑢′𝑘,𝑖 Δ𝑢′𝑘,𝑗 + 𝑢′𝑘,𝑗 Δ𝑢′𝑘,𝑖 ), 2

1 Δ𝜂𝑖𝑗 = Δ𝑢′𝑘,𝑖 Δ𝑢′𝑘,𝑗 . 2

where

(14) (15)

Through linearization, the weak integral form of Eq. (11) can be written as ∫Ω

𝑫 𝜺′ 𝑙 𝛿Δ𝜼dΩ +

∫Ω

𝑫 Δ𝒆𝛿Δ𝒆dΩ

𝑣 𝑣 𝑏 𝛿𝑥′ dΩ + 𝑡 𝛿𝑥′ dΓ − 𝑫 𝜺′ 𝑙 𝛿Δ𝒆dΩ. = ∫Ω kT 𝑖 𝑖 ∫Γ𝑡 kT 𝑖 𝑖 ∫Ω

( )T 𝒂(𝑧) = 𝑎2 (𝑧), 𝑎3 (𝑧), … , 𝑎𝑚 (𝑧) ,

(27)

( )T 𝒃(𝑧) = 𝑏(𝑧2) (𝑧), 𝑏(𝑧3) (𝑧), … , 𝑏(𝑧𝑚) (𝑧) ,

(28)

( ( ) ( ) ( ))T 𝒖 = 𝑢 𝑧1 , 𝑢 𝑧2 , … , 𝑢 𝑧𝑛 .

(29)

Then using the functional J presented by Bai et al. [29], we can obtain the unknown coefficient vector a(z) as following to ensure the functional J is minimum,

(16)

𝒂(𝑧) = 𝑨−1 𝑧 ( 𝑧 ) 𝑩 𝑧 ( 𝑧 ) 𝒖,

3. The IICVMLS method

(30)

where T

𝑨𝑧 (𝑧) = 𝑪 𝑧 𝑾 (𝑧)𝑪 𝑧 ,

In the IICVMLS method, the trial function can be written as 𝑢ℎ (𝑧) = 𝑢ℎ1 (𝑧) + 𝑖𝑢ℎ2 (𝑧) = (

𝑚 ∑ 𝑖=1

) 𝑧 = 𝑥1 + 𝑖𝑥2 ∈ Ω ,

T

𝑝𝑖 (𝑧)𝑎𝑖 (𝑧) = 𝒑T (𝑧)𝒂(𝑧),

𝑩 𝑧 (𝑧) = 𝑪 𝑧 𝑾 (𝑧), ( ) ⎡𝑏(𝑧2) 𝑧1 ⎢ (2 ) ( ) 𝑏 𝑧2 𝑪𝑧 = ⎢ 𝑧 ⎢ ⋮ ⎢ (2 ) ( ) ⎣𝑏 𝑧 𝑧 𝑛

(17)

where pT (z) = (p1 (z),p2 (z),..., pm (z)) is the complete basis function vector, and aT (z) = (a1 (z),a2 (z),..., am (z)) is the unknown coefficient vector. In the two-dimensional domain, the linear and quadratic complete basis function vectors are [41] ( ) 𝒑T (𝑧) = 1, 𝑧, 𝑧 , (18) ) 𝒑T (𝑧) = 1, 𝑧, 𝑧, 𝑧2 , 𝑧2 , 𝑧𝑧 .

(31)

(32) 3 ( ) 𝑏(𝑧 ) 𝑧1 ( ) 𝑏(𝑧3) 𝑧2 ⋮ 3 ( ) 𝑏(𝑧 ) 𝑧𝑛

( ) ⎡𝑤 𝑧 − 𝑧1 ⎢ 0 𝑾 (𝑧 ) = ⎢ ⎢ ⋮ ⎢ 0 ⎣

(

(19) 101

⋯ ⋯ ⋱ ⋯

0 ( ) 𝑤 𝑧 − 𝑧2 ⋮ 0

𝑚 ( ) 𝑏(𝑧 ) 𝑧1 ⎤ ( )⎥ 𝑏(𝑧𝑚) 𝑧2 ⎥ , ⋮ ⎥ (𝑚) ( )⎥ 𝑏𝑧 𝑧𝑛 ⎦

⋯ ⋯ ⋱ ⋯

⎤ 0 ⎥ 0 ⎥. ⎥ ⋮ ( )⎥ 𝑤 𝑧 − 𝑧𝑛 ⎦

(33)

(34)

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Substituting Eq. (30) into Eq. (26) yields 𝑢ℎ (𝑧) = 𝒗T (𝑧)𝒖 + 𝒃T (𝑧)𝑨−1 𝑧 (𝑧)𝑩 𝑧 (𝑧)𝒖 = 𝚽(𝑧)𝒖 =

𝑛 ∑ 𝐼=1

𝑲 (𝐿1) = ( ) Φ𝐼 (𝑧)𝑢 𝑧𝐼 ,

where 𝚽(z) is the interpolating shape function ( ) 𝚽(𝑧) = Φ1 (𝑧), Φ2 (𝑧), … , Φ𝑛 (𝑧) = 𝒗T (𝑧) + 𝒃T (𝑧)𝑨−1 𝑧 (𝑧)𝑩 𝑧 (𝑧),

T

∫Ω

𝑩 (𝐿0) 𝑫 𝑩 (𝐿1) dΩ+

T

∫Ω

𝑩 (𝐿1) 𝑫 𝑩 (𝐿0) dΩ+

T

∫Ω

𝑩 (𝐿1) 𝑫 𝑩 (𝐿1) dΩ, (48)

(35)

𝑲 (NL) =

(36)

∫Ω

𝑩 T 𝑺 𝑩 dΩ,

(49)

and the following singular weight function is used [40] )2 ⎧ 2 ( ⎪ 𝜌 1 − 𝑧−𝜌𝑧 | 𝐼| 𝑤(𝑧 − 𝑧𝐼 ) = ⎨ |𝑧−𝑧𝐼 |2 ⎪0 ⎩

|𝑧 − 𝑧 | ≤ 𝜌 𝐼| | , |𝑧 − 𝑧 𝐼 | ≥ 𝜌 | |

𝒇1 =

𝒇2 =

𝑛 ∑ 𝐼=1

( ) ∗ Φ𝐼 (𝑧)Δ𝑢′ 𝑧𝐼 = 𝚽(𝑧)Δ𝒖′ ,

(38)

[ ] ⎡Re[Φ𝐼,1 (𝑧)] ⎢Re Φ (𝑧) 𝑩 (𝑧) = ⎢ [ 𝐼,2 ] ⎢Im[Φ𝐼,1 (𝑧)] ⎣Im Φ𝐼,2 (𝑧) ⎡𝑆̂ ′ ⎢ 11 ′ 𝑣 ⎢𝑆̂12 𝑺= kT ⎢ 0 ⎢ 0 ⎣

𝐼=1

[ ] [ ] ] Re Φ𝐼 (𝑧) − Im Φ𝐼 (𝑧) , [ ] [ ] Im Φ𝐼 (𝑧) Re Φ𝐼 (𝑧)

̃ 𝑧) = (𝚽 ̃ 𝐼 (𝑧), 𝚽 ̃ 2 (𝑧 ), … , 𝚽 ̃ 𝑛 (𝑧)), 𝚽( [ ′ ] Δ𝑢 1 (𝑧𝐼 ) Δ𝒖′𝐼 = , Δ𝑢′ 2 (𝑧𝐼 )

(41)

(42)

(44)

(45)

𝑲 = 𝑲 (𝐿0) + 𝑲 (𝐿1) + 𝑲 (NL) ,

T

𝑩 (𝐿0) 𝑫 𝑩 (𝐿0) dΩ,

0 0 ′ 𝑆̂11 𝑆̂ ′ 12

0 ⎤ ⎥ 0 ⎥ ′ ⎥, ̂ 𝑆12 𝑆̂ ′ ⎥⎦

(54)

(55)

22

(56)

(57)

(58)

5. Numerical examples

where

∫Ω

⎤ ⎥ ⎥, ⎥ ⎦

Because the shape function of IICVMLS method in Eq. (36) possesses the interpolation property, the IICVEFG method can apply the essential boundary conditions directly without special techniques. The expressions of Eqs. (46) and (50) in the IICVEFG method are more concise than other non-interpolating complex variable meshless methods, which leads to higher numerical accuracy and efficiency for the large deformation of hydrogels.

(43)

Then according to Eq. (16), the final algebraic equations are obtained 𝑲 Δ𝒖′ = 𝒇 1 − 𝒇 2 ,

′ 𝑆̂12 ′ ̂ 𝑆22 0 0

[ ] −Im[Φ𝐼,1 (𝑧)]⎤ ⎥ −Im [ Φ𝐼,2 (𝑧]) , Re[Φ𝐼,1 (𝑧)] ⎥⎥ Re Φ𝐼,2 (𝑧) ⎦

𝑣 ̂′ ̂′ ̂′ T 𝑺̂ = (𝑆 , 𝑆 , 𝑆 ) , 𝑘𝑇 11 22 12 [ ] 𝑣 𝑏̄ 1 𝒃̄ = , 𝑘𝑇 𝑏̄ 2 [ ] 𝑣 𝑡̄1 𝒕̄ = . 𝑘𝑇 𝑡̄2

( ( ) ( ) ( ) ( ) ( ) ( ))T Δ𝒖′ = Δ𝑢′ 1 𝑧1 , Δ𝑢′ 2 𝑧1 , Δ𝑢′ 1 𝑧2 , Δ𝑢′ 2 𝑧2 , … , Δ𝑢′ 1 𝑧𝑛 , Δ𝑢′ 2 𝑧𝑛 .

𝑲 (𝐿 0 ) =

(51)

(53)

From Eq. (38), we can gain that ] 𝑛 [ [ ] [ ] ][ ( ) ] ∑ Re Φ𝐼 (𝑧) Δ𝑢′ 1 𝑧𝐼 −Im Φ𝐼 (𝑧) = [ ] [ ] ( ) Δ𝑢′ 2 (𝑧) Re Φ𝐼 (𝑧) Im Φ𝐼 (𝑧) Δ𝑢′ 2 𝑧𝐼 𝐼=1 (40) 𝑛 ∑ ∼ ̃ 𝑧)Δ𝒖′ , = 𝚽𝐼 (𝑧)Δ𝒖′ 𝐼 = 𝚽(

𝚽𝐼 (𝑧) =

T ̂ dΩ, 𝑩 (𝐿1) 𝑺

[ ] [ ] 𝐿1 1 Re [Φ𝐼, 1 (𝑧) ] + 𝐿2 1 𝐼𝑚[ Φ𝐼, 1 (𝑧)] ⎡ =⎢ [ ] 𝐿1 2 Re[ Φ𝐼, 2 (𝑧)] + 𝐿22 𝐼𝑚 [Φ𝐼, 2 (𝑧) ] [ ] ⎢ ⎣𝐿1 1 Re Φ𝐼, 2 (𝑧) + 𝐿1 2 Re Φ𝐼, 1 (𝑧) + 𝐿2 1 𝐼 𝑚 Φ𝐼, 2 (𝑧) + 𝐿22 𝐼 𝑚 Φ𝐼, 1 (𝑧) , [ ] [ ] −𝐿1 1 𝐼𝑚 [Φ𝐼, 1 (𝑧) ] + 𝐿2 1 Re[ Φ𝐼, 1 (𝑧)] ⎤ ⎥ [ ] −𝐿1 2 𝐼𝑚 [Φ𝐼, 2 (𝑧) ] + 𝐿22 Re [Φ𝐼, 2 (𝑧) ] [ ]⎥ −𝐿1 1 𝐼 𝑚 Φ𝐼, 2 (𝑧) − 𝐿1 2 𝐼 𝑚 Φ𝐼, 1 (𝑧) + 𝐿2 1 Re Φ𝐼, 2 (𝑧) + 𝐿22 Re Φ𝐼, 1 (𝑧) ⎦

Δ𝑢′ 1 (𝑧)

[

∫Ω

(50)

𝑩 (𝐼𝐿1) (𝑧)

where n is the number of nodes with domains of influence which cover point z, ( ( ) ( ) ( ))T ∗ Δ𝒖′ = Δ𝑢′ 𝑧1 , Δ𝑢′ 𝑧2 , … , Δ𝑢′ 𝑧𝑛 . (39)

∼

T ̂ dΩ + 𝑩 (𝐿0) 𝑺

̃ T 𝒕dΓ, 𝚽

(52)

Combining the weak integral form in Eq. (16) with the IICVMLS method, the IICVEFG method for the two-dimensional large deformation of hydrogel can be obtained. Decorate M nodes zI (I=1,2,...,M) in the problem domain Ω. In addition, their influence domains ΩI must cover the whole domain Ω. According to Eq. (35), for an arbitrary point z in the domain, the displacement increment Δu′(z) is written as

where

∫Ω

∫Γ𝑡

[ ] [ ] ⎡ Re Φ𝐼,1 (𝑧) − Im Φ𝐼,1 (𝑧) [ ] [ ] ⎢ 𝑩 (𝐼𝐿0) (𝑧) = ⎢ Im Φ𝐼,2 (𝑧) Re Φ𝐼,2 (𝑧) [ ] [ ] [ ] [ ] ⎢ ⎣ Re Φ𝐼,2 (𝑧) + Im Φ𝐼,1 (𝑧) − Im Φ𝐼,2 (𝑧) + Re Φ𝐼,1 (𝑧)

4. Disperse process

[

̃ T 𝒃dΩ + 𝚽

(37)

where 𝜌 is the radius of the influence domain of the point z. The above is brief introduction of the IICVMLS method. It is obvious that compared with other non-interpolating CVMLS approximations, the shape function in the IICVMLS method can satisfy the property of Kronecker 𝛿 function.

Δ𝒖′ (𝑧) = Δ𝑢′1 (𝑧) + iΔ𝑢′2 (𝑧) =

∫Ω

In this section, a simply constrained swelling deformation of hydrogel is simulated to validate the advantages of the IICVEFG method. Then the IICVEFG method is used to analyze the pattern transformation of hydrogels. Nine samples with different shapes and initial sizes of holes in a square hydrogel are designed to show their effects on the pattern transformation. The specific material parameters, boundary conditions, node distributions and other computing parameters are shown below.

(46)

(47)

102

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

x2

1

x1 O 1 Fig. 1. A square hydrogel with three constrained boundary and node distribution.

Fig. 3. Relative error of stretch 𝜆2 . 1600 1400

IICVEFG ICVEFG

CPU time (s)

1200 1000 800 600 400 200 0

7*7

9*9

11*11

13*13

15*15

Node distributions Fig. 2. Stretch 𝜆2 with different chemical potential 𝜇/kT.

Fig. 4. Computing time under different node distributions.

is a unique value which agrees well with the analytical solution. However, in the ICVEFG method it is necessary to choose appropriate penalty factors to apply the essential boundary conditions, and good numerical solutions can be obtained when the penalty factor 𝛼 is within the range of 102 ~ 105 . Fig. 3 displays the relative errors of the stretch 𝜆2 in both IICVEFG and ICVEFG methods. Compared with the ICVEFG method, the IICVEFG method has smaller error under different chemical potential. Fig. 4 shows the computing time of the IICVEFG method and ICVEFG method with the chemical potential 𝜇/kT = 0.0 under the node distributions 7 × 7, 9 × 9, 11 × 11, 13 × 13, 15 × 15, respectively. In the ICVEFG method, the penalty factor 𝛼 = 103 . It is clear that the CPU time using IICVEFG method is lower than that using the ICVEFG method. According to above discussions, the IICVEFG method has higher computing efficiency and accuracy than the ICVEFG method due to the essential boundary conditions can be applied directly and the use of the complete basis function.

5.1. Validation As shown in Fig. 1, a square hydrogel with three constrained boundaries is immersed in a solvent with chemical potential 𝜇/kT. The side length is 1mm and other material parameters are: Nv = 0.001 and 𝜒 = 0.1. This problem is regarded as a plane strain problem. 49 nodes are distributed uniformly in the square hydrogel and the 6 × 6 integral cells are used in the simulation process. In both the IICVEFG and ICVEFG methods, the scaling parameter is dmax = 3.5. Because this square hydrogel is only constrained in x1 direction, the hydrogel can swell freely along x2 direction and the analytical solutions of the swelling stretch 𝜆2 under different chemical potential 𝜇/kT can be obtained according to the following formula by setting the stress along x2 direction to be zero [22], ( ) ( ) 𝜆2 − 1 𝜒 𝜇 1 1 𝑁𝑣 𝜆2 − + ln + + − = 0. (59) 𝜆2 𝜆2 𝜆2 𝜆2 𝑘𝑇 2

In Fig. 2, the numerical solutions of the stretch 𝜆2 obtained from the IICVEFG method and non-interpolating ICVEFG method are compared with the analytical solutions under different chemical potential 𝜇/kT. When the node distribution and scaling parameter dmax are determined, the numerical solution of stretch 𝜆2 in the IICVEFG method

5.2. Sample design Fig. 5 shows the size and boundary constraints of a square hydrogel model which is soaked in a solvent with an alterable chemical potential 𝜇/kT. The dimension of this hydrogel model is 22 mm × 22 mm and the 103

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Nine samples are constructed by excavating holes in the square hydrogel mdoel as shown in Table 1. There are three types of holes with different shapes and sizes. Samples 1, 2 and 3 are antisymmetrical structures with four specially arranged elliptical holes, and other samples are symmetrical structures with four circles and squares, respectively. In these samples, the region between two adjacent holes is called the vertical pillar as shown in Table 1. In order to analyze the effects of shape of holes on the pattern transformation reasonably, the samples 1, 4 and 7 are chosen with same thinnest spaces of the vertical pillar and the width is 1 mm. Then the effects of initial size of holes in the samples with same shape of holes are analyzed by fixing the center of holes in the same place (like samples 1, 2 and 3).

x2

x1

22

O

5.3. Simulation scheme When using the IICVEFG method to simulate the transformation behavior of the samples in the environment of solvent and geometric constraints, these large deformation problems are considered as the plane strain problems. The scaling parameter is chosen as dmax = 3.5, and the space of the adjacent nodes is set to be same and uniform in different samples. 4 × 4 Gauss points are used for the Gauss quadrature in each background integration. Around the elliptical, circular and square holes, 40 nodes are arranged. The numbers of total distributed nodes and Gauss points are listed in Table 1.

22 Fig. 5. A square hydrogel model with four constrained boundaries.

material parameters are: Nv = 0.001 and 𝜒 = 0.0001. Four edges of the model are constrained by the rolling hinged supports and no external force boundary condition is applied.

Table 1 The designed samples with the corresponding total numbers of nodes and Gauss points.

vertical pillar b (a) Ellipse Sample 1 Sample 2 Sample 3

a(mm) 10 10 8

b(mm) 8 6 5

Nodes 1009 1213 1417

Gauss points 12,284 15,600 18,936

a

r (b) Circle Sample 4 Sample 5 Sample 6

r(mm) √ √20 √15 10

Nodes 989 1225 1437

Gauss points 12,312 15,624 18,956

a (c) Square Sample 7 Sample 8 Sample 9

a(mm) √ √20𝜋 √15𝜋 10𝜋

Nodes 941 1165 1441

Gauss points 12,114 15,196 18,876

104

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Fig. 6. Pattern transformations of sample 1 with different 𝜇/kT.

Sample1

Sample2

Sample3

Fig. 7. Pattern transformations of samples 1, 2 and 3 when 𝜇/kT = −0.9178.

Fig. 8. Pattern transformations of sample 4 with different 𝜇/kT.

5.4. Pattern transformation rule

Fig. 8 shows the deformation process of the sample 4 with increasing chemical potential. The circular holes still keep initial shapes with the chemical potential 𝜇/kT = −4.9678. Then when the chemical potential 𝜇/kT = −0.6678, the circles in the sample 4 are deformed into lozenge holes. This deformation mode is consistent with Okumura’s research [21]. With the chemical potential increasing continuously, the four sides of circular holes produce concave deformation and the circular holes become almost square holes when 𝜇/kT = −0.2678 as shown in Fig. 8. Fig. 9 observes the deformations of the samples 4, 5 and 6 with different initial size of circular holes under the same chemical potential 𝜇/kT = −0.9500. It is obvious that the vertical pillar produces a symmetrical deformation and the center line of the vertical pillar keeps straight. Compared with the samples 4 and 5, the deformation of vertical pillar

Fig. 6 shows the transformation of the sample 1 under the influence of the chemical potential 𝜇/kT. When the chemical potential 𝜇/kT = −4.9678, the elliptical holes have no obvious deformation. Then with the increase of the chemical potential, the elliptical holes become specific shapes as shown in Fig. 6, and the square hydrogel can produce a certain range of anticlockwise rotation due to the special arrangement of elliptical holes, which agrees well with the experimental description [16]. Fig. 7 reveals the pattern transformations of the samples 1, 2 and 3 with different initial sizes of elliptical holes when the chemical potential 𝜇/kT = −0.9178. We can see that the vertical pillar produces the buckling behavior. And in the sample 3, the deformation of vertical pillar is greater than that in the samples 1 and 2. 105

Y. Deng, X. He and L. Sun et al.

Sample4

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Sample5

Sample6

Fig. 9. Pattern transformations of samples 4, 5 and 6 when 𝜇/kT = −0.9500.

Fig. 10. Pattern transformations of sample 7 with different 𝜇/kT.

Sample7

Sample8

Sample9

Fig. 11. Pattern transformations of samples 7, 8 and 9 when 𝜇/kT = −0.9178.

in the sample 6 is greater, especially in the region near the circularxbrk holes. The deformation process of the sample 7 and the corresponding chemical potentials are displayed in Fig. 10. When the chemical potential increases from 𝜇/kT = −4.9678 to 𝜇/kT = −0.5678, the four sides of the square holes deform inward into the shapes as shown in Fig. 10. Such deformation behavior is similar to the experimental result of Yang et al [17]. Then with the chemical potential keeps increasing, the sections near the boundaries produce larger deformation than the sections near the center lines along x axis and y axis, because the distances between the square holes and the boundaries are greater than the distances between the square holes and the center lines. Under the chemical potential 𝜇/kT = −0.9178, Fig. 11 compares the deformations of the samples 7, 8 and 9 with different initial sizes of the

square holes. Similar to the hydrogel with circular holes, the vertical pillar yields a symmetrical deformation and the center line of the vertical pillar in these three samples also remains straight. Because the width of the vertical pillar increases from the samples 7–9, the deformation of four sides in the square holes of the sample 9 are larger than that of the samples 7 and 8. 5.5. Actuating feature It is seen from the above deformation modes that the samples with specific distributed elliptical holes will produce obvious rotation. The change of rotation angle with the increase of chemical potential in the samples 1, 2 and 3 is illustrated in Fig. 12. The midline along the x direction is chosen to analyze the rotation degree during the deformation 106

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Fig. 12. Rotation angle with different chemical potentials in the samples 1, 2 and 3.

Fig. 13. Rotation behavior of square hydrogel with different 𝜇/kT.

As we can see, with the increase of the chemical potential from 𝜇/kT = −4.9688 to 𝜇/kT = −0.2688, the circular holes become diamond holes and there is no obvious rotation. But when the 𝜇/kT = −0.1688 with larger deformation stretch 𝜆 = 1.4510, the area between holes produce buckling and the rotation behavior is produced.

process. The angle between the x coordinate and the line connecting the highest point and the base point in the midline is defined as the rotation angle 𝜃. As can be seen from Fig. 12, the deformation process can be divided into two stages [44]. (I) When the chemical potential changes within a range from 𝜇/kT = −∞ to 𝜇/kT = −3.9678, there is no obvious rotation. In this stage, the volume and shape of hydrogel have no change. (II) After 𝜇/kT = −3.9678, the increasing chemical potential leads to a remarkable change of the rotation angle in the square structure with elliptical holes, especially when the chemical potential is close to 0. The volume and shape of hydrogel are also changed. This stage can be regarded as a phase transition process [48,49]. As shown in Fig. 12, the rotation angles in samples 1, 2 and 3 are approximately equal to 30°, 20° and 18°, respectively. In the sample 1, the rotation angle 30° of point A agrees well with the experimental result obtained by Yang et al [16]. And the rotation angle decreases when reducing the initial size of elliptical holes. Compared with these three samples, other samples have no obvious rotation of the midline along the x direction because they are symmetrical structures. However, when the symmetrical hydrogel structure is designed reasonably with different distribution of holes and constrained boundary regions, it also can produce rotation behavior. A square hydrogel structure with circular holes is considered as an example shown in Fig. 13. The constrained boundary regions are changed to the surroundings of the holes compared with above samples, but the constrained way is still rolling hinged supports. The dimension of this structure is (8 × 10-3 mm) × (8 × 10-3 mm). There are nine intact circular holes distributed in the square hydrogel and the radius of the circular holes is 0.9 × 10-3 mm. The material parameters are Nv = 0.001 and 𝜒 = 0.00001. 840 nodes and 2304 Gauss points are uniformly distributed in the square hydrogel. In the IICVEFG method, the scaling parameter is dmax = 3.5.

5.6. Area shrinkage With the increase of chemical potential, the swelling stretch 𝜆 will increase and the holes in the square hydrogel will produce shrinkage behavior. Fig. 14(a) shows the relationship between the swelling stretch 𝜆 and the change of area of each hole in different samples. It is observed that when the initial area of each hole decreases from S = 20𝜋 mm2 to S = 10𝜋 mm2 , the corresponding swelling stretch 𝜆 is reduced, which means that the affordable chemical potential is reduced for the samples. As shown in Fig. 14(b), with the reduction of the initial size of holes, the percentage of area reduction increases in the samples with elliptical holes but decreases in the samples with circular and square holes. When S = 20𝜋 mm2 , comparing the samples 1, 4 and 7 which have same initial size of holes and vertical pillar but different shapes of holes, the sample 7 with square holes has the maximum reduction of area of 77%, and the sample 1 with elliptical holes has minimum reduction of area. For the sample 4 with circular holes, the reduction of area of 60% can be achieved. Before concluding, it is worth to point out that the square hydrogel can be designed properly with specific distributed holes to achieve the rotation and area shrinkage behaviors. By using the rotation behavior, the square hydrogel can be regarded as an actuator to obtain and control the expected motions. As for the area reduction, the square hydrogel has the potential to control flow as a soft valve. And the transparency of the 107

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

Fig. 14. Area reduction with the increase of swelling stretch 𝜆.

square hydrogel will change due to the area reduction of holes, as shown in Zhu’s experiment, which can be used in the governable window [15].

nicipality (JCYJ20160229165310679), and the National Water Pollution Control and Treatment Science and Technology Major Project (2017ZX07201001).

6. Conclusions References Through combining the IICVMLS method with the basic field theory of hydrogel, an effective IICVEFG method is proposed and the final discrete equilibrium equation is achieved for the large deformation analysis of hydrogels. In contrast to the non-interpolating meshless methods, the shape function in the IICVMLS method satisfies the property of Kronecker 𝛿 function and the IICVEFG method can apply the essential boundary conditions directly. Therefore, the expressions of matrices in the IICVEFG method are more concise and there is no need to choose appropriate penalty factor, resulting in higher accuracy and efficiency. Through analyzing the transformation of nine samples of a square hydrogel, the samples with different shapes and initial sizes of holes have diverse deformation patterns with the increase of chemical potential. When under the same chemical potential, compare the samples with same shape of holes, the vertical pillar produces larger deformation in the sample with smaller initial size of holes. The sample with special arranged elliptical holes can produce large rotation angle of 30° with low area reduction rate. While by reducing the initial size of elliptical holes, lower rotation angle and larger area reduction are obtained for the hydrogel. Being different from the samples with elliptical holes, the samples with circular and square holes have no rotation behavior, and when the initial size of holes decreases, the corresponding holes have less area reduction. By utilizing above findings, the hydrogel has the potential to actuate the soft machines, control the flow as a soft valve and make a transparency adjustable window if the researchers design reasonable shapes and initial sizes of holes in a square hydrogel.

[1] Ahmed EM. Hydrogel: preparation, characterization, and applications: a review. J Adv Res 2015;6:105–21. [2] De SK, Aluru NR. A chemo-electro-mechanical mathematical model for simulation of pH sensitive hydrogels. Mech Mater 2004;36:395–410. [3] Yang QS, Qin QH, Ma LH, Lu XZ, Cui CQ. A theoretical model and finite element formulation for coupled thermo-electro-chemo-mechanical media. Mech Mater 2010;42:148–56. [4] Lam KY, Li H, Ng TY, Luo R. Modeling and simulation of the deformation of multi-state hydrogels subjected electrical stimuli. Eng Anal Bound Elem 2006;30:1011–17. [5] Liu QM, Li H, Lam KY. Optimization of the cell microenvironment in a dual magnetic-pH-sensitive hydrogel-based scaffold by multiphysics modeling. Bioelectrochemistry 2019;129:90–9. [6] Wu T, Li H. Phase-field model for liquid-solid phase transition of physical hydrogel in an ionized environment subject to electro-chemo-thermo-mechanical coupled field. Int J Solids Struct 2018;138:134–43. [7] Chung HJ, Park TG. Self-assembled and nanostructured hydrogels for drug delivery and tissue engineering. Nano Today 2009;4:429–37. [8] Harmon ME, Tang M, Frank CW. A microfluidic actuator based on thermo responsive hydrogels. Polymer (Guildf) 2013;44:4547–56. [9] He T, Li M, Zhou JX. Modeling deformation and contacts of pH sensitive hydrogels for microfluidic flow control. Soft Matter 2012;8:3083–9. [10] Smela E. Conjugated polymer actuators for biomedical applications. Adv Mater 2003;15:481–94. [11] Singamaneni S, Bertoldi K, Chang S, Jang JH, Young SL. Bifurcated mechanical behavior of deformed periodic porous solids. Adv Funct Mater 2009;19:1426–36. [12] Singamaneni S, Tsukruk VV. Buckling instabilities in periodic composite polymeric materials. Soft Matter 2010;6:5681–92. [13] Zhang Y, Matsumoto EA, Peter A, Lin PC, Kamien RD. One-step nanoscale assembly of complex structures via harnessing of an elastic instability. Nano Lett 2008;8:1192–6. [14] Dortdivanlioglu B, Linder C. Diffusion-driven swelling-induced instabilities of hydrogels. J Mech Phys Solids 2019;125:38–52. [15] Zhu X, Wu G, Dong R, Chen CM, Yang S. Capillarity induced instability in responsive hydrogel membranes with periodic hole array. Soft Matter 2012;8:8088–93. [16] Yang D, Mosadegh B, Ainla A, Lee B, Khashai F, Suo ZG, et al. Buckling of elastomeric beams enables actuation of soft machines. Adv Mater 2015;27:6323–7. [17] Yang D, Verma MS, So J, Mosadegh B, Keplinger C, Lee B, et al. Buckling pneumatic linear actuators inspired by muscle. Advanced Materials Technologies 2016;1:1600055.

Acknowledgments The authors would like to acknowledge financial support from the Science and Technology Innovation Commission of Shenzhen Mu108

Y. Deng, X. He and L. Sun et al.

Engineering Analysis with Boundary Elements 113 (2020) 99–109

[18] Mullin T, Deschanel S, Bertoldi K, Boyce MC. Pattern transformation triggered by deformation. Phys Rev Lett 2007;99:084301. [19] Bertoldi K, Boyce MC, Deschanel S, Prange SM, Mullin T. Mechanics of deformation triggered pattern transformations and superelastic behavior in periodic elastomeric structures. J Mech Phys Solids 2008;56:2642–68. [20] Okumura D, Inagaki T, Ohno N. Effect of prestrains on swelling induced buckling patterns in gel films with a square lattice of holes. Int J Solids Struct 2015;58:288–300. [21] Okumura D, Kuwayama T, Ohno N. Effect of geometrical imperfections on swelling-induced buckling patterns in gel films with a square lattice of holes. Int J Solids Struct 2014;51:154–63. [22] Li DM, Zhang Z, Liew KM. A numerical framework for two-dimensional large deformation of inhomogeneous swelling of gels using the improved complex variable element-free Galerkin method. Comput Method Appl Math 2014;274:84–102. [23] Liu FB, Cheng YM. The improved element-free Galerkin method based on the nonsingular weight functions for inhomogeneous swelling of polymer gels. Int J Appl Mech 2018;10:1850047. [24] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Method Appl Math 1996;139:3–47. [25] Peng MJ, Li RX, Cheng YM. Analyzing three-dimensional viscoelasticity problems via the improved element-free Galerkin (IEFG) method. Eng Anal Bound Elem 2014;40:104–13. [26] Cheng H, Peng MJ, Cheng YM. Analyzing wave propagation problems with the improved complex variable element-free Galerkin method. Eng Anal Bound Elem 2018;100:80–7. [27] Li XL, Li SL. Analyzing the nonlinear p-Laplacian problem with the improved element-free Galerkin method. Eng Anal Bound Elem 2019;100:48–58. [28] Belytschko T, Lu YY, Gu L. Element free Galerkin methods. Int J Numer Meth Eng 1994;37:229–56. [29] Bai FN, Li DM, Wang JF, Cheng YM. An improved complex variable element free Galerkin method for two-dimensional elasticity problems. Chinese Phys B 2012;21:020204. [30] Liew KM, Feng C, Cheng YM, Kitipornchai S. Complex variable moving least squares method: a meshless approximation technique. Int J Numer Meth Eng 2007;70:46–70. [31] Cheng YM, Li RX, Peng MJ. Complex variable element free Galerkin method for viscoelasticity problems. Chinese Phys B 2012;21:090205. [32] Peng MJ, Liu P, Cheng YM. The complex variable element free Galerkin (CVEFG) method for two-dimensional elasticity problems. Int J Appl Mech 2009;1:367–85.

[33] Zhang LW, Deng YJ, Liew KM, Cheng YM. The improved complex variable element-free Galerkin method for two-dimensional Schrödinger equation. Comput Math Appl 2014;68:1093–106. [34] Cheng YM, Wang JF, Bai FN. A new complex variable element-free Galerkin method for two dimensional potential problems. Chinese Phys B 2012;21:090203. [35] Zhu T, Atluri SN. A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free Galerkin method. Comput Mech 1998;21:211–22. [36] Lancaster P, Salauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37:141–58. [37] Kaljevi I, Saigal S. An improved element free Galerkin formulation. Int J Numer Meth Eng 1997;40:2953–74. [38] Ren HP, Cheng YM. The interpolating element free Galerkin (IEFG) method for two-dimensional elasticity problems. Int J Appl Mech 2011;3:735–58. [39] Ren HP, Cheng YM. The interpolating element free Galerkin (IEFG) method for two dimensional potential problems. Eng Anal Bound Elem 2012;36:873–80. [40] Deng YJ, Liu C, Peng MJ, Cheng YM. The interpolating complex variable element-free Galerkin method for temperature field problems. Int J Appl Mech 2015;7:1550017. [41] Deng YJ, He XQ. An improved interpolating complex variable meshless method for bending problem of Kirchhoff plates. Int J Appl Mech 2017;9:1750089. [42] Hong W, Liu ZS, Suo ZG. Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load. Int J Solids Struct 2009;46:3282–9. [43] Galante S, Lucantonio A, Nardinocchi P. The multiplicative decomposition of the deformation gradient in the multiphysics modeling of ionic polymers. Int J Nonlin Mech 2013;51:112–20. [44] Hong W, Zhao XH, Zhou JX, Suo ZG. A theory of coupled diffusion and large deformation in polymeric gels. J Mech Phys Solids 2008;56:1779–93. [45] Liu ZS, Toh W, Ng TY. Advances in mechanics of soft materials: a review of large deformation behavior of hydrogels. Int J Appl Mech 2015;7:1530001. [46] Flory PJ. Thermodynamics of high polymer solutions. J Chem Phys 1942;10:51–61. [47] Flory PJ, Rehner J. Statistical mechanics of cross linked polymer networks II: swelling. J Chem Phys 1943;11:521–6. [48] Ding ZW, Liu ZS, Hu JY, Swaddiwudhipong S, Yang ZZ. Inhomogeneous large deformation study of temperature sensitive hydrogel. Int J Solids Struct 2013;50:2610–19. [49] Yang D, Jin LH, Martinez RV, Bertoldi K, Whitesides GM, Suo ZG. Phase transforming and switchable metamaterials. Extreme Mech Lett 2016;6:1–9.

109