#### Abstract

Analysis of literature data and our own experimental observations have led to the conclusion that, at high deformation rates, viscoelastic liquids come to behave as rubbery materials, with strong domination by elastic deformations over flow. This can be regarded as a deformation-induced fluid-to-rubbery transition. This transition is accompanied by elastic instability, which can lead to the formation of regular structures. So, a general explanation for these effects requires the treatment of viscoelastic liquids beyond critical deformation rates as rubbery media. Behaviouristic modeling of their behaviour is based on a new concept, which considers the medium as consisting of discrete elastic elements. Such a type of modeling introduces a set of discrete rotators settled on a lattice with two modes of elastic interaction. The first of these is their transformation from spherical to ellipsoidal shapes and orientation in an external field. The second is elastic collisions between rotators. Computer calculations have demonstrated that this discrete model correctly describes the observed structural effects, eventually resulting in a “chaos-to-order” transformation. These predictions correspond to real-world experimental data obtained under different modes of deformation. We presume that the developed concept can play a central role in understanding strong nonlinear effects in the rheology of viscoelastic liquids.

#### 1. Introduction

The general task in understanding the rheology of polymeric liquids is constructing models that allow one to describe the behaviour of a material using the minimal number of physical parameters. The first fundamental step along this path was a model of an individual viscoelastic chain placed in a viscous liquid [1], known in the global literature as the Rouse model. This model predicted the existence of a relaxation spectrum (with regular distribution of relaxation times) and gave insight into the correlation of the molecular parameters of a chain and its relaxation properties. Later, numerous modifications of this model were proposed and widely discussed. The results of this line of investigations were collected and generalized in the framework of the concept of conformational statistics of polymers [2]. Further development of the statistical approach for the rheology of concentrated solutions and melts was based on the concept of entanglements, which is a physical network with some characteristic lifetimes (relaxation times) of local contacts, or nodes, formed by macromolecules [3–5]. Numerous versions and modifications of more complex entanglement models, which rather successfully explained characteristic features of polymeric fluids, were also proposed.

A radical retreat from the viscoelastic chain model and the entanglement concept was proposed by Doi and Edwards in their “tube model” and the reptation-like movement of a single chain inside a tube [6]. The physical transparency of this model was the grounds for simple scaling relationships, which proposed simple relationships between molecular parameters and relaxation properties of polymers [7]. This approach was developed in numerous publications and finally describes qualitative and in many cases semiquantitative peculiarities of the mechanical properties of liquid polymers.

The state of the art in this field based on existing concepts of the dynamics of macromolecular chains committing reptation-like motions inside a tube and scaling ideas was outlined in [8], where the universality of many physical properties of polymers has been fairly stressed. This phenomenon “arises purely from the connectivity of their molecular chains rather than from the chemistry of their monomer units.” The possibility of describing the linear viscoelastic properties of complex systems consisting of interacting elementary units with many thousands of segments up to a range exceeding 12 decimal orders by means of a rather simple model was estimated as “actually rather remarkable” [8].

Meanwhile, it was stressed that despite the general qualitative agreement between theoretical predictions and experimental data, several aspects of quantitative descriptions remain unsatisfactory [9]. There have been different attempts to modify the basic theory by introducing additional factors, such as slip-spring elements, which constrain the reptation-like movement of a chain inside a tube and so on (e.g., [10–12]).

However, the principle point is that these model concepts work marvellously only in a linear region of viscoelasticity. In the transition to nonlinear behaviour, a lot of questions and problems appear. The first complex model of nonlinear viscoelasticity, known as the well-known Kaye-Bernstein-Kearsley-Zapas equation, explored the idea of separability (or factorization) of nonlinear viscoelastic properties into a “linear” part and a modifying function. This concept was then widely used by many authors; some versions of this approach include several modifications of the so-called Wagner models. The problem of nonlinear relaxation and its quantitative description has also been discussed by numerous authors. The correspondence of this approach with the Doi-Edwards model has also been examined. This line of studies is presented in several papers [10, 13–16] which cite and summarize the numerous investigations in this field.

The method of factorization in treating nonlinear viscoelastic properties is not general [17], and failure in applications of models of this class is also known, especially when macromolecules of complicated structure are studied [18, 19]. Different approaches to fast flows have been developed based on the stochastic simulation of entangled polymeric liquids [20].

Discussions of the problem of the nonlinear behaviour of viscoelastic fluids have always been limited by the frame of the continuum mechanics. Meanwhile, the mechanism of deformations at high shear rates remains not quite clear. Several publications in the last decade based on direct measurements demonstrated the heterogeneity of flow, with jumps in shear rates. In addition, direct observation showed that self-organization actually takes place in rotational flows of polymeric viscoelastic substances.

Discussions on nonlinearity and self-organization in the deformation of polymeric fluids tacitly assume that a deformed matter reminds fluid at any rate and the issue consists only of the choice of an appropriate constitutive equation. As such, many problems of instability, bifurcation, deformation-induced structure formation, and even fracture are considered in the framework of Giesekus or FENE models (e.g., [21]) or some other continuum models.

Meanwhile, the last statement is of crucial value because it is permissible to ask: is it true at arbitrary high deformation rate (or stress)?

To answer this question, it is appropriate to return to the early conceptual results of Vinogradov et al. [22] who described the “spurt” phenomenon, the transition from shear flow to slip at high enough stresses, as well as the filament behavior at high enough deformation rates [23]. Experimental results of this kind have been treated in a rather evident manner: complex viscoelastic fluids lose the ability to flow at high deformation rates and transit into an elastic (rubbery) state, or, to formulate it more accurately, a rubber-like behavioural state occurs [24].

The dominant component of full deformation in this state is elastic (recoverable) deformations, and real flow becomes negligible. This is analogous (but not equivalent) to the flow-to-rubbery transition that takes place as a function of frequency in oscillation measurements. An obvious reflection of this analogue is the well-known Cox-Merz rule, though it has a different understanding [25, 26].

The separation of the full deformation into elastic (reversible) and flow (irreversible) components is made by measuring the reaction of a sample after cessation of the external force. The partial recovery of the initial shape is observed and this is the elastic component and the other (nonrecoverable) part of the total deformation is attributed as the flow.

The domination by elastic deformations over irreversible flow at high deformation rates is illustrated in Figure 1 for extension. The situation for shearing is quite the same for the initial stage of deformation or high deformation rates, though for low deformation rates stationary flow takes place and irreversible deformations become dominant.

**(a)**

**(b)**

A flow-to-rubbery transition at high deformation rates (with absolute domination of elastic deformation and negligible flow) is especially evident in extension, where drawing takes place by elastic deformations and ends with the breakup of a rubber-like filament [27].

The fact is that the irreversible part of a deformation decreases and the elastic (rubbery) recoil increases along with increasing deformation rates, either in shear or in extension (Figure 1). The final act is either a spurt in shear (transition from flow to slip) [22] or a jet breakup in extension [24, 27]. It should be emphasized that the high rate extension of polymeric fluids is realized by rubbery deformations and irreversible deformations (flow) being negligible. The flow-to-spurt transition was first described for “Newtonian-like” monodisperse polymers [22] but the same phenomenon takes place for polydisperse non-Newtonian polymers (as in Figure 1(a)) because polydispersity leads to non-Newtonian flow [28] but does not exclude the spurt effect.

The concept of the fluid-to-rubber-like transition at high deformation rates seems noncontradictory and rather evident as the generally accepted consequence of different time-dependent modes of rheological behaviour. However, it is necessary to mention the other point of view on the high deformation rate transition. According to the tube model, increasing the deformation rate leads to disentanglement of initially entangled macromolecular chains and this results in a strong decrease in the apparent viscosity [29]. The decrease in apparent viscosity really takes place in the non-Newtonian part of a flow curve, but it proceeds simultaneously with increasing rubbery deformations and can continue inasmuch as a polymer does not lose fluidity by becoming a rubber-like medium.

Meanwhile the entanglement-disentanglement process at high deformation rates is worth special discussion. Macromodeling behaviour of entangled rubber strips as models of flexible macromolecules showed that disentanglement takes place at low velocities, while at high deformation rates entanglement knots are tightened, preventing interchain slippage (i.e., flow) [30]. Correlation of these results with the real behaviour of chains in entanglements depends on the ratio between the elastic energy stored in deformations and the Brownian movement energy. This criterion for a polymer melt for can be written as where is the shear stress corresponding to the deformation-induced flow-to-rubber-like transition (fluid-to-rubbery transition), is the shear modulus in the terminal zone (but not at the plateau), is the average molecular weight of a chain segment between neighbouring entanglements, is the density, is the universal gas constant, and is the absolute temperature.

This criterion can be rearranged by introducing the plateau modulus, :and therefore

The -criterion by virtue of its structure is an analogue of the well-known Peclét number, with the evident substitution of the stored elastic energy instead of the viscous dissipation energy.

Direct calculation shows that, at least in some cases, . This means that elasticity prevents Brownian slippage of macromolecules in entanglements. It leads to formation of chain bundles. Then, a chain must proceed on a much longer path to slip out of all entanglements, which can be treated as the increase in the apparent . Indeed, we would relate the Brownian movement to a much larger volume than a single chain length between neighbouring entanglements. So, one can expect that successful disentanglement would require many more acts of Brownian motion and consequently much higher values should be expected. This picture is also valid for polymer solutions.

Polymeric fluids at high deformation rates behave as rubber-like media. This is also illustrated by the experimental data presented in Figure 2, where dependencies for elastic (rubbery) deformations on the Weissenberg Number, Wi, for linear and slightly cured polyisoprenes are compared. It is evident that these dependencies are quite similar over a wide range of Wi values.

The general conclusion of the above discussed peculiarities of deformations in strong stress fields is that it seems necessary to assume that there is a limit to steady flow, because a polymeric fluid loses the possibility to flow and comes to behave in a rubber-like manner at the Weissenberg Number, Wi, of the order of several units. This fact is the starting point for the model developed below.

A rather provoking statement, “entangled liquids are solids,” sharply formulates the problems existing in this field [31]. This new concept answers the above statement in a quite unequivocal manner: “Yes, entangled matters are ‘solids’* at high deformation rates.*”

The main arguments that were the starting point for development of a new model of polymer behaviour are as follows:(i)Complex viscoelastic liquids (polymer solutions and melts as well as some other soft matters* at high deformation rates* or high values of the Weissenberg Number) behave like solid-like (rubbery) media and the possibility of their flow can be neglected. So, it is possible to neglect dissipative losses in comparison with elastic deformations.(ii)In strong stress fields, tightening of permanent knots (but not disentanglements) takes place with formation of bundles and inhomogeneous structure.(iii)The behaviour of complex fluids at high deformation rates cannot be examined treating a medium as continuum because in reality a transition from the homogeneous flow to heterogeneous deformations takes place, being a result of self-organization including the formation of periodic structures.

These general concepts were the basis for constructing a novel behaviouristic model of complex viscoelastic media behaviour at high deformation rates (stresses). According to this model, a medium at high deformation rates is treated as a set of discrete elastic elements that mechanically interact with each other. In a hypothetical nondeformed (by convention, “initial”) state, these primary elements are spherical particles, but they transform to an ellipsoidal shape under the action of external forces and mutual collisions. We suppose that this approach will allow us to advance into the domain of strongly nonlinear rheological behaviours of complex liquids, where standard ideas of continuum mechanics do not explain or describe the behaviour of these media.

So, this work is constructed in the following manner. First, a novel model will be presented which is based on the arguments cited above. The main point is considering polymer fluids in strong stress fields as elastic media with neglecting dissipative losses. The derivation of the base equations which are used along the whole text is explained in detail in the appendix. The final goal of this section was to receive the equation governing the behaviour of an elastic fluid with limited number of material parameters. Second, the consequences from the model are calculated and discussed. In this part of the text, such principle things as self-arrangements of structure elements and instability of deformations are the focus of discussion. Also, the model predictions are qualitatively compared with experimental data. In conclusion, the main concepts of a novel model and further possibilities of its development are summarised.

#### 2. Model

The image of the model presented in Introduction is shown in Figure 3. This model of behaviour corresponds to a complex liquid (viscoelastic polymer system or micellar colloid) being deformed at a high deformation rate. Figure 3(a) corresponds to the hypothetical nondeformed (by convention, “initial”) state that is achieved at the point of the deformation-induced liquid-like to rubber-like behaviour transition.

**(a)**

**(b)**

An element in the model (“granule”) is assumed to be large enough to neglect Brownian motion. It corresponds to the following limitation on the Péclet Number: where is viscosity, is shear rate, is a characteristic size, is the Boltzmann constant, and is absolute temperature.

Also, the condition should be satisfied that allows us to exclude knot disentanglement.

The most important assumption is the concept of purely elastic deformation of structural elements with negligible dissipation. Actually, it directly reflects the situation observed at a high rate of extension where (beyond some limit) full deformations appear totally elastic [23, 24, 27], as shown in Figure 1.

Upon further increases in deformation rates, events transpire that are illustrated in the left portion of Figure 3(a). The main effect is the formation of knots if entangled macromolecules connected by the sparse network of connecting chains. This topological effect leads to the formation of rubber-like domains and reputation modes are excluded in those domains. In the model under discussion, this situation is simulated by the ensemble of elastic ellipsoids (Figure 3(b)). The ellipsoids are located in the nodes of a squared lattice. Translation modes of movement are forbidden because we consider the situation when the flow is absent due to flow-to-rubbery transition at high deformation rates. The lattice is assumed to be rectangular though that is not of principle value. This is necessary only for the particle numbering.

It was supposed that the free volume in the system is large enough for particle rotations, so particles can collide with surrounding ones and deform elastically. It is assumed that deformations of particles as a result of mutual impacts are small in comparison with the global deformation in the sphere-to-ellipsoid transition. So, the ellipsoid shape of the particles is conserved. Third, the dominating orientation of the principle axes exists and depends on the deformation rate.

Only a single layer is presented in Figure 3, whereas a 3D lattice should be considered, in which elements of one layer interact not only with elements of that layer but also with neighbouring elements in upper and lower layers.

The equation of motion for the structural elements is determined by the Hamilton function for a system of particles positioned in the nodes of a uniform lattice with cells of a fixed size. This is possible due to neglecting viscous dissipation.

According to the model formulation, two items should be considered. The first is related to the energy of orientation interaction of a system of elastic ellipsoids. The centres of particle mass are localized in the nodes of the lattice and the particles have only a rotational degree of freedom. The shape of the ellipsoids is characterized by their eccentricity, , and their spacial position is determined by the vector coinciding with the direction of the main principle axis (), where is the eccentricity of an ellipsoid and are the components of the projectors of the eccentricity at the coordinate exes; if an elastic particle is not deformed, then .

Generally speaking, the deformation of structural elements (being rubber-like) can be large and may be described by any potential function proposed for rubbery deformations (e.g., by the Mooney-Rivlin potential function). However, to simplify the analysis, it is assumed that elastic deformations are described by a single value of the Hooke modulus, . The deformation behaviour of the particles was modelled by means of computer calculations in the* Comsol Multiphysics* medium by the finite-difference method. It allowed us to determine the energy of steric interactions of the model’s structural elements. The analytical approximation for the dependencies of the specific elastic energy on the elastic modulus, , eccentricity and the angle between principle axes of ellipsoids appears as follows (see the appendix for details): where is a scalar product of vectors and and , where the angle is expressed as .

The second item reflects the turning of particles under the action of the external dynamic field. This was also obtained by computer modeling. The obtained results for the specific energy of particle deformations due to rotating of ellipsoidal particles in a shear field can be approximated by the equation (see the appendix for details)where is the shear stress and is a unit vector showing the direction of shear.

Parameters *α* and *β* in (5) and (6) are numerical factors of the order of 1.

In formulating an equation describing particle rotary movements on a lattice, we accept the assumption of the much higher potential energy of the particles in comparison with the kinetic energy of rotation. This means that we consider only elastic interactions without exchange of momentum impulses. Then the Hamilton function can be written as the sumLet satisfy the classic equation of movement where is the Poisson brackets, which are the Lee derivative along the direction set in the vector space :Here is the Levi-Civita symbol.

In calculating the Poisson brackets, it is necessary to take into account the fact that the vector product is ( are unit vectors) and the values do not depend on direction. Then the particle evolution on the lattice can be written as the vector equation

The value in (10) characterizes the input of stress in the shear field, and the value presents the elastic potential of sterically interacting particles.

Equation (10) can be presented in the form of the differential finite-difference nonlinear Schrödinger equationwhere the variables and are connected by the relationsThe nonlinear Schrödinger equationis an integrated equation and its solutions describe the process of spreading of solitary waves or their superposition (multisoliton solutions). These solutions correspond to whirl structures in the hydrodynamics of viscous fluids. It is worth mentioning that the class of problems related to the nonlinear Schrödinger equation is rather wide and has numerous applications, for example, in discussing the stability of a plasma in a magnetic field, the flow of superfluid liquids, the behaviour of spin systems [32], and the impulse flow of a liquid in elastic tubes, for example, blood flow in arteries as pressure waves and conjugated waves of deformations of a liquid propagating along anelastic tube.

Indeed, the nonlinear properties of this equation are of fundamental value in hydrodynamics and presumably in the dynamics of elastic liquids too.

The following feature of the formulated equation is worth mentioning. Equation (7) does not contain a member reflecting the kinetic energy. So, it does not describe the behaviour of real ellipsoidal particles themselves because they have inertia, but momentless quasiparticles, able to be deformed and oriented by action of external field and due to interaction each with other. Essentially, the developed mechanical model of elastically interacting particles and the model of a spin system in a classical approximation are described by equivalent equations [33]. So, the results of the mathematical treatment of the equation for the evolution of a spin system are completely applicable for our mechanical model.

Equation (10) is formulated for the evolution of particles. Therefore, it is important to write down the correspondence between the former and subsequent states of a particle on a lattice in an explicit form:

Equation (10) contains two items corresponding to vector products, which describe the rotation of particles in an external field and in local interactions with neighbouring particles, respectively. The description of rotations is convenient to make by means of quaternions, that is, four numbers , where the index 0 relates to scalar and indices 1, 2, and 3 correspond to vector parts of a quaternion. The turn of vector in 3D space till the matching vector is given by the unitary transformation in the space , where a quaternion of rotation is described by the Cayley-Klein parameters of rotation:and the adjunct quaternion iswhere is the angle between vectors and () and the symbol shows an operation of multiplication of quaternions. It is worthwhile to compare the item and corresponding components of a quaternion of rotations for different angles between vectors in the plate . One can confirm the identity of these expressions if the condition is fulfilled. The physically meaningful range of angles (between the directions of orientation of neighbouring structural elements and the action of shearing) lies between 0 and . Consequently, items entering (10) at the limit of low excitations () correspond to the following quaternions of rotation:

Here, is the angle between the direction of the principle axis of an ellipsoid and the direction of shearing and is the angle between principle axes of neighboring ellipsoids.

Under the constraint of rotation angles and small amplitudes of excitation, it is possible to pass to the isomorphous quaternion presentation of rotation:where are vectors forming the 3D space and are quaternions in the 4D space .

Then, (10) can be presented in the following finite-difference form:

Analysis shows that th state of a particle can be obtained from ()th state by a combination of its consecutive three-dimensional rotations. Such local transformation of the vector space determined by the operator of the time shift is a nonlinear cubic mapping of the domain of previous states into the domain of current states. It can be proven by direct substitution of quaternions of rotation into (19).

Continuity and isotropism of turns for principle axes of structural elements of the model can be provided by the introduction of a special operator which describes the particle turns not for the whole angle but for its (small) part of the angle between the vectors. So, a new additional parameter, susceptibility, is introduced. It characterized an orientability of particles in the external and local fields. Quaternion of rotation, , providing the turn of the vector at th part of the angle along the direction of the field is given by the expression following directly from the solution of the power quaternion equation : where is a normalized vector and is a scalar part of the quaternion of rotation at an angle .

Then, bearing in mind the introduced parameter of susceptibility, the first member in (19) can be written as where is the operator of rotation for an anisodiametric particle in shearing.

The quaternion of rotation for a particle under the local field is determined in a similar manner for the angle equal to th of the complete angle of rotation :

Orientation interaction of particles should not depend on the order of arguments in noncommutative operator of rotation. Consequently, an operator describing rotation in the local field should be symmetrical. The simplest symmetrical operator of rotation is the linear combination of the two nonsymmetrical operators in which the order of arguments is changed:

It is also necessary to take into account the possibility of an arbitrary structural element in additional contact with elements in the upper and lower layers () as in Figure 4. Then, the equation for the evolution of a system of structural elements in the phase space of vectors is written as

The continuity and isotropy of the principle axes turns for structural elements are introduced by a special operator which describes the turn not only by the whole angle but also by a fraction of the angle between vectors. Thereby, an additional operator should be considered. This is susceptibility characterizing the ability of a particle to orient in either an external or in a local field.

A system under investigation is conservative. Then, the phase volume during the evolution should be constant and phase trajectories must lie inside a limited volume. The value is natural to take as a measure of the phase volume. Taking into account this normalization, (11) is reformulated as

Equation (23) models the behaviour of an elastic liquid being transferred into a rubbery state (at high deformation rates) that is characterized by the following features: rotary degrees of freedom can be realized, orientation is possible, inherent friction (energy dissipation) is excluded, and movement of particles along the lattice is impossible.

This model operates with four parameters: vector of the shear field , susceptibility to the action of shear stresses (where is a relative angle at the iteration step), the value of a local field , and the possibility of rotating in the field of neighbouring particles, .

The solution of this equation giving the deformation (or stress) and particle orientation distributions can be treated as the elastic and deformation behaviour of structural elements in a material.

#### 3. Discussion

The goal of this study consists not only in the investigation of the general and universal properties of (10) but also in attempts to find qualitative answers to the problem of elastic instability. First, is self-arrangement in a system of elastic structure elements possible? If the answer is positive, then, second, what is the utmost geometry of phase trajectories at a long time limit? Third, what is the characteristic size of the ordered domains? Fourth, where are these zones of instability in the phase space and what is their nature?

We will try to answer these questions based on simple model cases.

##### 3.1. Self-Arrangement of the Structure in the Absence of an External Shear Field

Equation (12) gives opportunities to model two limiting cases. The first supposes the absence of interparticle interactions, corresponding to . This case presents the situation of the behaviour of structural elements under the action of an external field only. This limiting case immediately leads to the problems in the continuum theory of elasticity. The second limiting situation presents the model of a system of particles behaving under orientation interactions (). In this case, self-arrangement with formation of regular structures becomes possible.

Let us discuss the latter case () as the most interesting for the goals of this study. The analysis made below will be related to the middle layer of a three-layer lattice.

The state of a system in this case is characterized by the surface enveloping the system of vectors called the Hashimoto surface [34]. The initial state of the structural elements on the lattice is assumed to be random:where is function generation random values of the vector length corresponding to the law of normal distribution with average value of and the half-width of . The function allows us to generate axial projections of vectors uniformly distributed from −1 to 1 and accordingly the uniform distribution of the vector orientation in space is given. A function norm realizes the quaternion normalization. Vector corresponds to every lattice node.

The initial state of the structural elements on the lattice is assumed to be random. According to the results of computer modeling, an increase in the number of iterations leads to an attractor, that is, some utmost configuration. An example of the phase state transformation from the initial state (a) to such an attractor (b) for a model case with parameters and is shown in Figure 5. The attractor under discussion is a pseudospherical surface with self-intersections. The lines of self-intersection are zones of bifurcation, where phase streams separate. Visual external boundaries of self-intersection are shown in Figure 5 by arrows.

**(a)**

**(b)**

The existence of bifurcations shown in the figure means that phase trajectories can change the direction of their movement and begin to move inside the different regions of an attractor under the action of force fields. Due to the existence of self-intersections, the attractor has an inner cell-like structure. The inner part of the cell is a simply connected surface; that is, the phase stream can be divided and localized inside a cell due to bifurcations. So, bifurcations create the necessary prerequisites for elastic instability of flow.

Calculations show that trajectories excited by small deformations are drawn together to the utmost surface. So, the process of self-arrangement, consisting of the formation of domains of cooperative behaviour of structural elements, takes place in a system of interrelated (interacting) rotators. It appears that peripheral regions in the phase space corresponding to strong excitations of particles (i.e., structural elements with maximum eccentricity) are the most regularly arranged. The central regions corresponding to slight excitations are less regular and the centre remains unorganized and even chaotic. A comparison of the initial and final distribution functions for is shown in Figure 6. It is clear that the initial distribution splits and two groups of structural elements appear: one is well ordered and the second remains disorganized. So, the self-arrangement of the structure accompanies the chaos-to-order transition.

Estimation of the space size of the cooperative behaviour of a system is based on the determination of coordinates, at which projections of pair correlation functions of on the coordinate axes reach zero (Figure 7).

Averaging the function for both directions (these functions are sometimes different for and directions) gives the maximum size of the ordered domains, which is equal to approximately 20 × 20 particles for the chosen parameters of the model (, ).

Separation in the cooperative behaviour of particle domains actually means that space localization of the movement modes for a system of particles takes place. It can be noted that this effect is a direct consequence of the geometry of an attractor, because simply connected domains appear as a result of its self-intersections.

##### 3.2. Instability and Structure Formation in Axial-Symmetrical Shear Fields

Let us consider the effect of mechanical shear fields with axial symmetry on the process of structural element ordering according to the model under discussion. This type of shearing takes place with rotational flows in standard measuring units, such as plate-and-plate, and sphere-and-plate rheometers, where torsion shear fields are realized.

As before, it is assumed that a local field is absent (). An example of the calculations will be discussed for the following values of the characteristic parameters: , .

The shear field in this case is given in the following manner:

This is an axial-symmetrical field with a pole in the centre of a lattice ; and are coordinates of the nodes of a lattice. The item 0.1 is added to exclude divergence in a pole.

Calculations show that after a few iterations the chaos in orientation in the lattice plane disappears and structural elements become oriented along the concentric force lines in the lattice plane (Figure 8). No anomalies in the orientation and deformation behaviour of a system of noninteracting particles are observed.

**(a)**

**(b)**

Now let us add the radial component to the axial-symmetrical field in shearing:

This case is interesting because in compression of a medium (e.g., when lower and upper parts of the working unit are converging at positioning this unit into the working position) the radial components of the flow appear.

Computer simulation shows that the particles orient along the spiral force lines in shearing and the limiting status reached after several dozen iterations (Figure 8(b)). Unlike the previous case the symmetry is broken. An asymmetric deformation field (part of particles at the lattice is compressed, light domains, while the other part is stretched, dark domains) is shown in this figure. It means that prerequisites for cavitations or spurt appear in the vicinity of the pole of the sphere-plate pair in the string shear fields with radial component when stresses exceed the cohesive or adhesive strength.

##### 3.3. Superposition of Shearing and Interaction between Structural Elements (General Case)

Superposition of shearing and local interactions creates spiral structures. The origin and evolution of such structures are shown in Figure 9 for a single-layer model with the shear field given by (24) and the following parameters:

**(a)**

**(b)**

**(c)**

As seen, in this case (as opposed to the situations discussed above) is not equal to zero.

Figure 9 shows that spiral waves of compression and extension stresses scattering to the periphery are generated near a centre. A triplet of alternating extension-compression domains appears. The centre of the spiral structure has a singularity in the form of two extremes corresponding to the maxima of the extension and compression work. This situation is similar to the action of shearing with a radial component. So, zones with negative deformations are conjugated with the appearance of a radial component in a shear field. This component is the result of oriented self-congruent interactions of particles. As the process of self-arrangement is associated with bifurcations, the resulting shear field becomes unstable. In this situation, one can predict that structural elements will be oriented in an average spiral force field with superimposed fluctuations of the field characteristics, which are commensurate with the size of the domains containing ordered particles. Results of computer modeling confirm this supposition.

##### 3.4. Some Experimental Observations

In this section, we will present some experimental observations which illustrate and confirm the results of computer simulation based on the proposed model. These data were obtained in our previous experimental studies devoted to the behaviour of polymer melts, polymer blends, and slightly filled polymer melts in strong shear flows created in deformation between cone and plate or spherical surface and plate [35, 36]. In the first case, we were dealing with homogeneous shear flow and in the second one the shear rate varied along the radius.

In light of the above described computer simulation, regular structure formation in the strong flow of a single-phase polymer melt is the most interesting phenomenon. Indeed, structures of the described type have actually been observed in flows of pure polymer melts above some critical shear rate (Figure 10(a)) and in the flow of melts loaded with hard nondeformed particles (Figure 10(b)) [36]. In both cases one can see typical circles and/or spiral structures.

**(a)**

**(b)**

The experiment carried out in the flow between spherical and plane surfaces is especially interesting. The dotted line in Figure 10(b) shows the boundary of the formation of regular structures. Estimation shows that this effect becomes possible at the Weissenberg Number of the order of 5-6 that corresponds to the well-known condition of the flow-to-rubbery transition.

The effect of particle ordering in shear flows was previously observed with planar shearing (e.g., [37]) and for flows of micellar colloids in a rotational instrument [38]. Our data presented in Figure 11 for rotational flow demonstrate the other type of self-organization related to elastic instability in strong flows.

The rings formed due to high deformation rate are positioned in a quite regular manner. The exact rule of positioning depends on details of the shearing device geometry. These regular ring structures appear as the final stage in spiral transformation.

It is important to remember that the initial distribution of solid particles in the melt was quite homogeneous. The self-assembly of particles into regular rings is the result of shearing. Figure 11 shows the frames of a movie demonstrating the gradual transformation of the homogeneous distribution into regular rings (from left to right and from upper to lower frames). This is similar to Figure 11, right, but was made for a different object and the photos were generated using polarized light.

Figure 12 presents the* qualitative* graphic illustration of the distribution of the filler particles (used as a marker) along the boundaries of “elastic tubes,” the triplet of spiral sleeves, based on the results of computer simulation (Figure 8) and compared with the experiment. It is reasonable to suppose that the experimentally observed spiral or ring distribution of particles is the result of normal stresses in nonhomogeneous elastic medium. The particles are localized in the zone of minimal stresses and play the role of a kind of indicator of the normal stress profile.

A similar morphological picture was directly observed in strong flows of a two-phase polymer blend [35]. Component separation took place and a low-viscosity LC-polymer created the structure in the form of a regular circular torus. There are two effects: first, the stretching of large separated drops and, second, migration of small drops in radial direction initially distributed homogeneously in the polymer matrix, feeding the formed circular morphology.

The consecutive stages of the structuring and its transformation by twinning are also shown in Figure 13.

**(a)**

**(b)**

Figure 13 demonstrates the periodic structure which appeared as a result of the spurt effect in a strong flow of poly(isobutylene). The twinning phenomenon typical of unstable situations is demonstrated in frame (a) where arrows show the start of bifurcations.

The structural periodicity directly corresponds to the periodic stress distribution shown in Figure 14, where Figure 14(a) is the interference image and Figure 14(b) demonstrates the structure seen under transmitting light. The light zones correspond to normal compressing stresses and the dark zones to stretching ones that qualitatively corresponds to the results of calculation presented in Figure 8(b).

**(a)**

**(b)**

The elastic instability in strong flows leads to the radical change in the mechanism of deformation. This might be a quasisolid rotation of rather large structural elements (“grains”) of a medium, which were observed for concentrated emulsions [38, 39]. We also think that “grain” as a structural element can be rather huge in order for Brownian movement to be neglected.

The issue concerning the nature of “particles” used in the above-developed model or “grains” is still open for discussion. In the frame of this study, this is only a convenient phenomenological model, though they can have a physical sense.

In conclusion, it is worth mentioning that instability in the flow of elastic liquids visually is similar to inertial instability at a rotating surface (in a boundary layer). To compare pictures presenting the wave formation in a boundary layer of a viscous liquid on a rotating disk [40] and photos of circular structures shown in Figures 11 and 13, it would be rather difficult to distinguish the difference, though the mechanism of each phenomenon is quite different. Moreover, the mechanism of bifurcation shown in Figure 13 as doubling of spiral is of general meaning for any liquid leading to turbulence of either an inertial or an elastic type (as studied by us). Possibly, there is the classical order-to-chaos transition.

#### 4. Conclusion

Analysis of literature and our own experimental data show that high deformation rates in polymer systems lead to a flow-to-rubber-like transition in their behaviour. This phenomenon consists of a pronounced dominance of elastic deformation over irreversible flow, either in shear or in extension. It is the physical limit for flow of elastic liquids. As a result, instability and different effects demonstrating the formation of ordered structures were observed. A proposed concept describing behaviour of this type is based on the refusal from the continuum field approach and developing a discrete phenomenological model that consists of elastic rubber-like rotators. This is a behaviouristic model of elastic liquid behaviour in strong stress fields. The movement of the structural elements in this model is determined by their elastic deformations in shearing (transition from a spherical to ellipsoidal shape of particles and their orientation) and by elastic interactions of these elements due to their collisions. The computer simulation demonstrated that this model really describes the consequence of elastic instability as the chaos-to-order transition and creation of experimentally observed regular structures.

#### Appendix

#### Derivation

The equation for interaction in the system of elastic ellipsoidal particles is obtained via the Hamilton function.

At first, let us calculate the energy of pair interaction for elastic ellipsoids located at the sites of the uniform coordinate lattice. It is assumed that the centre of mass of ellipsoids is located at the lattice nodes and particles have only rotation degrees of freedom. The ellipsoid space location is determined by the vector coinciding with the direction of the principle axes and the shape is characterized by eccentricity (modulus of the vector ). Particles in their contact can expose large deformations and therefore the Mooney-Rivlin potential with constant and was used.

Computer simulation of deformations of particles contacting with each other was carried out in the* Comsol Multiphysics* media by the finite-difference method. Figure 15 demonstrates two consecutive stages of deformation accompanied by changes in the relative position of ellipsoidal particles. The distribution of the elastic energy density is shown by colors.

**(a)**

**(b)**

Figure 16 presents the results of calculations for the energy of space interaction of particles, , obtained by the integration of the elastic energy density throughout the whole particle volume and subsequent averaging for a unit particle based on the Mooney-Rivlin and (for comparison) Hooke models.

As is seen, the calculated dependence in Figure 16 is linear for both types of the elastic potential. So, the simulation results can be approximated by the following equation:where is the elastic modulus of particles (for the Mooney-Rivlin model, ), is the scalar product of vectors and , and is a coefficient of the order of unity.

The second basic equation used in the model under discussion is the elastic energy density for an ellipsoidal particle under external shearing field. For this, an ellipsoidal region was marked out inside a cube. In the initial state, the major axle was oriented perpendicular to the shear stress, , applied to the cube surface. Figure 17 illustrates some shots of the shearing process. During shearing, the ellipsoidal region is deformed and rotates at the shear direction. Then the averaged value of the elastic energy density localized in the ellipsoid volume, , was calculated. Figure 18 presents the dependence of normalized by the shear stress on the angle between the shear direction and the major axle of ellipsoid, . This dependence appears linear if we use the value as the argument.

**(a)**

**(b)**

So, computer simulation showed that the elastic energy density of an ellipsoidal particle under shearing is expressed as where is the unit vector of the shearing field, is the scalar product of the vectors and , and is a coefficient of the order of unity.

It is worth mentioning that the expression (A.2) for is written not as the standard invariant squared form of the components of the stress tensor (as is usually done, e.g., for shearing , where is the shear modulus) but through the fitting function including the shear stress and the scalar product of the vectors and ( is the vector characterized by the direction of the principle axis and the ellipsoid eccentricity). It also follows the invariance of the obtained expression for in relation to the coordinate transformation (their rotation and transposition).

Such an equivalent presentation of the elastic energy in shearing allows for simulation of the behavior of an elastic particle under shear without resorting to the calculation of the stress tensor components. This is especially important in connection with the transition from continual to the discrete method for description of a grainy elastic medium.

Equations (A.1) and (A.2) are the base for building the model under discussion in whole.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The financial support for this work from the Russian Science Foundation (Agreement no. 14-23-00003 of August 2014) is gratefully acknowledged.