Finite element analysis of auxetic obstacle deformation and fluid flow in a channel

Journal of Non-Crystalline Solids 355 (2009) 1387–1392

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/locate/jnoncrysol

Finite element analysis of auxetic obstacle deformation and fluid flow in a channel T. Strek a, P. Kedziora b, B. Maruszewski a,*, A.A. Pozniak b, K.V. Tretiakov b, K.W. Wojciechowski b,* a b

Poznan University of Technology, Institute of Applied Mechanics, ul. Piotrowo 3, 60-965 Poznan, Poland Institute of Molecular Physics, Polish Academy of Sciences, ul. M. Smoluchowskiego 17, 60-179 Poznan, Poland

a r t i c l e

i n f o

Article history: Available online 17 June 2009 PACS: 62.20.-x 62.20.dj 47.10.ad 47.11.Fg

a b s t r a c t The influence of the Poisson’s ratio of a viscoelastic obstacle on the flow in a two-dimensional channel was studied by a finite element computer simulation. The arbitrary Lagrangian–Eulerian (ALE) formulation for deformed mesh was applied to handle the moving boundaries of the computational domain. An anomalous deformation of the obstacle was observed when the Poisson’s ratio was negative and close to its minimum value. It was also found that the maximum velocity of the fluid in the system decreased with the decreasing Poisson’s ratio. Ó 2009 Elsevier B.V. All rights reserved.

Keywords: Auxetics Fluid-solid interaction Mechanical properties

1. Introduction Materials for which the Poisson’s ratio [8] is negative were first manufactured by Lakes [7] and later named auxetics by Evans [6]. Typically, the Poisson’s ratio of isotropic materials is positive [8]. Although a negative Poisson’s ratio has not been observed in most everyday materials, the classical theory of elasticity states that it is possible for isotropic three-dimensional materials to exhibit the Poisson’s ratio within the range of 1  m  þ1=2 [8] whereas two-dimensional isotropic systems can exhibit Poisson’s ratios within the range of 1  m  þ1 [15]. In anisotropic materials Poisson’s ratios can have any positive or negative value [13]. Recently, increasing interest has been shown in the field of materials with a negative Poisson’s ratio [11] which are interesting from both an intellectual point of view as well as in terms of potential applications. Studies of simple model systems exhibiting a negative Poisson’s ratio or including auxetic parts should help recognize possible unconventional properties of such systems that may result in finding new potential applications of auxetics. The problem studied in this paper concerns a horizontal flow of an incompressible viscous fluid along a channel with an auxetic obstacle formed by a thin plate normal to its walls (Fig. 1). The main aim of our considerations is to find how the obstacle behaves under the flow. For the sake of simplicity, a version of the problem is considered within a plane stress state [3,4,9].

* Corresponding authors. E-mail addresses: [email protected] (T. Strek), [email protected] (B. Maruszewski), [email protected] (K.W. Wojciechowski). 0022-3093/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2009.05.032

Exact analytical treatment of the problem’s mechanical model, even in its simplest version considered in this paper, is very difficult. Fortunately, solid–fluid interactions can be modeled by approximate finite element method (FEM) techniques [16]. In particular, the arbitrary Lagrangian–Eulerian (ALE) method can be used to follow the FEM mesh geometry deformation resulting from the motion of the computational domain boundaries [4,5]. The technique is very useful for a simulation of solid–fluid mechanics. We applied this approach to improve the accuracy of the mathematical description of our problem. This requires a few words of comment. In the special case of the Lagrangian method, the mesh motion follows the motion of the solid material. However, when the material motion is more complicated, as in a fluid flow case, that method is not appropriate. In such a situation the Eulerian method, where the mesh is fixed, is often used even though this method cannot account for the moving boundaries. Therefore, we used the combined ALE technique to describe the fluid flow with detailed tracking of the obstacle boundaries together with the obstacle deformation. 2. Mathematical model Since the fluid flow along a channel with an obstacle concerns two interacting media, a mathematical description of the mechanics of the entire structure has to consist of two parts. The first part concerning the flow domain is described by an initial-boundary value problem based on the Navier–Stokes equation. The second part concerns the solid obstacle which is governed by the viscoelastic equation of motion with proper boundary and initial conditions.

1388

T. Strek et al. / Journal of Non-Crystalline Solids 355 (2009) 1387–1392

The initial condition for (1) is based on the assumption that the fluid at t ¼ 0 is motionless (v ¼ 0) and p ¼ 0. The equation of motion for the solid obstacle reads:

qs

  @2u @u  r  r ¼ G; þ n @t @t 2

ð2Þ

where qs is the obstacle density, u is the displacement vector, nð@u Þ @t is the damping term, r is the stress tensor and G is the body force.

Fig. 1. The geometry of the studied system.

As has been mentioned above, the fluid flow in the channel is described by the Navier–Stokes equations:

q0

  @v þ v  rv ¼ rp þ r  S þ F; @t

ð1Þ

r  v ¼ 0; where F is the volume force affecting the fluid and v is the fluid q0 is  the fluid density, p is the pressure, velocity,  S ¼ g rv þ ðrv ÞT is the extra stress tensor coming from fluid viscosity and g is the dynamic viscosity. In the studied case the force affecting the fluid has been omitted, i.e. F ¼ 0. The boundary conditions resulting from the stated problem are the following: 

v ¼ 0 at the lower and upper channel surface because there is no



min ¼ 4m0 hy ðhy  1Þ for y 2 h0; hi, where h is the obstacle height

slip there, (inlet parabolic profile condition),  ðpI þ SÞ  n ¼ p0 n, where p0 ¼ 0 is the outlet pressure,  v ¼ 0 at the clamped obstacle surfaces (no slip condition), at the deformed obstacle surfaces (u is the displacement  v ¼ @u @t vector of the obstacle point).

Table 1 The parameters used in simulations. Quantity Symbol

Unit

Value

q0 g qs

kg/m3 Pa s kg/m3 kPa – s1 S m/s

1000 103 7850 20 0.99, 0.7, 0.3, +0.3 1 103 1.0

E

m adM bdK

m0

Table 2 Mesh data. Quantity

Mesh A

Mesh B

Mesh C

Global number of mesh points Global number of mesh elements

492 952

772 1491

2332 4528

Fig. 2. Total displacement in a deformed auxetic obstacle for Poisson’s ratio m ¼ 0:99 with meshes: (A) A, (B) B, and (C) C from Table 2.

1389

T. Strek et al. / Journal of Non-Crystalline Solids 355 (2009) 1387–1392

Since the material of the obstacle has been assumed to be isotropic, the following constitutive equation writes:

r ¼ De ¼ kðr  uÞI þ 2le;

ð3Þ T

where I is the unit matrix, e ¼ 12 ðru þ ðruÞ Þ is the strain tensor and k and l denote the Lamé constants as follows:

k¼

Em ; ð1 þ mÞð1  2mÞ

l¼

E ; 2ð1 þ mÞ

where E is Young’s modulus, elasticity matrix,

2 6 6 6 6 E 6 D¼ ð1 þ mÞð1  2mÞ 6 6 6 4

1m

m

m m

1m

m m

m

0

0

0 0

0 0

0 0

1m

0

0

0

12m 2

0

0

0

0

12m 2

0

0

0

0

ð4Þ

m is the Poisson’s ratio and D is the

3 0 0 7 7 7 0 7 7: 0 7 7 7 0 5 12m 2

ð5Þ and in the same notation as in (5) the stress and strain read:

Fig. 3. Total displacement (A) and first principal strain (B) in a deformed auxetic obstacle (m ¼ 0:99) after time t = 0.02 s.

Fig. 4. Velocity field of fluid in a deformed domain at the steady state: (A)

m ¼ 0:99 (B) m ¼ 0:7 (C) m ¼ 0:3 (D) m ¼ 0:3 for v0 = 1 m/s.

1390

T. Strek et al. / Journal of Non-Crystalline Solids 355 (2009) 1387–1392

2

rx

3

6r 7 6 y7 7 6 6 rz 7 7 r¼6 6s 7 6 xy 7 7 6 4 syz 5

2

ex

3

6 ey 7 7 6 7 6 6 ez 7 7 e¼6 6 c 7: 6 xy 7 7 6 4 cyz 5

sxz

2l k Em2 l ¼ l; k ¼ ¼ : k þ 2l 1  m2 ð6Þ

cxz

f ðx; y; tÞ ¼ 1 H

Z

H

f ðx; y; z; tÞdz;

ð7Þ

3

r x 6 7  r ¼ 4 ry 5; sxy

2

¼ D

 1  þ ðru  ÞT ; ru 2

ð11Þ

we arrive at the final equation of motion in terms of (7) [3,4]:

qs H

   @2u @u @u  r  c  ru þ c  bdK r þ adM qs H ¼ 0; 2 @t @t @t

ð12Þ

where is c the flux matrix:

" c¼

E 1m2

H

0

#

0 E H 2ð1þmÞ

ð13Þ

:

0

where f stands for each variable which should be averaged according to the plane stress mode. As the results one has,

2

Neglecting now the body force in (2), assuming the Rayleigh damping in the obstacle [10,14] and its small deformations in the form:

e ¼

Looking at the considered model geometry (Fig. 1) one can note that when H tends to infinity, in each plane x–y along z-direction the same mechanical situation may be assumed. This means that all z components of the stress are equal to zero what encourages one to describe the problem within a two-dimensional (2D) plane stress state. For finite H, Eqs. (2)–(5) are transformed to that state by using the following averaging operation [9]:

ð10Þ

ð8Þ 1

m

E 6 4m 1 1  m2 0 0

0 1m 2



3

0 7 5:

Two damping coefficients have been introduced in (12): adM is the mass damping parameter and bdK is the stiffness damping parameter [10,14]. H is the width of the obstacle and channel in @ @ ex þ @y ey . z-direction (Fig. 1(A)), then r ¼ @x The boundary conditions for (12) result from the fact that all free obstacle boundaries are loaded by the fluid as follows:



r  n ¼  pI þ g rv þ ðrv ÞT



 n;

ð14Þ

ð9Þ where n is the normal vector to the obstacle boundary. This load represents a sum of pressure and viscous forces. The edge load

Moreover [9]:

Fig. 5. Total displacement in a deformed obstacle at the steady state: (A)

m ¼ 0:99 (B) m ¼ 0:7 (C) m ¼ 0:3 (D) m ¼ 0:3 for v0 = 1 m/s.

T. Strek et al. / Journal of Non-Crystalline Solids 355 (2009) 1387–1392

Fig. 6. First principal strain in a deformed obstacle near the fixed boundary: (A)

has been defined as a force per area taking the thickness into account. The initial conditions for (12) read: the displacement u ¼ 0 and there is no initial stress, strain and body load on the obstacle. 3. Numerical simulations Following Fig. 1, the detailed geometry of the problem is the following:  the flow channel is h ¼ 100 lm high (y-direction) and l ¼ 300 lm long (x-direction),  the flank surfaces (upper and lower) of the obstacle are clamped,  the obstacle dimensions: 5 lm width (x-direction), 50 lm height (y-direction),  the obstacle is placed 70 lm away from the channel left boundary,  the system width is H ¼ 107 m (z-direction). The remaining data have been collected in Table 1. The flow and deformation of the structure were investigated within the time interval between 0 and 0.02 s. The problem was solved using the COMSOL code with the solver DASPK [4,1,2] and the direct UMFPACK linear system solver for FEM matrices [4] with 4528 Lagrange quadratic triangular mesh elements. It should be noted that the total displacement and first principal strain values in all the considered examples are small comparing to the obstacle width. Since the obstacle deformation is very small we

1391

m ¼ 0:99, (B) m ¼ 0:7, (C) m ¼ 0:3, (D) m ¼ 0:3 for v0 = 1 m/s.

(artificially) scaled up the deformation magnitude on the figures to see the deformation details in the vicinity of the plate’s end. To confirm the reliability of the results we performed the finite element analysis with different meshes. The number of mesh points and the number of mesh elements are presented in Table 2. It can be seen in Fig. 2(A)–(C) that the results are practically indistinguishable for all the three fine meshes. Figs. 5 and 6 show in detail only the upper part of the deformed obstacle. All Figs. 3–6 were drawn based on the FEM calculation results using 1586 mesh elements (mesh B). Fig. 3(A) and (B) present the deformed obstacle after flow for t ¼ 0:02 s. Fig. 3A presents the total displacement in a deformed auxetic (m ¼ 0:99) obstacle. Fig. 3(B) presents the first principal strain of the obstacle. It is easy to show that the time-dependent flow with a constant inlet parabolic velocity profile reaches a steady-state after time t = 0.02 s. Figs. 4–6 present the numerical results for the flow with maximum inlet velocity v 0 ¼ 1 m/s. Fig. 4 presents the steady state velocity field and the deformed obstacles with different Poisson’s ratio value (m ¼ 0:99; 0:7; 0:3 and 0:3). Fig. 5 presents a total displacement in the deformed obstacle. Fig. 6 presents the first principal strain in the considered obstacle. Let us still consider how Young’s modulus values are connected with the elastic properties of the obstacle both for classical and auxetic states. In Fig. 7, we have shown a situation where that modulus increases for the material with the Poisson’s ratio m ¼ 0:99.

1392

T. Strek et al. / Journal of Non-Crystalline Solids 355 (2009) 1387–1392

Fig. 7. First principal strain of an auxetic obstacle vs. Young’s modulus (A) E = 5  104 N/m2, (B) E = 105 N/m2, (C) E = 106 N/m2, (D) E = 107 N/m2.

4. Conclusions

Acknowledgement

A finite element analysis of fluid–solid interactions was applied to study the deformations of an auxetic obstacle in a 2D channel flow. It was found that the displacement of the obstacle points coming up to the clamped edges along the plate vanished faster than for the case of a non-auxetic material (cf. Figs. 2, 3(A), 5(A) and (D)). Figs. 3(B) and 6(A) indicate that there is practically no deformation of the auxetic obstacle in the regions directly close to its clamped edges contrary to the situation for the non-auxetic material shown in Fig. 6(D). A surprising feature of the auxetic plate can be seen in Figs. 2, 3(A) and 5(A) where the displacement vector firstly turns oppositely during the auxetic obstacle deformation (growing away from the clamped edges) and then it is in the direction of the fluid flow contrary to the non-auxetic case (Fig. 5(D)), where such a turn was not observed and the displacement vector had always the same direction as the fluid flow. Similar behavior was observed when analyzing the deformation of a three-dimensional plate [12]. This effect is most significant in materials of low values of Young’s modulus. Finally, the obstacle with a negative Poisson’s ratio influences the fluid velocity field in such a way that the values of velocity in the vicinity of the channel walls are lower compared to the classical (non-auxetic) case (Fig. 4). This means that the probability of a turbulence or even cavitation decreases when an auxetic obstacle is used.

A part of this work was supported by Grant N202 070 32/1512 (2007-2010) of the Ministry of Sciences and Higher Education. References [1] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, 2nd Ed., Elsevier, New York, 1989. SIAM, 1996. [2] P.N. Brown, A.C. Hindmarsh, L.R. Petzold, SIAM J. Sci. Comput. 15 (1994) 1467. [3] Cosmol Structure Mechanics Module, FEMLAB 3.1 User’s Guide, Comsol, AB, 2004. [4] Comsol Multiphysics 3.3 User’s Guide, Comsol AB, 2006. [5] J. Donea, A. Huerta, J.-Ph. Ponthot, A. Rodrígues-Ferran, in: Erwin Stein, René de Borst, Thomas J.R. Hughes (Eds.), Encyclopedia of Computational Mechanics, John Wiley, 2004. [6] K.E. Evans, M.A. Nkansah, I.J. Hutchinson, S.C. Rogers, Nature 353 (1991) 124. [7] R. Lakes, Science 235 (1987) 1038. [8] L.D. Landau, E.M. Lifshits, Theory of Elasticity, Pergamon, London, 1986. [9] R.B. Hetnarski, J. Ignaczak, Mathematical Theory of Elasticity, Taylor & Francis, New York, London, 2004. [10] J.F. Semblat, J. Sound Vib. 206 (5) (1997) 741. [11] Ch. Smith, K.W. Wojciechowski, Phys. Stat. Solidi (b) 245 (2008) 486. [12] T. Strek, B. Maruszewski, J. Narojczyk, K.W. Wojciechowski, J. Non-Cryst. Solids 354 (2008) 4475. [13] T.C.T. Ting, T.Y. Chen, Q. J. Mech. Appl. Math. 58 (2005) 73. [14] T. Trombetti, S. Silvestri, J. Sound Vib. 292 (2006) 21. [15] K.W. Wojciechowski, Phys. Lett. A137 (1989) 60. [16] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, 5th Ed., vols. 1–3, Butterworth-Heinemann, Oxford, 2000.