Bubble Deformation and Growth Inside Viscoelastic Filaments

Article pubs.acs.org/IECR

Bubble Deformation and Growth Inside Viscoelastic Filaments Undergoing Very Large Extensions J. Papaioannou, A. Giannousakis, Y. Dimakopoulos, and J. Tsamopoulos* Department of Chemical Engineering, Laboratory of Computational Fluid Dynamics, University of Patras, Patras 26500, Greece ABSTRACT: Motivated by the probe experiment for characterizing the adhesion strength of polymeric materials, we studied the axisymmetric extensional flow of a viscoelastic liquid filament that contains one or three initially spherical bubbles along its axis of symmetry. The filament is confined between two parallel and coaxial disks of the same radius and its initially cylindrical outer surface is surrounded by air. The flow is induced by the axial extension of the upper disk under constant velocity, while the lower disk remains stationary. The rheology of the liquid is described by the exponential Phan−Thien and Tanner viscoelastic model. A quasi-elliptic transformation combined with domain decomposition and local mesh refinement is employed to discretize the domain. The mixed finite element Galerkin method for the mass and momentum balances combined with the EVSS-G method and SUPG weighting for the constitutive equation are used to obtain their weak forms. All these are necessary for the successful simulation of much larger filament extensions and of liquids with higher elasticity than those reported earlier. The evolution of the filament and the bubble(s)-free surfaces depends on the interplay of the viscous, elastic, and capillary forces. It was found that, when the ratio of the elastic and the capillary forces is small compared to the viscous forces, the bubble attains very large deformations. The onset of cusps was observed at the bubble poles for intermediate values of capillarity and elasticity. The force varies with the filament extension (or time) in a way reported in experiments measuring the strength of adhesive materials: Its initial increase up to a maximum is followed by a plateau, and finally it drops to zero when the adhesive tends to fail. Increasing the filament elasticity delays the development of stresses (and the applied force).

1. INTRODUCTION Pressure-sensitive adhesives (PSAs) are the most common type of adhesives found in various consumer products. Today’s PSAs are based on block copolymers (polystyrene, polyisoprene, or polyacrylates) and natural rubber. Typically, they must have an elastic modulus up to a certain threshold value, high viscosity to dissipate energy upon debonding and undergo strain-hardening to avoid forming residue on the solid surface. They are designed to stick to almost any surface at room temperature by applying only light pressure, without any chemical reaction, solvent evaporation, or UV irradiation.1 To verify the extent to which a particular PSA has these characteristics, its relevant mechanical properties must be carefully measured via specifically chosen methods. To this end, the peel and the probe tests are used.2 According to the latter method, a cylindrical flat-ended probe comes into contact with a thin layer of the adhesive which is placed on another flat substrate. After a small compression on the PSA is applied, the probe is removed axially at a constant velocity, while the required force to do so is measured. This force initially increases almost linearly with time, and subsequently after attaining a maximum, it reaches a distinctive plateau and finally approaches zero as the PSA film breaks apart. Commonly, the maximum force corresponds to the appearance of tiny bubbles either in the bulk of the filament due to cavitation during material decompression or on its contacts with the two solid surfaces due to gas entrapped in even minute crevices of the substrates.3 The plateau lasts longer in those PSAs that are able to sustain high peel forces and is associated with bubble growth which turns the bulk of the PSA into several fibrils (fibrillation process). This assists the polymer to average the stresses in its volume, thus avoiding rapid debonding and material failure, which is observed in structural © 2014 American Chemical Society

glassy adhesives. On the contrary, lateral bubble growth and coalescence must be avoided, especially on the solid surface, because it leads to detrimental de-adhesion through crack propagation. Cavity growth and coalescence is decelerated by increasing the elasticity and strain hardening of the adhesive.4−6 Clearly, understanding bubble growth and coalescence in such materials is essential to understanding their mechanical properties.7 Less frequent, is the formation of fingers on the three-phase contact line between adhesive, substrate, and surrounding air, which grow following a mechanism similar to the one leading to the Taylor−Saffman type instability.8 Thus the study of the dynamics of such phenomena can illuminate further the debonding mechanisms of PSAs and is very important for designing the composition of high quality adhesives. Experimental state-of-the-art work in the area focuses on the development of methods for precisely controlling the roughness of both the probe and the adhesive and on systems with a well characterized rheology in order to link stickiness with molecular architecture.9 Apparently, understanding and simulating the deformation of a filament containing bubbles and undergoing stretching is a very demanding task. On the contrary, the decoupled problems of a single bubble growth in a viscoelastic material or of a filament undergoing axial deformations have been studied quite extensively. For example, Pearson and Middleman10 have Special Issue: John Congalidis Memorial Received: Revised: Accepted: Published: 7548

October 3, 2013 March 10, 2014 March 16, 2014 March 16, 2014 dx.doi.org/10.1021/ie403311n | Ind. Eng. Chem. Res. 2014, 53, 7548−7569

Industrial & Engineering Chemistry Research

Article

and interaction with each other all situated on the solid substrate or in the bulk of the stretched complex liquid and how they affect the mechanical or tacky properties of the adhesive, have only recently been addressed.31 It is the purpose of this and forthcoming reports to resolve some of these questions through a systematic investigation that relies on stateof-the-art numerical algorithms applied on geometries of relevance to the available probe experiments and to use relevant constitutive models. The present study focuses on the calculation of bubbles growing along the axis of symmetry of a stretched viscoelastic filament. Our goal is to report simulation results achieving (a) large extensions of filaments containing bubbles, (b) work with realistic constitutive models and values of the Deborah number, and (c) accurately compute the deforming interfaces and resolve the stress boundary layers arising there. To this end, we will use two different discretization methods: in the first one, the entire volume of the liquid in the filament will occupy a single domain which will be discretized using the quasi-elliptic grid generation technique. This method has been used successfully already in problems with various complex and evolving geometries.29,32 In the second discretization method, the fluid domain will be decomposed into multiple domains, each one surrounding each bubble and the last one covering the remaining liquid.33 In principle, the latter method should be able to more accurately resolve the stress boundary layers formed not only on the solid−liquid interfaces, but also on the bubble-liquid interfaces, because it controls more closely the structured grid in each domain. This method could also be easier to adopt in three-dimensional flows, which are of interest in the probe test. In sections 2 and 3, we present the problem to be solved and the method of solution, respectively. In section 4, we discuss the results of our simulations and present a complete parametric analysis of the uniaxial extension of viscoelastic filaments with different initial aspect ratios and bubble positions. Finally in section 5, we present our conclusions and describe future directions.

examined the collapse of a bubble in a Maxwell viscoelastic medium and concluded that elasticity accelerates it. Ting11 has addressed how the initial bubble growth rate is affected when polymer is added to a Newtonian solvent. Papanastasiou et al.12 have studied the bubble collapse in a viscoelastic medium modeling its elastic nature with the integral constitutive equation K-BKZ, while neglecting the inertial terms in the momentum equations. Finally, Bousfield et al.13 have investigated the transient deformation of a gas bubble in a stretched viscoelastic medium and concluded that large compressive stresses are responsible for the faster collapse of the viscoelastic liquid. The dynamics of liquid bridges, that is, filaments confined between parallel and coaxial disks, when one of the disks undergoes oscillations has been studied theoretically by Chen et al.,14 Tsamopoulos et al.,15 Chen and Tsamopoulos,16 and experimentally by Mollot, et al.17 to determine viscosity and surface tension of the liquid. Filament stretching rheometers have been proposed to determine the extensional viscosity of polymers.18 However, their contact with the upper and lower disks and the no-slip boundary condition applied there inevitably introduce shear in the extending filament, except for its midsection, and spatially and temporarily varying extension because of the nonuniform decrease of the filament radius along the axial direction. To reduce them, one could use longer filaments, but this would increase the importance of gravitational sagging and, hence, complicate the measurements. All these necessitate adjustments in the pulling velocity of the upper disk, which, in this test, varies exponentially with time and using filaments with initially small aspect ratio, Λ, defined as the distance between the disks over their radius. Additionally careful analysis is required to determine the true extensional viscosity and other characteristics of the material as well as the shape of the filament depending on the viscoelastic model used.19−21 Sometimes the filament develops nonaxisymmetric shapes and fibrillation at its contact with the solid disks.22−24 Simulations of the probe experiment which is relevant for PSA characterization, that is, filament stretching containing bubbles, are quite limited. For example, Foteinopoulou et al.25 have investigated the growth of a single bubble in Newtonian and viscoelastic filaments subjected to stretching at constant velocity, as in the probe experiment. They have used an algebraic mapping to a fixed domain (for this technique see, for example, Poslinski and Tsamopoulos26). In this way, they have computed accurately the deforming interfaces and have presented results for a variety of capillary numbers. Because of this mapping however, they were limited to liquids of low elasticity, filament aspect ratios, Λ ≥ 1, and reached rather low filament extensions (typically up to 200%). They also studied the onset of cavitation with criteria proposed by Joseph27 and found that elasticity delays its appearance. Foteinopoulou et al.28 have also studied the growth of one or more bubbles, but in a Newtonian filament, using the quasi-elliptic meshgeneration technique advanced by Dimakopoulos and Tsamopoulos.29 In this way, they obtained converged solutions even at very small aspect ratios and large filament extensions. More recently, Sujatha et al.30 have simulated a single bubble growing in the middle of a Newtonian or highly viscous elastic filament, which undergoes stretching because of an exponentially increasing axial velocity of the upper plate. They have compared results obtained by either an entirely arbitrary Lagrangian-Eulerian (ALE) or a so-called ALE/VOF method. Other more complicated issues, such as multibubble growth

2. PROBLEM FORMULATION We consider the isothermal and incompressible extensional flow of a viscoelastic liquid that exhibits a relaxation time λ* and has a total dynamic viscosity μ* = μp* + μs*, where μp* is the viscosity of the polymer and μs* is the viscosity of the Newtonian solvent. We assume that the fluid is incompressible with a constant density ρ* and the fluid-air interfaces have a constant interfacial tension σ*. Symbols followed by an asterisk denote dimensional quantities. Similar to our previous approaches25,28 for stretching filaments containing a single or multiple bubbles, the liquid filament is considered to have initially a perfect cylindrical outer surface of radius Rc*, and be confined between two solid and coaxial disks of radius also Rc*, at an initial distance H*0 from each other. Extending the liquid up to the radial edges of the confining disks helps pinning the contact line there. Inside the volume of the liquid one or more gaseous bubbles exist, which are initially spherical, with radius R*b,i, i = 1, 2, ... and at a distance h*0,i above the lower disk. To retain the assumption of axial symmetry, the center of each spherical bubble must lie along the axis of symmetry of the filament. A schematic of the system considered with an arrangement of three bubbles is shown in Figure 1. The filament undergoes deformation by pulling the upper disk with a constant velocity U*, while the lower disk, onto which the filament adheres permanently, remains fixed. The origin of the 7549

dx.doi.org/10.1021/ie403311n | Ind. Eng. Chem. Res. 2014, 53, 7548−7569


Article

In the corresponding probe experiment, the size of the liquid column is small enough and the viscous and elastic forces are large enough so that gravitational body forces and inertia are negligible. Hence, Re = 0, Bo/Ca = 0, although keeping either one does not increase the computational difficulties. This simplification makes the filament symmetric about its instantaneous midplane, if the bubbles arise symmetrically in it. A large number of models have been proposed over the years to describe the rheology of PSAs.34,35,31 Quite recently however, it has been found that the two-mode PTT36 model of the exponential type describes the rheology of a large class of such materials quite satisfactorily.37 Moreover, Yao et al.38 have shown that one or two additional modes of the Giesekus model, which like the PTT model is weakly strain hardening, improve only quantitatively the predictions of the single mode model in filament stretching simulations, and primarily during the very early and the late stages of stretching. Hence, for simplicity, the single mode PTT model will be used in the present study. Its dimensionless form is ∇

Y (τ p)τ p + De τ p = 2(1 − β)γ ̇

where γ ̇ = (∇̲ u̲ + ∇̲ u̲ T)/2 is the rate of strain tensor, while the symbol ∇ over the viscoelastic stress represents the upperconvected Maxwell derivative defined as

Figure 1. Schematic of the filament geometry containing three bubbles.

cylindrical coordinate system is placed at the center of the lower disk. Consequently the height of the filament increases linearly with time, H*(t*) = H*0 + U*t*, and the bubbles following the dynamics of the filament deform, expand, and translate along the filament axis. The presence of these multiple deforming interfaces and their coupling with the nonlinear flow field and constitutive equation render the solution of this problem a very difficult task. We scale all lengths with the radius of the lower bubble, R*b,1 and velocities with the velocity U* of the upper disk. In addition, time is scaled with R*b,1/U*, while the viscous scaling μ*U*/R*b,1 is used for the pressure and the stresses. Thus the dimensionless parameters that arise are the Reynolds number, Re = ρ*U*Rb,1 * /μ*, the Capillary number, Ca = μ*U*/σ*, the Bond number, ρ*g*R*b,1/σ, the Deborah number, De = λ*U*/ R*b,1, and the ratio of the solvent viscosity over the total viscosity, β = μs*/μ*. The initial filament configuration is set by the following geometric ratios Rc = Rc*/Rb,1 * , H0 = H0*/Rb,1 * , Rb,i = R*b,i/R*b,1, h0,i = h*0,i/R*b,1, and Λ = H*0 /R*c , which is the initial filament aspect ratio. In the case of a single bubble we use h in place of h0,1. The system of equations governing the motion of a viscoelastic fluid consists of the mass and momentum balances. Following the aforementioned scaling the resulting dimensionless mass and momentum equations are (1) ∇̲ · u̲ = 0 Re

D u̲ Bo = −∇̲ P + ∇̲ ·τ + e̲ z Dt Ca

∇

τp =

Dτ p Dt

− τ p·∇̲ u̲ − ∇̲ u̲ T·τ p

(5)

The function Y(τ p) distinguishes between the two PTT fluid models. Here the exponential PTT model is used, according to which Y(τ p) = exp[εDe/(1 − β)tr(τ p)]. This viscoelastic constitutive model was derived under the assumption of a transient network created by temporary cross-links of the polymeric segments. The relaxation time of the polymer is related to the inverse probability that a strand will be destroyed. Since a single relaxation time is used in this study, all strands are statistically the same. Besides the liquid viscosity μ* and the relaxation time λ*, two more parameters control the model predictions. The first one, ξ, is related to the nonaffine motion of the polymeric chains and intermolecular slip. Here ξ is set to zero and the model is the so-called affine PTT model. The second one is the dimensionless parameter ε, which introduces an upper limit in the elongational viscosity. When this parameter increases, both shear and elongational viscosity decrease. It usually takes values quite smaller than 1. Clearly the PTT model reduces to the OLDROYD-B model by setting ε = 0 and to the UCM model by setting β = 0 as well, both of which predict constant shear viscosity and extensional viscosity that increases without bound for De ≥ 0.5. The governing equations are discretized using the finite element method. To solve accurately and efficiently various viscoelastic flows at even high Deborah numbers,39 we introduced the elastic-viscous split stress (EVSS) formulation. This method consists of splitting the polymeric part of the extra stress tensor into a purely elastic and a viscous part

(2)

where u̲ and P stand for the axisymmetric velocity vector and the isotropic pressure respectively, while D/Dt = ∂/∂t + u̲ ·∇̲ is the material derivative and ∇̲ is the gradient operator. Moreover, τ is the extra stress tensor, which is split into a purely viscous part due to the Newtonian solvent 2βγ̇ and a polymeric contribution τ p. τ = τ s + τ p = 2βγ ̇ + τ p

(4)

τ p = Σ + 2(1 − β)γ ̇

(6)

The success of this scheme resides on the fact that the elliptic nature of the momentum equations is ensured even in the absence of solvent, β = 0. Brown et al.40 proposed a modification of this formulation (EVSS-G) according to which an independent interpolation of the components of the

(3) 7550



Article

Pi ,idealVi γ = const,

velocity gradient tensor is introduced in order to satisfy the compatibility in the approximation between elastic stress and velocity gradients in the constitutive equation. The corresponding (additional) equation that must be solved is

G = ∇̲ v̲

Thus after reformulating the momentum and constitutive equations using the EVSS-G idea, we obtain

DF̲ = u̲ , Dt

∇

Y (τ p)Σ + De Σ + 2De(1 − β)D − 2(1 − β)(1 − Y (τ p))D = 0̲

t ̲ · u̲ = 0, (9)

where D = (G + G )/2 . In addition the definition of the upper-convected Maxwell derivative becomes τp =

Dτ p Dt

− τ p·G − GT ·τ p

2Hn̲ Ca

n̲ d t̲ − R2 ds

(19a, b)

The model is completed assuming that initially the liquid and the gas in the bubble(s) coexist under the conditions of hydrostatic and thermodynamic equilibrium. For the velocity field u̲(r,z,t = 0) = 0̲, the stresses τ = 0, while for the initial pressure in the liquid and the gas, we set the ambient pressure to zero (datum pressure) Pamb = 0, and then we have,

(11)

initial pressure in the filament: P(t = 0) = Pamb +

1 R cCa

(20)

initial pressure inside each bubble: PG, i(t = 0) = Pamb +

1 2 + R cCa R b, iCa

(21)

3. NUMERICAL IMPLEMENTATION To numerically solve the governing equations and accurately simulate the growth of one or more bubbles confined in a stretching filament following its large deformations, we chose the mixed finite element method combined with an elliptic grid generation scheme for the discretization of the highly deforming physical domain. This method was combined with local mesh refinement to resolve the flow with a structured grid, where it was necessary. Moreover, due to the continuously deforming domain primarily in one direction, a remeshing procedure was developed in order to add/subtract elements to/ from specific regions of it. 3.1. Elliptic Grid Generation. We use a set of partial differential equations the solution of which generates a boundary fitted discretization in the space that is occupied by the liquid, which undergoes large anisotropic deformations. This method was advanced by Dimakopoulos and Tsamopoulos29 and was successfully applied in various complex flows.28,32,41,42 Here we will only present our adaptation of its essential features to the current problem and combine it with a domain decomposition method to generate improved discretizations. Further details on the important issues of the method are available in ref 29. With this scheme the time dependent (physical) domain is mapped onto a fixed in time computational one, through the one-to-one transformation:

(13)

To avoid dealing with the second order derivatives that arise in the form of the mean curvature, we split it into two parts 2Hn̲ =

(18a, b)

n̲ · u̲ = 0, t ̲ n̲ · σ = 0.

I is the unit tensor, PG is the gas pressure, and 2H is twice the mean Gaussian curvature of the corresponding interface defined as ∇̲ s = (I − n̲ n̲ ) ·∇̲

(17a, b)

On the axis of symmetry, flow symmetry holds, (10)

where the left-hand side (LHS) is the term that arises in the line integral upon applying the divergence theorem on the residual of the momentum equation and is substituted by the right-hand side (RHS) of eq 11. Moreover, n̲ is the unit normal that is always directed outward from the bubble or the outer interface of the filament and is defined as zs er − rs ez n̲ = rs2 + zs2 (12)

2H = −∇̲ s · n ,

n̲ · u̲ = 0

t ̲ · u̲ = 0, n̲ · u̲ = 1

In terms of the boundary conditions, along the free surfaces the velocity field should satisfy a normal force balance between surface tension, viscous stresses, and pressure in the liquid and pressure in the inviscid gas − n̲ ·( −PI + τ ) = PG n̲ ±

(16)

On the upper disk, the no slip and the constant axial velocity are used,

T

∇

F̲ = r e̲ r + z e̲ z

where F̲ is the position vector. On the surface of the lower disk, the usual no slip and no penetration conditions are imposed,

(8)

∇

(15)

where Vi denotes the instantaneous volume of each bubble, while γ = 1.4. All liquid/gas interfaces move according to the kinematic condition,

(7)

∇̲ P − ∇̲ ·Σ − 2∇̲ ·γ ̇ = 0̲

i = 1, 2, ...

(14)

where the first term describes the change of the tangential vector along the free surface, while R2 is the second principal radius of curvature of the interface (in the azimuthal direction). We denote here the unit normal that points always from the liquid to the gas with n .̂ Each one of the two types of interfaces (filament/surrounding gas or bubbles/filament) should be treated differently when applying eq 11. More specifically, (a) at the filament/ambient gas interface, PG = Pamb = 0, n = n ̂ and the minus (−) sign is used in front of the capillary term. (b) at the bubble/filament interface, PG,i = Pi,ideal, n = −n ̂ and the plus (+) sign is used instead. The instantaneous gas pressure Pi,ideal of each bubble is calculated by the adiabatic gas law: 7551


Industrial & Engineering Chemistry Research J

(r , z ) → (ξ , η)

Article

two methods in terms of their accuracy, computational cost, and difficulty in implementing them. The mesh consists of two segments (see Figure 2a), one around the bubble and a second one covering the rest of the control volume. The basic elements of this method include the following: (a) The first segment (no. 1 in Figure 2a) is created by using polar spherical coordinates, (ρ ∈ [ρ1 = 1, ρ2], θ ∈ [0, π]),

(22)

The mesh is generated first in the fixed domain, which is chosen to be a perfect, fixed, and completely filled with liquid cylinder. The inverse mapping generates the “deformable” mesh, which closely follows the deformations of the liquid volume. The mapping is based on the solution of the following system of partial differential equations, ⎛ ⎞ rη2 + z η2 ⎟ (1 ) + − ∇̲ ·⎜ε1 2 ε 1 ⎟∇̲ η = 0 2 ⎜ r z + ξ ξ ⎝ ⎠

(23)

∇̲ ·∇̲ ξ = 0

(24)

where ε1 is a parameter that adjusts the orthogonality of the resulting mesh; its value is chosen by trial and error and here it is set to 0.1. The mesh in the “deformable” domain changes, with a velocity that is not necessarily equal to the local fluid velocity, and therefore this method belongs to the group of ALE (Arbitrary-Langrangian-Eulerian) methods. The initial physical domain is the filament with a perfect cylindrical outer surface containing the bubble or the bubbles of a perfect spherical shape. To construct the mesh in physical space starting from that in the computational domain and following the single domain method, the equation of a prolate ellipse is imposed where the bubble(s) are located. The boundary of the ellipse is mapped onto a line segment of the same length with the ellipse major axis and at the same position on the axis of symmetry in the computational domain. Subsequently, the ellipse is inflated continuously in the direction of its minor axis only, until it becomes a sphere. The location of this boundary and all other boundaries, together with a condition that uniformly distributes the interfacial nodes are used as boundary conditions for eqs 23 and 24; hence, the intermediate meshes are generated in physical space and finally the mesh to initiate the dynamic simulations, as described in Foteinopoulou et al.28 Having generated the initial mesh, the flow equations are solved in the entire domain. Then, the kinematic condition, eq 16, is used to advance the liquid/gas boundaries of the physical domain in time. Especially for the nodes at the poles of the bubble(s), the kinematic condition is reduced to the following essential conditions for the mesh in the z direction ∂Zil = uz|ξ = 0, η = Zi1,0 , ∂t

∂Ziu = uz|ξ = 0, η = Ziu,0 ∂t

Figure 2. Schematic of the mesh in the domain decomposition method with mesh refinement at the free surfaces for a viscoelastic liquid with De = 1, Ca = 10, β = 0.01, ε = 0.03, and Λ = 1 that contains one gas bubble placed initially in the midplane of the filament. (a) The two mesh segments that will be joined to form the global mesh at the beginning of the simulations (at t = 0, only half of the domain is shown) and (b) the two mesh segments and the global mesh at t = 9.88.

where typically ρ2 = 1.5, and then changing the coordinate system into cylindrical (r, z) with the usual transformations. Local mesh refinement is also seen next to the bubble interface. This will be described in section 3.2. (b) The second segment (no. 2 in Figure 2a) is generated similarly with the entire domain in the single domain method, starting with an outer cylindrical surface and an inner surface, which follows the equation of the prolate ellipse. The latter surface is again inflated to obtain a spherical shape of radius, ρ2, that is, the outer radius of segment no. 1. (c) Since the two segments have been constructed in the same coordinate system and they have the same spherical interface, they are connected to form the initial physical domain. During the dynamic simulations and while the mesh is regenerated, the common boundary of the two segments is a priori unknown. It is defined as a curve drawn by expanding the bubble surface from the previous time step tn along the normal vector n̲. The coordinates of the desired boundary curve (bc) are thus equal to

(25)

where i is an index that refers to the ith bubble and the superscripts l and u denote the lower or the upper pole, respectively. The boundary location of the lower stationary disk (z = 0) remains fixed, while the boundary location of the upper disk is advanced following its constant velocity, and in dimensionless form it is given by z = H(t) = Ho + t. All new boundary locations of the physical space together with a node equidistribution condition are again used as boundary conditions when solving eqs 23 and 24 to generate the mesh in the new physical space. This entire process is repeated in every time step. The modification of this procedure when using the domain decomposition method is described below. 3.2. Domain Decomposition Method. This method was developed for a filament containing a single bubble in order to more accurately solve for the stress and velocity fields near the deforming bubble. It also gives us the capability to compare the 7552



Article

rbc(t n + 1) = rbs(t n) + d ·( n̲ ·er)

(26a)

z bc(t n + 1) = z bs(t n) + d ·( n̲ ·ez)

(26b)

3.4. Local Mesh Refinement. To preserve the high accuracy of computations and resolve the flow in regions of particular interest, we employed the h-refinement method. The h-method was proposed by Szabo and Babuska43 who subdivided the elements in which the error measure was larger than a prescribed tolerance. Later, Tsiveriotis and Brown44 applied a nonconforming splitting of rectangular elements in a free boundary problem and argued that local mesh refinement is essential in cases where elliptic grid generators are used, because this technique relaxes the requirements on the mapping equations and provides a greater flexibility on the handling of the grid. Quite recently, Chatzidai et al.33 have presented a detailed analysis and results of an alternative conforming splitting of triangular elements to refine their finite element meshes. The benefits of this method, as explained there, are 2-fold: the considerable reduction of the overall computational cost and the increase in the accuracy of the calculations. In the current study, the stretching filament dynamics cause in some cases elongation of the filament and the bubble interface several times their initial arc length. Moreover strong boundary layers appear in the vicinity of these regions, which need fine discretizations in order to be resolved. Thus, we employed the method discussed in Chatzidai et al.,33 as follows: (a) A uniform mesh of rectangular elements is initially created. (b) The region of interest is selected and is refined by splitting each element into four equal elements. (c) A transition zone of triangular elements is introduced in between in order to make the elements of the two different regions conforming. (d) All remaining rectangular elements are subdivided into two triangular ones, because triangles follow closer the large distortions of the deformable mesh. The local mesh refinement technique was developed for both of the domain methods employed herein. However it can be shown that one can fully take the advantage of this method if the domain decomposition method along with remeshing is used. Figure 3 shows a refined mesh in the case of the single domain method for the same viscoelastic liquid as in Figure 2b and a close time instant t = 9.99. In this case, it can be seen that a large number of refinement elements is placed on the axis of

where (rbs,zbs) are the coordinates of the bubble free surface (bs) from the previous time step and d = ρ2 − ρ1 is a constant that defines the distance between the two curves. The previously defined coordinates are then interpolated to optimize the distribution of the nodes along the boundary. (d) The flow equations are mapped directly from the unit triangular element onto each triangle of the physical domain without a transformation to the computational mesh. A high quality mesh that is generated by this method is illustrated in Figure 2b. Here we show the two separate segments, no. 1 and no. 2, and their union, no. 3, form the entire extended filament. This mesh discretizes the volume of a viscoelastic liquid with De = 1, Ca = 10, β = 0.01, ε = 0.03, and Λ = 1 at t = 9.88, that contains one gas bubble placed initially in the midplane between the disks. 3.3. Remeshing Procedure. Despite the high quality of the initial mesh, the large deformation of the filament, primarily in the axial direction, results in large deformation of the mesh in physical space. To deal with the large aspect ratios of the elements, a new remeshing procedure was developed along with the domain decomposition method during which elements are added to and/or subtracted from the “deformable” local mesh domains. Time integration and all computations are halted, and complete mesh reconstruction is enforced, when the filament deformation exceeds a certain level. To determine the number of elements to be added, the following criterion is examined: || Nz || > k max(dle(trm))

(27)

where ∥·∥ stands for the Euclidean norm of Nz = dle(t ) − dle(trm), which is the difference between the length of the eth element at the current time step (tn+1) in the z-direction and the corresponding one that the same element had initially or just after the previous remeshing (t = trm), while k is a factor that determines the number of elements to be added. This is a nice way to fill locally the additional space in the axial direction, due to the extension of the filament. However, the regions above or below the bubble may be compressed owing to the faster axial deformation of the bubble in comparison to the motion of the upper disk. If this is the case, the criterion still holds, but now k determines the elements to be subtracted from these regions. The remeshing procedure is summarized in the following steps: (a) As the initial shape of the segment that surrounds the bubble we choose a hollow spherical shell, while for the other segment we choose a cylinder of height equal to that of the physical domain at the current time instant. (b) Initially, the segment surrounding the bubble is deformed until its inner surface coincides with the current shape of the bubble. Its outer surface is determined by expanding the bubble surface along its normal vector. The node positions of the outer surface are used as essential conditions for the inner boundary of the second segment. (c) All the interfaces attain the shape that they had just before the remeshing using quadratic interpolation with nodes that are equidistributed along the boundary curve. In this way, the physical domain is created using the new node positions of both segments. (d) The flow variables (velocities, pressure, stresses, and velocity gradients) are finally calculated on the new nodes of the physical domain from the old ones using search and finite element interpolation techniques. n+1

Figure 3. Schematic representation of the global mesh in the single domain method for the case of the viscoelastic liquid under the conditions given in Figure 2 at t = 9.99. 7553



Article

symmetry, a region of no special interest. This happens because the axis of symmetry along the bubble interface is mapped onto the same ξ = constant curve. Therefore, if we place refined elements on the bubble interface, the same refinement is required for the entire axis of the domain, a fact which increases the total computational time and memory. On the other hand, in the domain decomposition method, Figure 2b, that excess of nodes is easily avoided by construction, because the mesh around the bubble can be isolated and treated separately. In this case the bubble and the axis of symmetry belong to the ξ = constant curve of two different domains. 3.5. Mixed Finite Element Method. The discretization of the flow and mesh equations is performed using the mixed finite element/Galerkin method. Triangular elements are used, as already described. We approximate the velocity and position vectors with 6-node Langrangian basis functions, φj, and the pressure with 3-node Langrangian basis functions, ψj, in accord with the LLB43 condition. Different approximations have been used in the past for the stresses and the velocity gradient tensor. Szady et al.45 used bilinear basis functions for the velocity gradient, while biquadratic ones were used for the velocity. The stresses according to them were approximated with bilinear basis functionsthe same functions as the velocity gradient tensorto ensure element compatibility when u̲ = 0. Here, we use 3-node linear basis functions for both tensors, as for the pressure. Applying the divergence theorem on the weak form of the momentum balances we obtain

∫Ω ψj(G − ∇̲ u̲) dΩ = 0

Finally, the constitutive equation that is a hyperbolic one is solved by the Streamline Upwind Petrov-Galerkin (SUPG) method46 ∇

∫Γ [ n̲ ·(−PI + Σ + 2γ̇)]φj d Γ = 0

− 2(1 − β)(1 − Y (τ p))D]Ψj dΩ = 0

Ψj = ψj +

⎛ ⎝

+L

rη2 + z η2 rξ2 + zξ2

∫Γ

∂φj ∂η

(28)

(29)

∫Ω ∇ξ·∇φj dΩ = 0

(34)

4. RESULTS AND DISCUSSION We performed numerical simulations in order to study the transient stretching of viscoelastic filaments containing one or three bubbles. The flow is induced by the pulling of the upper plate with a constant velocity. The motion of this plate results in the sagging of the filament free surface toward the axis of symmetry and the deformation of the bubble(s). The inertial terms in the Navier−Stokes equations have been neglected in all simulations, since the materials considered here are known to have high viscosity values (up to 107 Pa s), resulting in practically zero Reynolds numbers. Thus the characteristics of the flow, that is, the velocity, pressure, and stress fields, as well as the free surface shapes and the force on the upper plate are quantities that depend upon the interplay of capillary, viscous, and elastic forces. Hence we present a complete parametric analysis to describe the characteristics of such flows, with

⎞ + (1 − ε1)⎟∇η ·∇φj dΩ ⎟ ⎠

rη2 + z η2 dη = 0

δt n + 1 u̲ ·∇̲ ψj 2

where δtn+1 is the current time step. 3.6. Solution Procedure. The resulting set of discrete equations is integrated in time with the implicit Euler method. An automatically adjusted time step is used for that purpose which ensures the convergence and optimizes code performance (see Dimakopoulos and Tsamopoulos).41 The nonlinear set of algebraic equations is solved iteratively in each time step using a Newton−Raphson/Nonlinear Gauss-Seidel iteration scheme, which decouples the solution of the previously defined set of equations and splits the procedure into four subproblems. The momentum equations are solved on the nodes of the physical domain, determined from the previous time step until convergence, and only the appropriate information from the computed flow variables is subsequently used for the solution of the spatial equations. For the domain decomposition method the two different mesh segments evolve in time separately, and then they are connected to form the physical domain. It should be noted here that the constitutive equation is solved separately. The iterations of the Newton−Raphson method are terminated using 5 × 10−10 as tolerance for the absolute error of the residual vector. This is an effective method because it results in significantly smaller matrices, which are easier to handle. The Jacobian matrices that are generated are stored in Compressed Sparse Row (CSR) format, and the linearized system is solved by Gaussian elimination using PARDISO, a robust, hybrid matrix solver.47,48 The initial time step for all the simulations was δt0 = 10−5. The codes were written on Fortran 90 and were run on a workstation with 2 Quad Xeon CPU at 2.4 GHz in the Laboratory of Computational Fluid Dynamics.

where dΩ and dΓ are the differential volume and surface area, respectively. The boundary term in the momentum equation is split into as many parts as the boundaries of the physical domain and the relevant boundary condition is applied on each one. The term of the moving interface integral involving the stress tensor is substituted with the mean curvature of the interface and is treated as explained in the previous section. At the remaining boundaries the momentum balances are replaced by the essential conditions imposed therein. The weak form of the mesh generation equations are derived similarly by applying the divergence theorem,

∫Ω ⎜⎜ε1

(33)

where the function Ψi is formed from the finite element basis function for the elastic stress tensor according to

whereas the weak form of the mass balance is

∫Ω ψj∇̲ · u̲ dΩ = 0

∇

∫Ω [Y (τ p)Σ + DeΣ + 2De(1 − β)D

∫Ω [∇̲ φj·(−PI + Σ + 2γ̇)] dΩ −

(32)

(30)

(31)

where the penalty parameter is L = O(103 −105) and the line integral is along the free surfaces. Moreover the weak form of eq 7 for the velocity gradient tensor is given by 7554



Article

in each simulation is not constant but varies as explained previously, the parameter which will be examined is the maximum allowed time step, which is also the time step during most of the simulation. Results for the case of three bubbles and for four different maximum allowed time steps are shown in Figure 5. Here the common mesh that we used is M1. The

respect to the Deborah and the capillary numbers and the geometric ratios characterizing the filament. 4.1. Code Validation. First, we verified that the solution with the single-domain mesh converged with mesh refinement and time step choice and that the new code could reproduce the earlier results.25,28 Additionally, we checked the performance of adapting the time step, given the variation of the solution between time steps. Typical results of these tests for the case of three bubbles are shown in Figure 4. The solution

Figure 5. Convergence test for a filament with three bubbles. Evolution with time of τrr,p, P, τrz,p, and τθθ,p at the lower pole of the upper bubble for various maximum allowed time steps. The dimensionless parameters are those in Figure 4.

Figure 4. Convergence test for a filament with three bubbles. Evolution with time of τrr,p, P, τrz,p and τθθ,p at the lower pole of the upper bubble for the three different meshes shown in Table 1. The dimensionless parameters are Re = 0, De = 1, Ca = 10, ε = 0.03, β = 0.01, Λ = 2.5, H0 = 10, Rb1 = Rb2 = Rb3 = 1, h0,1 = 2, h0,2 = 5, h0,3 = 8.

solution is again monitored at the lower pole of the upper bubble. Plainly, all variables are in perfect agreement irrespective of the time step examined. Next, we compared the results obtained either with the single-domain or the multidomain methods. To this end, we considered the case of a single bubble placed initially at the midplane of the filament and we compared the interfaces (of the filament or the bubble) in terms of the Euclidean L2 norm of the relative difference of either the radial or the axial position vector. This norm is defined as

for the pressure and three components of the elastic stress tensor, τrr,p, τrz,p and τθθ,p, over time is shown corresponding to the lower pole of the upper bubble. We chose these variables and this point for our tests because the bubble poles are regions where the stress components may exhibit singular behavior by increasing abruptly over time for certain parameter values, and their accurate calculation is a critical test for our numerical code. The three different spatial discretizations, which are used, are given in Table 1. It should be noted here that the numbers

L 2 = ||(DD − SD)/DD ||

where DD stands for the solution with domain decomposition at every point on the interface, and SD stands for the solution at the same point using the single domain method and, if necessary, interpolation using the arc-length of the surface to get it at the same location as in DD. The results are shown in Tables 2−7. Table 2 presents the case of a viscoelastic liquid filament with (Re = 0, εPTT = 0.03, β = 0.01), De = 1, Ca = 10, and Λ = 1. The three parameters in parentheses are kept constant, except for a brief study on the effect of εPTT and β on our predictions in the end of section 4.2.1.1. The mesh we used and the time instants at which remeshing took place with DD are indicated in Table 3. In this case, the maximum difference in time instants at which the corresponding shapes were calculated is 0.002 dimensionless units. Nevertheless, very good agreement between the two methods is observed, especially at small and intermediate times. Only at the third time instant, when shape deformations are large, does the norm become quite larger. This occurs because the mesh in the single domain method has become coarser at these very large deformations, because remeshing is not used with it. On the other hand, in the domain decomposition method we continuously add elements in the interior of the domains and on the free surfaces thus retaining the accuracy of the solution at high level.

Table 1. Different Meshes Used for a Filament with the Conditions Given in the Caption of Figure 4 with Increasing Number of Nodes in Both Directions meshes

no. of nodes in rdirection

no. of nodes in zdirection

no. of elements in r-direction

no. of elements in z-direction

total number of nodes

M1 M2 M3

21 31 41

201 301 401

10 15 20

100 150 200

4221 9331 16441

(35)

of nodes given in this table are before the local refinement procedure, which is employed subsequently to ensure optimal discretization along the interfaces, that is, additional mesh nodes are added locally as follows: Around the bubble surface and along the axis of symmetry these meshes were refined three times at the radial locations: Rc/4, Rc/8, and Rc/20. Also next to the outer filament surface they were refined once at the radial location 7 Rc/8. It is clear that the three solutions for each variable are in very good agreement, certainly until the variable begins to show singular behavior at large times, which here takes place at t ≈ 12. Subsequently, we checked the solution independency with respect to the time step. Since the time step 7555



Article

Table 2. L2 Norm of the Difference between the Solutions Obtained with the Single Domain (SD) and Domain Decomposition (DD) Methods for the Radial and Axial Coordinates of the Bubble and Outer Filament Interface for Ca = 10, De = 1, Λ = 1 filament

bubble time

r-direction

z-direction

r-direction

z-direction

2.248 ± 0.001 8.250 ± 0.001 15.250 ± 0.001

2.671 × 10−4 2.644 × 10−3 2.187 × 10−2

3.667 × 10−4 3.342 × 10−4 1.405 × 10−4

5.346 × 10−4 3.823 × 10−3 6.844 × 10−3

6.446 × 10−4 1.765 × 10−3 1.835 × 10−3

Table 3. Number of Elements along Either the Bubble or the Filament Outer Surface for Ca = 10, De = 1, Λ = 1a time →

for all times

method →

SD

DD

DD

DD

DD

on the bubble interface on the filament interface

320

240

288

472

624

160

t = 10−5 t = 2.248 t = 8.250

320

344

544

Table 5. Number of Elements along Either the Bubble or the Filament Outer Surface for the Conditions Ca = 100, De = 100, Λ = 1a

t = 15.250

for all times

time →

712

a

The time instants are indicated at which the elements with DD increase due to remeshing. The initial number of nodes and triangular elements in the entire liquid volume are (for SD) 37349 and 18240, respectively, and (for DD) 42941 and 21100, respectively.

t = 10−5 t = 3.2485

t = 5.4985 t = 7.7485

method →

SD

DD

DD

DD

DD


320

240

344

400

448

160

320

424

488

536

a

The time instants are indicated at which the elements with DD increase due to remeshing. The initial number of nodes and triangular elements in the entire liquid volume is (for SD) 37349 and 18240, respectively, and (for DD) 42941 and 21100, respectively.

Table 4 compares the two methods for a viscoelastic liquid with De = 100, Ca = 100, and Λ = 1, and the mesh details are given in Table 5. Here the maximum difference between time instances is smaller, 0.001, but the norms retain very satisfactory levels, especially for the filament interface whose curvature does not present steep variations. The results of the last case we examined are given in Table 6. Here the aspect ratio is 5 times smaller, that is, Λ = 0.2 with De = 1 and Ca = 1, and the mesh details are given in Table 7. Again the calculated interfaces are clearly very close. 4.2. Single Bubble at the Midplane. 4.2.1. Effect of Physical Parameters on Flow Field and Bubble and Filament Deformations. Turning into our new results, we note that this work extends the results of Foteinopoulou et al.,25 which were limited to relatively low Deborah numbers (only up to 5) and small filament deformations, whereas commercial PSAs are known to be characterized by high relaxation times, resulting in Deborah numbers well above 10. Here we achieved with our new robust codes converged solutions for much higher filament deformations for the same Deborah numbers, as well as for much higher Deborah numbers, up to 100. All subsequent results were obtained using the single domain method. 4.2.1.1. Large Aspect Ratio, Λ = 1. First we will discuss cases in which the disk radius and height are 4 times larger than the bubble radius. All calculations with one bubble and Λ = 1 were conducted with a mesh having in its bulk 161 and 41 nodes in the axial and radial directions, respectively. Given that the aspect ratio has been decreased by a factor of 2.5 from that in Figure 4, this mesh has the same density as M3 in Table 1. Subsequently, the local mesh refinement technique was

employed, as follows: Around the bubble surface and along the axis of symmetry these meshes were refined 3 times at the radial locations: Rc/5, Rc/10, and Rc/20. Also next to the outer filament surface they were refined once at the radial location 9 Rc/10. The kinematics of the flow is included in the Deborah number through the presence of the pulling velocity in its definition. To examine the effect of elasticity thoroughly, it is useful to separate the elastic from the kinematic effects. This is achieved by introducing the elasto-capillary number, Ec, defined as Ec =

De λ*σ * = * Ca μ*R b,1

(36)

This dimensionless number compares elastic and capillary to viscous forces, while the pulling velocity is now absent from it. Bubble growth and shape are important factors in the adhesive strength of the material. First, the effect of Ec on bubble growth is examined by plotting the bubble deformation parameter D versus time for various Ec, Ca, and De numbers and the same aspect ratio Λ = 1 in Figure 6. We focus on the case in which the bubble is placed initially at the midplane of the filament, that is, at the height h = H0/2. The bubble deformation parameter D is defined as

D=

A−B A+B

(37)

Table 4. L2 Norm of the Difference between the Solutions Obtained with the Single Domain (SD) and Domain Decomposition (DD) Methods for the Radial and Axial Coordinates of the Bubble and Outer Filament Interface for Ca = 100, De = 100, Λ = 1 filament

bubble time

r-direction

z-direction

r-direction

z-direction

3.2485 ± 0.0005 5.4985 ± 0.0005 7.7485 ± 0.0005

7.978 × 10−4 2.636 × 10−4 1.377 × 10−3

1.185 × 10−4 4.127 × 10−5 1.172 × 10−4

1.403 × 10−5 9.045 × 10−5 2.691 × 10−4

9.409 × 10−5 3.693 × 10−5 1.119 × 10−4

7556



Article

Table 6. L2 Norm of the Difference between the Solutions Obtained with the Single Domain (SD) and Domain Decomposition (DD) Methods for the Radial and Axial Coordinates of the Bubble and Outer Filament Interface for Ca = 1, De = 1, Λ = 0.2 filament

bubble time

r-direction

z-direction

r-direction

z-direction

4.399 ± 0.001 12.199 ± 0.001 16.659 ± 0.001

9.793 × 10−4 8.995 × 10−4 1.195 × 10−2

4.943 × 10−4 2.491 × 10−4 2.632 × 10−3

1.881 × 10−4 1.807 × 10−4 3.113 × 10−3

1.248 × 10−4 8.040 × 10−5 2.832 × 10−3

equation. The lowest deformation occurs for the largest Ec number (Ec = 100). The deformation parameter in this case retains rather low values throughout and the changes in the shapes of the bubbles do not exceed 11% of their initial state. There is of course an intermediate case where the bubbles seem to deform nearly linearly with time at least for a major part of the simulation. The curve for De = 5 and Ca = 0.5 approximates better the linear behavior, and so represents the case where a reversal occurs for the other curves from concave to convex. In Figure 7a, snapshots of the instantaneous filament and bubble shapes are presented for various Ec, obtained while keeping the capillary number constant at Ca = 100 and changing De. All the other parameters are kept the same and are given in the figure caption. The filament shapes are in scale, showing that its initial radius and height are 4 times the bubble radius. The very large capillary number is expected to allow large bubble deformations. Two snapshots are shown for each Ec. The first one is chosen for comparison purposes, and corresponds to the same dimensionless time for all three cases presented in the graph; that is, the deformations are all equal (100% in Figure 7a) because the dimensionless pulling velocity is always unity. The second one is the final converged solution in each case. In the first two snapshots the elasticity number is rather low Ec = 0.01, so viscous forces are dominant. The bubble even at early times (small deformations) acquires a prolate ellipsoidal shape, while it gradually gets elongated in the axial direction following closely the motion of the upper disk. Moreover, it preserves its symmetry with respect to the instantaneous midplane of the filament approaching slowly the lower stationary plate. This simulation has reached very high deformations, up to 747%. In this final stage, the bubble is narrowly confined by the surrounding liquid, which, in turn, is extremely thinned because of mass conservation and the filament turns into a fibril. In the second group of snapshots Ec = 0.1, one can observe that the bubble, in the same filament extension as in the previous case, 100%, is somewhat larger, attaining though a similar ellipsoidal shape. In the last case for this figure, Ec = 1, the bubble is elongated, while in the end of the simulations and in the midplane of the filament it approaches the axis of symmetry, so much that it is anticipated that subsequently it would collapse in its middle forming two equal bubbles. We stopped the simulation at the point shown, since convergence could not be obtained at the next time step. Calculating the collapse of the bubble would require resolving the created singularity, describing the filament dynamics containing two bubbles from this point on, and, more importantly, defining an initial state for the two bubbles somewhat arbitrarily. So this was not pursued. We examined further, the results for Ec = 0.1, because here the simulations diverge earlier and the bubble is wider than in the other 2 cases with the same Ca. By comparing the contours of the stress components between these three cases with Ca = 100, we find that those with Ec = 0.1 increase faster and to larger values and that the two normal stress components start having wiggles

Table 7. Number of Elements along Either the Bubble or the Filament Outer Surface for Ca = 1, De = 1, Λ = 0.2a time →

for all times

t = 10−5

method →

SD

DD

DD

DD

DD

80

96

188

304

360

80

160

200

316

372


t = 4.399 t = 12.199

t = 16.659

a

Here the initial number of nodes and triangular elements in the entire liquid volume is (for SD) 35270 and 17280, respectively, and (for DD) 58237 and 28800, respectively.

Figure 6. Evolution of the bubble deformation parameter D in time for a single bubble initially placed at the midplane of the filament, (h = 2 and H0 = 4) for a PTT liquid with ε = 0.03, Λ = 1, and various Ec and Ca numbers.

where A and B is the major and the minor axes of the bubble, respectively. Larger positive D values denote a highly extended bubble shape in the axial direction. In general, bubble deformation is affected primarily by the elastocapillary number, less by the capillary number, and even less by the Deborah number. This becomes evident in Figure 6, where the curves are grouped together if they share the same Ec, while this grouping is weaker if they share the same Ca and almost ceases to exist if they share the same De. One can clearly see that decreasing Ec, that is, increasing the importance of viscous forces, causes the deformation to increase even at short times. For a group with the same Ec the larger deformation occurs for the largest Ca, that is, when the capillary forces get weaker, and for the largest De number, that is, the highest relaxation times. The smaller capillary forces allow the free surface of the bubble to deform more readily in response to the forces applied by the surrounding liquid, while the higher elasticity of the liquid allows the elastic stresses, which tend to resist deformation in the liquid, to develop slower. This should have been expected, since De multiplies the time-derivative in the constitutive 7557



Article

extension of the bubble, which for Ec = 10; Ca = 10 acquires a lemonlike shape. In general, the larger the De number is, the slower the elastic stresses grow. In this very early stage, the capillary forces dominate the elastic ones, forcing the bubble to retain its initial shape, while higher polymeric stresses develop first at the poles, where the bubble interface is locally elongated earlier. Two other groups of shapes for Ca numbers (1, 0.1) are shown in Figures 7 panels c and d, respectively, for a variety of Ec. What we observe here is the gradual restriction of the bubble to its spherical shape with increasing Ec or decreasing Ca. The final cases in these two panels, Ec = 100, are the limiting ones that we examined with our simulations. Now the bubbles remain spherical despite the squeezing applied by the filament’s outer surface, which tends to form a catenoid due to the very strong capillary forces. These cases are the ones shown in Figure 6 to have small time-variations of the deformation parameter D. It is interesting to examine the bubble shapes with respect to Ec. Observing the snapshots chosen to correspond to the same dimensionless time in Figures 7(a−d) (the snapshot in the left for each case), we see that in general, there are three different possible bubble shapes for the aspect ratio, Λ = 1. An elongated bubble of large volume with rounded edges is created for small Ec (Ec = 0.01, Ec = 0.1), a smaller or less elongated one, possibly with sharp poles, arises for intermediate values (Ec = 1, Ec = 10), and finally a small nearly or completely spherical one arises for large Ec (Ec = 100). Clearly, it is the interplay of elastic and capillary forces versus viscous forces that lead to these different shapes. Small Ec results from small elastic and/or capillary forces. For such parameter values the bubble deforms and elongates easier and faster. Large Ec results from large elastic and/or capillary forces. The latter keeps the bubble spherical from the beginning of the filament extension, while the former is activated much later (the larger De multiplies the timederivative in the constitutive model) and reinforces the latter. The cusps at the bubble poles that appear for intermediate Ec are caused by the earlier increase of the elastic forces (at intermediate De) and the stronger squeezing of the bubble by the outer filament surface (at intermediate Ca as well), while it is still close enough to the plates. This combination increases the axial velocities and velocity gradients at the bubble poles, forming the cusp. Indicative velocity fields for such a case are shown in Figure 8. Cusp formation in the buoyancy driven rise of a bubble in a viscoelastic fluid is a well-known and documented phenomenon as reported in Hassager,49 Liu et al.,50 Pilz and Brenn,51 Leal et al.,52 and Malaga and Rallison.53 Here it is predicted for the first time in a filament undergoing extension. On the other hand, the outer surface of the filament is not affected as much by these variations of the physical parameters. This becomes clearer if the outer surfaces are plotted in the same figure for various Ec (not shown here). We note however, that this surface approaches a spherical segment, when capillary forces increase and elastic forces decrease, whereas in the other end of the examined values for both forces, they become flatter in their middle portion with a sharp radial increase closer to the plates. Pressure sensitive adhesives are characterized by high elongational viscosity and very low “solvent” concentration. This is the reason we generally keep the ePTT model parameters, εPTT or β, to a small value. For completeness, we present briefly here the effect of changing either εPTT or β from their basic values. Figure 9 depicts the effect of εPTT on cases b2 and a3 presented in Figure 7. In both cases, the first set of

Figure 7. Snapshots of the instantaneous bubble−liquid and liquid−air interfaces for various Ec numbers and (a) Ca = 100, (b) Ca = 10, (c) Ca = 1, and (d) Ca = 0.1. The remainder of the dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, Λ = 1, H0 = 4, h = 2.

near the bubble surface, which cause the earlier convergence failure of our calculations. Their increase to larger values is caused by the larger De with respect to the case with De = 1 and their faster increase is caused by the smaller De with respect to the case with De = 100. Moreover, the increased elastic forces keep the outer filament surface flatter in comparison to the other two cases where it is more curved. The flatter interface allows the bubble to be thicker for De = 10. Figure 7b shows another family of shapes with smaller capillary number, Ca = 10, and with increasing De by a factor of 10, as previously. We observe here the onset of cusps at the two bubble poles at deformations above 100%. The larger the De number is, the earlier the cusps are formed, so that for Ec = 0.1, that is, for the lowest elasticity, such cusps are formed at very large deformations, where the filament has become an elongated fibril. Moreover, the higher the De, the smaller the 7558



Article

Figure 8. (a) Axial and (b) radial velocity contours for De = 10, Ca = 10 at 166% deformation. The remainder of the dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, Λ = 1, H0 = 4, h = 2.

Figure 9. Snapshots of the instantaneous bubble−liquid and liquid−air interfaces for (a) Ca = De = 10 and (b) Ca = De = 100 and different εPTT values. The remainder of the dimensionless parameters are Re = 0, β = 0.01, Λ = 1, H0 = 4, h = 2.

Figure 10. Snapshots of the instantaneous bubble−liquid and liquid− air interfaces (a) Ca = De = 10 and (b) Ca = De = 100 and different β values. The remainder of the dimensionless parameters are Re = 0, ε = 0.03, Λ = 1, H0 = 4, h = 2.

shapes is at 75% extension and the second set is when the code ceases to converge. We observe that for Ca = De = 10 increasing εPTT by a factor of 40 produces bubble and filament shapes that are almost identical, the only (minor) effect is that the code diverges earlier. Similarly, for Ca = De = 100 increasing εPTT up to 0.1 has a minor effect, but increasing it further to 0.2, which is fairly large for PSAs, produces qualitatively different results: At 75% elongation, the bubble is quite smaller and the filament shrinks more. Additionally, when the code fails to converge, the bubble is more elongated and tends to form cusps at its two poles. The physical explanation for this is as follows: Since Ca = 100, capillary forces are very small compared to viscous forces, but, when εPTT = 0.2, the ePTT model predicts strong shear and elongational thinning, which decreases the effective Ca, which, in turn, prevents the bubble from expanding, as it did for the two smaller εPTT values. The effect of increasing β by up to 50 times from its base value of 0.01 is depicted in Figure 10. The same two cases of Figure 7 are examined. Very close

examination of the 75% elongation for Ca = De = 10 reveals that this very large increase in β results in a slightly more elongated bubble, but has no effect on the outer surface of the filament. On the other hand, the much stronger Newtonian component of the solvent allows the computations to proceed to larger extensions, finally producing thinner bubbles with sharper cusps. Increasing β for Ca = De = 100 decreases the bubble size at 75% elongation for the following reason: Increasing the contribution of solvent viscosity not only decreases the elastic forces (see eq 4), but also increases shear and elongational thinning through exponential factor of the PTT model. Apparently the latter effect dominates leading to smaller viscosity and capillary number and, finally, smaller bubbles. Larger extensions, here shown even before simulations fail to converge, lead to similarly elongated bubbles without cusps on their poles. 4.2.1.2. Small Aspect Ratio, Λ = 0.2. In the probe-tack test of PSAs, the filament aspect ratio is much smaller than unity. To examine such cases, we used Λ = 0.2, that is, now the disk 7559



Article

radius is five times its height, which, in turn, is four times the radius of the bubble. The bubble is again placed initially at the midplane of the filament. After various convergence tests, all calculations with one bubble and Λ = 0.2 were conducted with a mesh having in its bulk 81 nodes in both the axial and radial directions. Subsequently, the local mesh refinement technique was employed, as follows: Around the bubble surface and along the axis of symmetry these meshes were refined twice at the radial locations: Rc/6 and Rc/16. Also next to the outer filament surface they were refined once at the radial location 7 Rc/8. These result in 8640 triangular elements and 111022 nodes in the filament. A sequence of snapshots just as in the previous case of Λ = 1 is presented, keeping in each triplet the same Ca and changing De, resulting in a different in each case Ec. In Figure 11a the Ca

Newtonian one. As the De, and hence the Ec, increases by a factor of 10 or 100 the bubble becomes quite wider with nearly straight side surfaces for the same filament extension. In general, the bubble rapidly expands at early times in both the radial and axial directions. This is caused by the fact that a small upward translation of the upper plate, δh, here draws upward a comparatively much larger volume of liquid πδhR2c , which in combination with the no slip condition on the fluid at the disks forces the bubble to expand substantially in order to conserve the liquid mass. Similarly using ideas from lubrication theory, one can say that the upward motion of the upper plate, induces in the liquid a much smaller pressure closer to the symmetry axis than on its free surface. This sharp pressure decrease in its surroundings allows the bubble to expand radially and grow axially forming a thin layer of liquid between its surface and the upper/lower disk or the free filament surface. For the same large Ca number a larger bubble and with straighter side surfaces arises, when De increases. These changes in the bubble shape are caused by the higher elasticity, which forces the outer filament surface to remain closer to its initial straight shape for most of its length. Again increasing the De number will lead to the slower growth of the stresses in the liquid. Thus the liquid filament is more easily elongated and exhibits smaller (negative) pressures throughout the domain. For the first two cases of Ec numbers (0.01, 0.1) the simulations reached extremely large deformations, 3960% and 4125%, respectively. Most likely they were terminated because of large element distortions. For Ec = 1, they did not converge just after the 500% deformation given in this figure. For this case of Ca = 100, the rather low capillary forces remain subdominant to the grown stresses in the liquid. Conservation of mass causes the outer surface of the filament to move rapidly toward the axis of symmetry, causing the filament walls to become extremely thin in its middle section as well as near the two plates. In the second case that is shown in Figure 11b, the capillary number is 10 times larger, Ca = 10. We consider first the intermediate case with Ec = 1 and observe that the bubble becomes largest in comparison to the other two cases for the same filament extension of 500%. The case with Ec = 0.1 shows that decreasing fluid elasticity allows capillarity to dominate and enforce larger curvature on the filament free surface. This in turn pulls inward the liquid film and the bubble surface around its midplane. In the other end of the spectrum shown in this triplet, Ec = 10, the larger De(= 100) delays the development of the elastic stresses, again giving the bubble a slightly smaller size and higher curvature on its side surface than in the middle case. At the maximum filament extensions we reached with our simulations for Ec = 0.1, 1, the bubble has nearly collapsed in its midplane, whereas for Ec = 10 the filament has nearly collapsed at the same location. Both events force the simulations to terminate earlier than in Figure 11a, with weaker capillarity. Figure 11c illustrates the final triplet of cases with Ca = 1. We observe that for the first two cases with Ec = 1 and Ec = 10, the bubble grows in both directions in a similar manner. At the intermediate filament extensions (150%), the bubble attains a nearly cylindrical shape with rounded corners while at the final extensions reached it collapses at its midplane. As the filament is being stretched, the strong capillarity and mass conservation force the outer free surface to acquire a large curvature, which in turn pushes the liquid inward. Clearly, this does not hold for the third case with Ec = 100. Now the bubble retains its original small and spherical shape in the extending liquid. Instead, it only translates upward along the axis of the filament remaining

Figure 11. Snapshots of the instantaneous bubble−liquid and liquid− air interfaces for various Ec numbers and (a) Ca = 100, (b) Ca = 10, and (c) Ca = 1. The remainder of the dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, Λ = 0.2, H0 = 4, h = 2.

number is set to 100. Just as in the previous section we show first the same filament deformation, here at 500%, for all three cases for comparison purposes. In the first case, Ec = 0.01, and in the first snapshot the bubble attains an almost cylindrical shape with rounded corners and is slightly squeezed in around the midplane by the highly inwardly curved free surface of the filament. Here the capillary and the elastic effects are rather small and the behavior of the viscoelastic liquid approaches the 7560



Article

in its midplane, due to the absence of liquid inertia and gravity. This behavior manifests the domination of the capillary forces on the bubble from the onset of the flow, despite the very large elastic forces (De = 100), which are very much delayed. 4.2.2. Axial Force on the Upper Plate. A very important parameter for the characterization of PSAs is the axial tensile force applied on the upper plate. Its dependence on the strain (or equivalently on time) can offer very useful conclusions for the adhesive properties of the material and the debonding mechanisms of a sample. The important interplay between the time-variation of the applied force and the debonding mechanisms is reported in Lakrout et al.,3 Gay and Leibler4, Brown et al.,5 and Crosby et al.6 When capillary forces on the outer filament surface are strong enough, they will affect significantly the pressure inside it and, hence, the applied force to extend it. To set to zero the applied force prior to filament extension, we need to account for the hydrostatic pressure caused by capillarity. Hence the force applied on the upper plate is calculated directly from the flow field by integrating the normal stresses at z = H(t) and by adding to it the hydrostatic component, πRc/Ca: F̲ = 2π

∫0

Rc

e̲ z ·( −PI + Σ + 2D)|H(t )r dr + π

Rc e̲ z Ca (38)

The force is made dimensionless by the viscous stress scale and the length with the bubble radius, as discussed in section 2. In general, the time-variation of the force exhibits a maximum at early times, followed by an abrupt decrease and a plateau at intermediate times, finally approaching zero when the filament tends to break. The larger the area below this curve is, the larger the adhesion energy during the tack experiment is. Results will be given for some representative scenario of the previously studied cases of large and small aspect ratio. Each figure corresponds to one triplet of cases with three different De and the same Ca, given in Figures 7 and 11. We first present results for Λ = 1 and Ca = 100, in Figure 12a. The force for De = 1 starts from its initial zero value, grows sharply, and exhibits a strong maximum at early times, and then it decreases monotonically reaching asymptotically a value close to zero at large times, when strain has exceeded 700%. This bell-type curve is typical for adhesive materials. As De increases the maximum in the curve decreases and seems to arise at later times. Unfortunately, the evolution at later times is not possible to monitor with our present code. As previously explained, when De increases, the elastic stresses grow slower and extensional thinning is induced in the vicinity of the disks, which decreases the applied force. The applied force changes only a little when Ca decreases by a factor of 10, so we next present the case with Λ = 1 and Ca = 1 in Figure 12b. Of course, with increasing capillarity, one would expect the required force for extension to increase. This is evident if we compare the maximum force for De = 1, which is ∼26.7 when Ca = 100 and increases to ∼40.8 when Ca = 1. By increasing the De to 10 and then to 100, we can see now that these curves have, respectively, a plateau or a maximum which however is weaker than that for the smaller De, as explained previously. These curves for the force correspond to Figure 7c, where it was seen that the bubble tends to maintain its initial spherical shape. When the capillary force gets much stronger, Ca = 0.1, all three cases of De described in Figure 7d exhibit a clear maximum with values in the range of 85−125, with time variations resembling those already described. Clearly, the

Figure 12. Effect of elasticity on the evolution of the force calculated on the upper plate over time for various Ec numbers and (a) Ca = 100 and Λ = 1, (b) Ca = 1 and Λ = 1, and (c) Ca = 10 and Λ = 0.2. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 4, h = 2.

required force increases more abruptly and the effect of varying De decreases when Ca gets smaller. Capillarity prevents the filament and the bubble surfaces from extending, acting similar to a membrane resisting extension and, thus, increase the resistance of extending the filament. Comparing the results in Figure 12 panels a and b, we observe that simulations with different Ca and De but corresponding to the same Ec do not necessarily result in similar force profiles. 7561



Article

shape, and it is deduced that here it finally shrinks below its initial value. 4.3. Single Bubble below the Midplane. It has been discussed elsewhere3 that the onset of bubble cavitation occurs mostly near the lower substrate or onto it, and the study of this particular phenomenon is of great importance. So next we discuss cases with large or small aspect ratios, in which the bubble is initially placed close to the lower substrate. 4.3.1. Large Aspect Ratio, Λ = 7/4. In Figure 14a bubble shapes for a viscoelastic liquid with Ec = 0.01, Ca = 100 are shown for five time instances, t = 1 × 10−5, 1, 3, 5, and 7. The initial filament aspect ratio is Λ = 7/4, resulting in an initial filament height of 7, given that its radius is 4, while the bubble center is at h = 1.5. After various convergence tests, all calculations with 1 bubble and Λ = 7/4 were conducted with a mesh having in its bulk 241 and 81 nodes in its axial and radial directions, respectively. Subsequently, around the bubble surface and along the axis of symmetry this mesh was refined twice at the radial locations: Rc/5 and Rc/10. Also next to the outer filament surface it was refined twice at the radial locations 4 Rc/5 and 9 Rc/10. These result in 44880 triangular elements and 570397 nodes in the filament. The bubble loses its spherical symmetry at early times and gets elongated in the positive z direction. Although it translates, its upper pole is prevented by the thinning filament from ever reaching the plane of symmetry, which at the last instant shown is located at H = 10.5, while its lower pole is only slightly displaced toward the lower disk. The low capillarity and the immediately effective elastic force allow the bubble to deform quite freely, subjected simultaneously to the squeezing motion of the surrounding liquid, which is largest in the upper part of the bubble due to the strongest filament shrinking at the plane of symmetry, while its lower end is “protected” by the lower plate. As a result, the bubble evolves into a teardrop-like shape. In Figure 14b bubble shapes at the same five time instances are shown as in the previous case for a viscoelastic liquid, but with Ca = 10 and so Ec = 0.1. Now, the stronger capillary force prevents the larger deformations and hence the extensive axial growth of the bubble. Consequently, its lower pole translates less toward the lower plate and its upper part is less elongated, although the region near the upper pole is more pointed. The reason for this is the following: The larger capillary force changes the outer filament surface from being nearly flat for most of its height and sharply increasing toward the disk edge near the disk to become more rounded along its entire height and attain a smaller radius at the filament middle plane, which locally squeezes the bubble more. In Figure 14c bubble shapes are given in six time instances, t = 1 × 10−5, 1.00, 1.99, 3.50, 4.99, and 5.56 for a viscoelastic liquid with increased De = 10, but the same Ca = 10. The bubble deviates early from its initial spherical shape and its upper pole translates upward, while the lower one remains almost stationary. This axial elongation of the upper part of the bubble is accompanied by its radial contraction. In the last time instant one can see the very sharply pointed upper pole. Clearly this cusp must have been generated by the increased elastic forces, which pull the bubble stronger upward, although they arise with a longer delay as explained already. Increasing the elasticity further to De = 100, still keeping Ca = 10, the bubble shapes change drastically. The bubble initially translates along the axis of symmetry and only at later times does it get slightly deformed, mostly at its upper part. The cusp does not arise up to the point we obtained converged solutions. The reason for this is that the very large

The force increases dramatically when the aspect ratio decreases to Λ = 0.2. This is shown in Figure 12c for Ca = 10. Especially when the elastic forces are in effect from the beginning of the extension, De = 1, the maximum force is as high as ∼2585. Subsequently, the force falls smoothly to lower values and asymptotically approximates a nonzero, but much smaller value at long times. As explained in Foteinopoulou et al.,28 capillary effects become unimportant for very small aspect ratios. Instead the viscous and elastic forces dominate the flow kinematics. In fact for Newtonian fluids in the absence of bubbles and neglecting capillarity, it was shown that the only contributing factor to the exerted force is the pressure and that it varies as Λ−3. So when Λ = 0.2 (Figure 12c), the force is expected to increase by a factor of 125, over its value when Λ = 1 (Figure 12a). The smaller increase we predict here is attributed to the fluid elasticity and the presence of the bubble in the filament. When De increases to 10 and then to 100 the entire curve shifts to much smaller values, because of significant shear thinning in the liquid and the delayed application of the elastic forces. Finally, all three curves are an example of a relatively good adhesive, since the force does not go directly to zero but is stabilized in an apparent plateau. 4.2.3. Volume of Bubble. Another interesting point that we examined with our simulations is the evolution of the bubble volume in time for a sequence of parameters. This dependence is shown in Figure 13 for various (Ec,Ca) pairs and the aspect

Figure 13. Time evolution of the volume of a single bubble placed initially at the filament midplane for various Ec and Ca numbers. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 4, h = 2.

ratio of Λ = 1. One can see that the curves except for the first case present qualitatively the same shape: they increase attaining a maximum and then most of them descent to lower values. In the first case the bubble, at least for the time interval the simulation proceeded, as shown in Figure 7a for Ec = 0.01 is elongated with a volume that changes linearly in time. In the other cases one can see that as Ec increases, the variation of the bubble volume decreases, and the curves are monotonically flattened. For the same Ec the volume in general is higher for the higher elasticity and the weaker capillarity. We observe that the case represented with the solid curve for Ec = 100 is a limiting case where the volume slightly increases, and then it gradually decreases toward low values, lower than the initial one, that is, 1.33π. In this case the bubble retains its spherical 7562



Article

Figure 14. Bubble interface profiles at five time instances, t = 10−5, 1, 3, 5, and 7 for a viscoelastic liquid with (a) Ec = 0.01 and Ca = 100, (b) Ec = 0.1 and Ca = 10, (c) at time instances, t = 10−5, 1, 1.99, 3.5, 4.99, and 5.56 and Ec = 1 and Ca = 10 and (d) at time instances, t = 10−5, 1, 2.99, 5, and 6.65 and Ec = 10 and Ca = 10. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 7, h = 1.5, and Λ = 7/4.

Figure 15. Contours of the (a) axial and (b) radial velocity at t = 4.99 for a viscoelastic liquid with Ec = 0.1 and Ca = 10. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 7, h = 1.5, and Λ = 7/4.

De allows the elastic forces to increase at much later times than the ones we depict in Figure 14d. It would be interesting to examine the above cases from the kinematic point of view. In Figure 15 panels a and b we present

contours of the axial and radial velocity components, respectively, for the viscoelastic liquid with Ec = 0.1 (De = 1 and Ca = 10) at t = 4.99 or ∼71.4% increase in the filament height. The maximum value of the axial velocity is 1, that is, the 7563



Article

Figure 16. A close up at the region of the bubble with contours of the radial τrr,p and the axial τzz,p normal stress at two time instances t = 3.99 and t = 5.49 with De = 10 and Ca = 10.

5.49 for a viscoelastic liquid with De = 10 and Ca = 10. The smoothness of the contour lines even close to the singular (cusp) points attest to the high quality of our simulations. In Figure 16a and in the lower half of the bubble τrr,p is small and positive and thus the bubble interface expands slightly. In its upper part it attains small compressive values, except for its upper pole, where it turns sharply positive. There it takes its largest positive value 0.362. As for τzz,p shown in Figure 16c, it takes small negative values all around the bubble, except again at the upper pole, where it sharply increases to 0.471. In the second time instant t = 5.49, Figure 16 panels b and d, the cusp in the upper pole is well formed. We can see that the filament interface has come closer to the bubble squeezing it mainly at its upper part. There the compressive radial stresses have become more negative, while on the upper pole they sharply increase in magnitude to −12.323. In the same region the axial normal stress takes very large values of the same order. This combination of radially compressive and axially extensive large normal stresses creates the cusp in the bubble. These locally large stress components grow faster in material with lower elasticity, resulting in larger bubble deformations for the same capillarity. This can be clearly observed by plotting the axial stress component along the bubble interface at t = 3.99 in Figure 17. Indeed, the curve with the lower De(= 1) presents a similar monotonicity with the one for the intermediate De(= 10), for which we have seen that the bubble shapes are also similar, but still it is above it. In both these cases, the stress values present a local maximum around the equatorial plane of

pulling velocity of the upper plate, while the minimum value is −0.0028, very close to the lower pole of the bubble. The interval between these two extremes is divided by 20 contour lines. We observe that these lines are almost straight and equidistant for most of the axial extent of the filament, intersecting the axis of symmetry normally. Hence layers of fluid elements are lifted upward due to the motion of the upper plate. However, in the region where the bubble exists, the contour lines deviate from being straight. The bubble rearranges its shape and gets elongated following closely the motion of the adjacent liquid. The fact that there is an almost linear dependence of the axial velocity on the axial coordinate and the upper disk is being pulled with constant velocity implies that the axial velocity at the same distance from the bottom disk will decrease monotonically with time. The contours of the radial velocity are shown in Figure 15b. Its maximum value is 0.001 and its minimum one is −0.150. Ten contour lines are equally distributed between these two values. Small and negative radial velocities appear at the filament interface. For the upper part of the filament they increase toward zero on the axis of symmetry. Moreover, it admits zero values on the upper and the lower plate surface due to the no slip condition. As for the region around the bubble, the radial velocity takes also negative values, which are in magnitude larger on its upper part than the ones in the lower part, for reasons explained in relation to Figure 14b. The elastic stresses τrr,p and τzz,p, affecting the bubble shapes, are illustrated in Figure 16, at two time instances t = 3.99 and 7564



Article

other hand, the upper pole continuously rises and follows the motion of the upper plate, while it replenishes the relatively large amount of liquid that initially existed above it. One can see that at large times the curved sides of the bubble become almost straight and extend about 9 radii from the filament axis. The shapes of the free surfaces and the pressure field in the liquid are given in Figure 19, at times t = 5 and t = 15 for a viscoelastic liquid with De = 10 and Ca = 100. We observe that the pressure takes negative values throughout the domain. Starting from the filament outer free surface and close to the midplane, it takes small values close to zero, which decrease as we move inward and up to about the middle of the filament width where a saddle point arises. Thereafter, pressure increases radially as we approach the bubble surface where it gets nearly zero, but decreases axially as we approach either one of the disks. Fifteen contour lines divide the interval between the pressure maximum −0.038 and its minimum −0.357 for t = 5 and the same number of contour lines is used in the next time instant. The negative pressures shown are responsible for the radial expansion of the bubble interface. It seems that the low pressure overcomes the elastic stresses in this extensional type flow. Moreover capillarity is weaker than the pressure in this case, so the filament outer free surface acquires shapes of rather low curvature. Indeed, in the second snapshot it is seen that the filament free surface approaches a straight line, except for its part close to the disks, where it bends to reach the edge of the plates due to the no slip condition. From t ≈ 6 onward, the liquid in the vicinity of the two interfaces and away from the region of the bubble poles follows an almost pure axial elongation and retains an annular shape. Pressure is minimized very close to the disks and at a distance from the axis of symmetry which equals about 2/3 the disk radius. In Figure 20 we present a comparison between two cases of viscoelastic liquids, De = 1 and De = 10 for the same capillary number Ca = 100 and the aspect ratio Λ = 1/3. Bubble interfaces are shown in three time instances, t = 1.51, 2.40, and 4.49. At the first time instant, the deviation of the interfaces between the two materials is very small. Its maximum is about 7% at the equatorial plane, while they nearly coincide at their poles. This deviation increases with time, becoming 11% at t = 2.40 and 16% at the last time instant. As already explained, the larger De causes a delay in the growth of the elastic stresses, the absence of which allow the bubble to expand easier in response to the surrounding negative pressures in the liquid. 4.4. Axisymmetric Case with Three Bubbles. Finally, we present results of our simulations with three bubbles in the volume of the filament. Now the mesh M3, given in Table 1, is used throughout. In Figure 21 we show the case of a viscoelastic liquid with De = 1 and Ca = 10 at a relatively large aspect ratio, Λ = 2.5, that is, the dimensionless filament height is 10. The bubbles are placed symmetrically with respect to the filament midplane. Initially, the centers of each spherical bubble are located at the heights h0,1−3 = 2, 5, and 8 from the filament bottom and are aligned along the axis of symmetry. The filament deformation in this figure is 135%. Along with the interfaces we show pressure, radial velocity, and axial velocity contours, respectively. Following the filament extension, the bubbles are elongated and form cusps at their poles that are away from the two disks. Thus, the lower and upper bubbles form cusps in their upper and lower poles, respectively, attaining the shape of a tear drop and remaining symmetric with respect to the instantaneous midplane. The middle bubble forms cusps on both of its poles. These shapes correspond very

Figure 17. Curves for different Ec but Ca = 10 of the axial stress component τzz,p along the interface of the bubbles vs z at t = 3.99.

each bubble and very sharp increase close to the upper pole. The third curve for De = 100 is qualitatively different; it presents two local maxima close to the two poles of the bubble, while the values for the stress are negative all along the interface, and present a rapid decrease locally at the poles. Here the capillarity on the bubble interface dominates the relatively low stresses at this stage of the flow, while at later times (not shown here) the curve for the axial stress attains the same monotonicity with the other two curves. 4.3.2. Small Aspect Ratio, Λ = 1/3. When the aspect ratio is small, the flow field is completely different from that in the previous case. Here we examine the case where the center of the bubble is located on the axis of symmetry and initially at the axial position h = 2, while the distance between the disks is Ho = 6, which with a filament aspect ratio of Λ = 1/3 results in a disk radius of Rc = 18. The time-evolution of the bubble interface is shown in Figure 18 at six time instances. The bubble from its initial spherical shape turns first into an oblate ellipsoid. At these early times the lower pole of the bubble was translated slightly toward the lower plate and resides below the axial position z = 0.4 for all subsequent time instances. On the

Figure 18. Bubble interface profiles at six time instances, t = 10−5, 1, 2, 5, 10, and 15 for a viscoelastic liquid with Ec = 0.1, Ca = 100, and Λ = 1/3. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 6, h = 2. 7565



Article

Figure 19. Contours of the pressure at two time instances (a) t = 5 and (b) t = 15 for a viscoelastic liquid with Ec = 0.1, Ca = 100 and Λ = 1/3. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 6, Rb = 1, h = 2.

Figure 20. Comparison of bubble interfaces at three time instances t = 1.51, 2.4, and 4.29 for two viscoelastic liquids: De = 1 and De = 10. The other dimensionless parameters are Re = 0, Ca = 100, ε = 0.03, β = 0.01, H0 = 6, Rb = 1, h = 2 and Λ = 1/3.

Figure 21. Contours of (a) the pressure, (b) the radial velocity, and (c) the axial velocity at time t = 13.5 for a viscoelastic filament with De = 1 and Ca = 10 with three bubbles in its volume. The other the dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 10, h0,1 = 2, h0,2 = 5, h0,3 = 8, and Λ = 2.5.

nicely with those in Figure 7b for a bubble placed in the middle of the filament and Figure 14b,c for a bubble placed in the lower end of the filament. Interestingly, in all cases Ca = 10 and De = 1 or 10. The pressure in the filament presents a small variation with the radial coordinate and a stronger one with the axial coordinate. Its values remain rather low, of the order of −0.01, except close to the poles where cusps are formed. There they increase by more than 2 orders of magnitude, because of the singularity arising there. Between the two types of cusps, the one belonging to the middle bubble attains a somewhat smaller pressure. As explained for the single bubbles, cusps arise when the bubbles do not extend too much (as for intermediate Ca cases) and elasticity is large enough and sufficient time is allowed for it to become effective (as for intermediate De). These large values of the pressure balance the increased elastic stresses in the same regions (not shown here for conciseness). The radial and axial normal stress components are strongly compressive and extensive, respectively, around the cusps subjecting the poles to a radially inward and stronger axial motion. The radial velocity presents a minor axial variation along the filament outer free surface, always assuming negative values corresponding to its inward direction. Away from the

disks, this radial velocity first increases in magnitude toward the filament axis or the bubbles’ surfaces and then decreases to zero at the axis of symmetry. Moreover, it increases locally around the cusps at the bubble poles. In the third part of this figure the axial velocity shows a pure dependence on the axial direction only ranging from the value of zero at the lower plate to the value of unity, that is, the pulling velocity of the upper plate. The elongational nature of the flow transformed the initial cylindrical filament with spherical bubbles into a thin fibril. Despite the presence of the bubbles the contour lines are nearly straight, intersecting the axis of symmetry perpendicularly. On the other hand, in Figure 22 the shapes of the filament and the bubbles and of the flow field for a viscoelastic liquid with a larger Ca(= 100), but the same De = 1 is quite different. The rather weak capillary forces allow the bubble to deform more freely and acquire more elongated shapes, while the bubble poles now are quite curved. The smaller capillarity allows the outer filament surface to change more abruptly from the radial position of 4 at the disk surface and to remain flatter 7566



Article

With our robust numerical tools we managed to successfully simulate the filament stretching of viscoelastic liquids characterized by very high De (up to 100), while the magnitude of the capillary force applied on the liquid/air interfaces varied widely. In the case of a large filament aspect ratio and a bubble initially located at the midplane of the filament, the smaller Ec leads to larger bubble deformations. For the same Ec, the weaker capillarity allows the bubble to deform more. At moderate Ec, the formation of cusps at the bubble poles was observed. These cusps are formed by the interplay of the increased degree of capillarity, which tends to suppress the bubble shape, via squeezing the filament outer surface, and the elastic stresses which must be large enough and develop fast enough and initially at the bubble poles. Moreover, the lower the elasticity is, the larger the adhesion energy becomes; that is, the force vs time presents a higher peak, for the same degree of capillarity. With increasing capillarity, the force required to pull the upper disk is increased as expected. When the bubble is located closer to one of the plates and at large aspect ratios it quickly loses its spherical symmetry attaining less elongated shapes than the previous case. The onset of cusps was also seen here and attributed to the large in magnitude radial compressive (negative) and axial normal stresses. For the case of a filament with small aspect ratio and the bubble located at the midplane, the radial and axial growth of the bubble is observed quite early. Just as in the Newtonian case reported in our previous study this is due to the negative pressures that develop in the flow field. Now we managed to reach very large deformations (over 1000%). Similar is the situation when the bubble is placed initially below the midplane. For the same degree of capillarity the larger De gives rise to larger bubble expansion. The case of three bubbles existing in the filament was also examined. In general, we found that, for the same parameter values, the bubble shapes depend on the bubble location and the filament aspect ratio, but less on the number of bubbles in the filament as long as they are not close enough to interact strongly with each other.

Figure 22. Contours of (a) the pressure, (b) the radial velocity, and (c) the axial velocity at time t = 13.5 for a viscoelastic filament with De = 1 and Ca = 100 with three bubbles in its volume. The other dimensionless parameters are Re = 0, ε = 0.03, β = 0.01, H0 = 10, h0,1 = 2, h0,2 = 5, h0,3 = 8, and Λ = 2.5.

than in the previous case away from the disks. This change of the outer surface causes less squeeze of the bubbles inside the filament and delays the formation of the cusps. The pressure varies in both directions, but takes smaller values than the previous case all over the domain. However, it increases sharply but to a lesser extent than before when the bubble poles are approached where cusps were formed in Figure 21. This means that cusps could arise even in this case if the simulations proceeded further. The radial velocity in this case presents a weaker variation in the axial direction. Ten contour lines divide the interval between its maximum (zero) value at the axis of symmetry and its minimum, −0.062 at the free surface. The flatter outer surface of the filament for this case of Ca gives rise to this weak variation. For the same radius in the cylindrical filament, a fluid element that is located next to one of the bubbles moves radially faster toward the bubble than a fluid element which resides at the level of the bubble pole. In Figure 22c one can observe that the contour lines present a slight radial variation close to the bubbles only not too different from the one in Figure 21c.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been cofinanced by the European Union (European Social FundESF) and Greek national funds through the Operational Program “Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF). Research Funding Program: Heracleitus II of GSRT.

5. CONCLUSIONS We studied the axisymmetric extensional flow of viscoelastic filaments containing one or three bubbles in their volume and placed between two coaxial disks. We examined the cases of large and small aspect ratios, and we tested different positions of the bubble along the axis of symmetry. A domain decomposition method was employed in order to compare and validate our numerical results with respect to those obtained via a single domain method. Very good agreement was found in terms of bubble and filament shapes during large filament deformations. In the former case a new remeshing procedure was developed, during which finite elements are added in the domain, a fact which ensures optimal discretizations in the physical space. Moreover, various tests that were performed ensure mesh convergence and correct time-step size adaptivity.

■ ■

DEDICATION This work is dedicated to the memory of John P. Congalidis, a noble human being, a top engineer, and a real friend. REFERENCES

(1) Creton, C. Pressure sensitive adhesives: An introductory course. MRS Bull. 2003, 28 (06), 434. (2) Zosel, A. Adhesive failure and deformation behaviour of polymers. J. Adhes. 1989, 30, 135. (3) Lakrout, H.; Sergot, P.; Creton, C. Direct observation of cavitation and fibrillation in a probe tack experiment on model pressure-sensitive-adhesives. J. Adhes. 1999, 69, 307.

7567



Article

(4) Gay, C.; Leibler, L. Theory of tackiness. Phys. Rev. Lett. 1999, 82, 936. (5) Brown, K.; Hooker, J. C.; Creton, C. Micromechanisms of tack of soft adhesives based on styrenic block copolymers. Macrom. Mater. Eng. 2002, 287, 163. (6) Crosby, A. J.; Shull, K. R.; Lakrout, H.; Creton, C. Deformation and failure modes of adhesively bonded elastic layers. J. Appl. Phys. 2000, 88, 2956. (7) Gent, A. N.; Wang, C. Fracture mechanics and cavitation in rubber-like solids. J. Mater. Sci. 1991, 26, 3392. (8) Gay, C. StickinessSome fundamentals of adhesion. Integr. Comp. Biol 2002, 42, 1123. (9) Multi-Scale Modelling of Interfacial Phenomena in Acrylic Adhesives Undergoing Deformation (MODIFY); European Collaborative project NMP-2008-2.5-2; Project Coordinator, Vlasis Mavrantzas, Contract No. 228320; June 2009−May 2012, 2012. (10) Pearson, G.; Middleman, S. Elongational flow behavior of viscoelastic liquids. Part I. Modeling of bubble collapse. AIChE J. 1977, 23, 714. (11) Ting, R. Y. Effect of polymer viscoelasticity on the initial growth of a vapor bubble from gas nuclei. Phys. Fluids 1977, 20, 1427. (12) Papanastasiou, A. C.; Scriven, L. E.; Macosko, C. W. Bubble growth and collapse in viscoelastic liquids analyzed. J. Non-Newton. Fluid Mech. 1984, 16, 53. (13) Bousfield, D. W.; Keunings, R.; Denn, M. M. Transient deformation of an inviscid inclusion in a viscoelastic extensional flow. J. Non-Newton. Fluid Mech. 1988, 27, 205. (14) Chen, T. Y.; Tsamopoulos, J.; Good, R. Capillary bridges between parallel and non-parallel surfaces and their stability. J. Colloid Interface Sci. 1992, 151 (1), 49. (15) Tsamopoulos, J.; Chen, T.-Y.; Borkar, A. Viscous oscillations of capillary bridges. J. Fluid Mech. 1992, 235, 579. (16) Chen, T.-Y.; Tsamopoulos, J. Nonlinear dynamics of capillary bridges: Theory. J. Fluid Mech. 1993, 255, 373. (17) Mollot, D.; Tsamopoulos, J.; Chen, T.-Y.; Ashgriz, N. Nonlinear dynamics of capillary bridges: Experiments. J. Fluid Mech. 1993, 255, 411. (18) McKinley, G. H.; Sridhar, T. Filament-stretching rheometry of complex fluids. Annu. Rev. Fluid Mech. 2002, 34, 375. (19) Spiegelberg, S. H.; Ables, D. C.; McKinley, G. H. The role of end effects on measurements of extensional viscosity in filament stretching rheometers. J. Non-Newtonian Fluid Mech. 1996, 64, 229. (20) Yao, M.; McKinley, G. H. Numerical simulation of extensional deformations of viscoelastic liquid bridges in filament stretching devices. J. Non-Newtonian Fluid Mech. 1998, 74, 47. (21) Yao, M.; Spiegelberg, S. H.; McKinley, G. H. Dynamics of weakly strain-hardening fluids in stretching devices. J. Non-Newtonian Fluid Mech. 2000, 89, 1. (22) Spiegelberg, S. H.; McKinley, G. H. Stress relaxation and elastic decohesion of viscoelastic polymer solutions in extensional flow. J. Non-Newtonian Fluid Mech. 1996, 67, 49. (23) Bach, A.; Rasmussen, H.; Longin, P.; Hassager, O. Growth of non-axisymmetric disturbances of the free surface in the filament stretching rheometer: experiments and simulation. J. Non-Newtonian Fluid Mech. 2002, 108, 163. (24) Rasmussen, H. K.; Hassager, O. Three-dimensional simulations of viscoelastic instability in polymeric filaments. J. Non-Newtonian Fluid Mech. 1999, 82, 189. (25) Foteinopoulou, K.; Mavrantzas, V.; Tsamopoulos, J. Numerical simulation of bubble growth in Newtonian and viscoelastic filaments undergoing stretching. J. Non-Newton. Fluid Mech. 2004, 122, 177. (26) Poslinski, A.; Tsamopoulos, J. Inflation dynamics of fluid menisci inside a mold cavity: I. Deformation driven by small gas pressures. Chem. Eng. Sci. 1991, 46 (1), 215. (27) Joseph, D. D. Cavitation and the state of the stress in a flowing liquid. J. Fluid Mech. 1998, 366, 367. (28) Foteinopoulou, K.; Mavrantzas, V.; Dimakopoulos, Y.; Tsamopoulos, J. Numerical simulation of multiple bubbles growing

in a Newtonian liquid filament undergoing stretching. Phys. Fluids 2006, 18, 042106−1. (29) Dimakopoulos, Y.; Tsamopoulos, J. A quasi elliptic transformation for moving boundary problems with large anisotropic deformations. J. Comput. Phys. 2003, 192, 494. (30) Sujatha, K. S.; Matallah, H.; Webster, M. F.; Williams, P. R. Numerical predictions of bubble growth in viscoelastic stretching filaments. Rheol. Acta 2010, 49, 1077. (31) Dimakopoulos, Y.; Papaioannou, J.; Tsamopoulos, J. Filament stretching of pressure sensitive adhesives with gas inclusions: A 3D hyperelastic-viscoelastic model and numerical simulations, 23rd ICTAM, 19−24 August 2012, Beijing, China. (32) Papaioannou, J.; Karapetsas, G.; Dimakopoulos, Y.; Tsamopoulos, J. Injection of a viscoplastic material inside a tube or between two parallel disks: Conditions for wall detachment of the advancing front. J. Rheol. 2009, 53, 1155. (33) Chatzidai, N.; Giannousakis, A.; Dimakopoulos, Y.; Tsamopoulos, J. On the elliptic mesh generation containing multiple inclusions and undergoing large deformations. J. Comput. Phys. 2009, 228, 1980. (34) Nase, J.; Lindner, A.; Creton, C. Pattern formation during deformation of a confined viscoelastic layer: From a viscous liquid to a soft elastic solid. Phys. Rev. Lett. 2008, 101, 074503−1. (35) Carelli, C.; Déplace, F.; Boissonnet, L.; Creton, C. Effect of a gradient in viscoelastic properties on the debonding mechanisms of soft adhesives. J. Adhes. 2007, 83, 491. (36) Thien, N. P.; Tanner, R. I. A new constitutive equation derived from Network theory. J. Non-Newton. Fluid Mech. 1977, 2, 353. (37) Creton, C. Personal communication, 2012. (38) Yao, M.; McKinley, G. H.; Debbaut, B. Extensional deformation, stress relaxation and necking failure of viscoelastic filaments. J. NonNewtonian Fluid Mech. 1998, 79, 469. (39) Rajagopalan, D.; Armstrong, R. C.; Brown, R. A. Finite element methods for calculation of steady, viscoelastic flow using constitutive equations with Newtonian viscosity. J. Non-Newtonian Fluid Mech. 1990, 36, 159. (40) Brown, R. A.; Szady, M. J.; Northey, P. J.; Armstrong, R. C. On the numerical stability of mixed finite-element methods for viscoelastic flows governed by differential constitutive equations. Theoret. Comput. Fluid Dyn. 1993, 5, 77. (41) Dimakopoulos, Y.; Tsamopoulos, J. Transient displacement of a Newtonian and viscoplastic liquid by air ion complex tubes. J. NonNewton. Fluid Mech. 2007, 142, 162. (42) Karapetsas, G.; Tsamopoulos, J. Steady extrusion of viscoelastic materials from an annular die. J. Non-Newtonian Fluid Mech. 2008, 154, 136. (43) Szabo, B.; Babuska, I. Finite Element Analysis; Wiley: New York, 1991. (44) Tsiveriotis, K.; Brown, R. A. Solution of free-boundary problems using finite element /Newton methods and locally refined grids: Application to analysis of solidification mictrostructure. Int. J. Numer. Meth. Fluids 1993, 16, 827. (45) Szady, M. J.; Salamon, T. R.; Liu, A. W.; Bornside, D. E.; Armstrong, R. C.; Brown, R. A. A new mixed finite element method for viscoelastic flows governed by differential constitutive equations. J. Non-Newton. Fluid Mech. 1995, 59, 215. (46) Brooks, A. N.; Hughes, T. J. R. Steamline upwind/PertovGalerkin formulations for convection dominated flows with particular emphasis on thr incompressible Navier-Stokes equations. Comput. Appl. Mech. Eng. 1982, 32, 199. (47) Schenk, O.; Gärtner, K. Solving unsymmetric sparse systems of linear equations with PARDISO. J. Future Generat. Comput. Syst. 2004, 20, 475. (48) Schenk, O.; Gärtner, K. On fast factorization pivoting methods for symmetric indefinite systems. Electron. Trans. Numer. Anal. 2006, 23, 158. (49) Hassager, O. Negative wake behind bubbles in non-Newtonian liquids. Nature 1979, 279, 402. 7568



Article

(50) Liu, Y.; Liao, T. Y.; Joseph, D. D. A two-dimensional cusp at the trailing edge of an air bubble rising in a viscoelastic liquid. J. Fluid Mech. 1995, 304, 321. (51) Pilz, C.; Brenn, G. On the critical bubble volume at the rise velocity jump discontinuity in viscoelastic liquids. J. Non-Newtonian Fluid Mech. 2007, 145, 124. (52) Leal, L. G.; Skoog, J.; Acrivos, A. Motion of gas bubbles in a viscoelastic liquid. Can. J. Chem. Eng. 1971, 49, 569. (53) Malaga, C.; Rallison, J. M. A rising bubble in a polymer solution. J. Non-Newtonian Fluid Mech. 2006, 141, 59.

7569


Bubble Deformation and Growth Inside Viscoelastic Filaments

Recommend Documents