Mechanics of Axially Moving Structures at Mixed Eulerian-Lagrangian Description

We discuss a series of methods of the mathematical modelling of large deformations of axially moving strings, beams and plates. Both uni-axial and looped trajectories of motion are considered, which allows the application of these methods to such practically important problems as rolling mills or belt drives. Based on the principles of Lagrangian mechanics, we transform the variational formulations of structural mechanics to problem-oriented exact kinematic descriptions with mixed spatial and material coordinates. The discretization of an intermediate domain results in a consistent non-material ﬁnite element formulation with particles of continuum ﬂowing across the mesh. This allows avoiding numerically induced oscillations in the solution, while keeping the discretization ﬁne where necessary, e.g. in the regions of contact.


Introduction
Axially moving structures are frequently used in engineering as an intrinsic part of many technical solutions. It is not surprising that they have long attracted the attention of mechanical engineers. Practical relevance, technical difficulties in maintaining the desired regime of motion, non-trivial and even sometimes counter-intuitive behavior are coupled with challenges, intrinsic to the theoretical investigation of such systems.
The mathematical modeling is traditionally based on a spatial (Eulerian) kinematic description, see review papers by Chen (2005); Marynowski and Kapitaniak (2014). Considering the unknown displacements, forces and other mechanical fields as functions of a fixed coordinate in the axial direction simplifies the analysis, as the problem needs to be solved in a fixed control domain only and the boundary Yury Vetyukov Institute of Mechanics and Mechatronics, Vienna University of Technology, Getreidemarkt 9, 1060 Vienna, Austria e-mail: yury.vetyukov@tuwien.ac.at conditions are applied at fixed points. The common difficulty is that the basic equations of structural mechanics are generally available in the Lagrangian form, where the variables are observed in material points. The present contribution is focused on consistently solving this issue by transforming the equations to new problem-specific coordinates.
Below, in Sect. 12.2, we demonstrate the transformation of the equations to the spatial form in a simple linear case, when the relation between the material and the spatial coordinates is known in advance. In nonlinear problems, many authors adopt certain kinematical simplifications and derive the mathematical model in the spatial frame from scratch (Ghayesh et al, 2013;Kong and Parker, 2005), or rely upon other known equations from the literature. Thus, the models of Wickert (1992); Mote (1966) are often applied. While numerous fascinating and practically important results were obtained in this way, accurate modelling of large deformations of axially moving structures requires exact treatment of geometrically nonlinear effects both in the elastic response as well as in the inertial properties with an established structural mechanics theory in the background. Despite growing interest of researchers (Humer, 2013;Humer andIrschik, 2009, 2011;Eliseev and Vetyukov, 2012;Pechstein and Gerstmayr, 2013;Vetyukov et al, 2017b), developing efficient and reliable techniques of transforming the general equations of motion to a new spatial form remains a challenging problem.
In statics of flexible solids, modelling of a discretized system is convenient with the principle of minimality of the potential energy. For a given finite element (or global Ritz) approximation one simply integrates the total strain energy and energy of external forces as functions of generalized coordinates, and the derivatives of these functions constitute the equations. An extension to dynamics is traditionally accomplished with Lagrange's equations of motion of the second kind; the kinetic energy becomes a quadratic form of generalized velocities. The approach is particularly straightforward when combined with the material description of the kinematics of deformation and motion. Considering mechanical fields as functions of coordinates in the reference configuration is advantageous because an elastic body keeps "memory" of its undeformed state. Moreover, dealing with the same material volume ensures the validity of the variational principle and the corresponding equations. Modern applications pose challenging problems of computational modelling of fluid-structure interaction, material forming processes, etc. Similar to fluid mechanics it becomes more efficient to observe the processes at given points in space, which gave birth to the so-called Arbitrary Lagrangian-Eulerian (ALE) formulations. This family of methods features control volumes, which are moving in a problem-oriented manner relative to both the spatial actual state as well as the reference configuration (Donea et al, 2004). In their traditional form, ALE methods imply accounting for the advection of material in the local forms of the constitutive and balance equations. Although their applications in structural mechanics are presently sparse, we can mention the works by Hong and Ren (2011);Yang et al (in press). The authors of the mentioned papers make use of a redundant set of degrees of freedom in a finite element model with additional constraints. This flexible formulation results in a differential-algebraic system of equations, which needs to be integrated over time.
In this chapter, we systematically discuss the application of a novel kinematic description of thin axially moving structures, which undergo finite deformations superimposed upon the gross nominal motion. The main results were partially presented by the author previously in Vetyukov (2017); Vetyukov et al (2016). The presented work in progress is currently actively developed towards the consistent model of dry friction with stick and slip phenomena in the zones of contact of the axially moving continua with pulleys and drums.

Linear Waves in a Moving String
We begin with the simple example of waves, running in a moving string. A particle of a string is identified by a material (Lagrangian) coordinate s. We assume small transverse deflection w and constant tension force and write for a dynamic process over time t with the wave velocity c. The mechanical field w is considered at a given material point s, and the wave equation features a material time derivative ∂ t , which is typical for the Lagrangian formulation. In mathematical physics, one conventionally considers (12.1) with boundary conditions. Our string is moving axially between two spatially fixed points, and the boundary conditions need to be posed in time varying material points s 1,2 (t). Moving boundary conditions make the problem essentially more difficult and challenging for solving by either the established methods of mathematical physics or numerical techniques of finite differences or finite elements. Assuming a constant velocity v of the gross axial motion, we find simple relations between the material coordinate s and the spatial one x, along which the string is moving, see Fig. 12.1: In the following, we will mainly deal with the Eulerian form of equations, in which the unknown fields are observed at a given spatial coordinate x. Derivatives with respect to time and space need to be transformed from the material to the spatial form: We introduce a shortened notation: a prime means a derivative with respect to x and a dot denotes a time derivative at given x. Applying Eqs. (12.3) twice, we transform the wave equation (12.1) to the Eulerian form and arrive at an equivalent boundary value problem with the new pair of variables x and t: The boundary conditions are imposed at fixed points x = 0 and x = l, the latter being the length of the control domain. The equation appears to be more complicated, but the boundary value problem may be solved either analytically in the form of infinite series (van Horssen and Ponomareva, 2005), or, even easier, numerically using the method of finite differences. The example solution for the following values of parameters and the smooth initial state c = 1, l = 1, v = 0.5; t = 0 :ẇ = 0, w = w 0 (1 + cos(6π(x/l − 1/2)))/2, l/3 < x < 2l/3 0, x < l/3 or x > 2l/3 is shown for three instances of time in Fig. 12.1; SI system of units will be used throughout the Chapter. The chosen large sample value of the amplitude w 0 lies certainly far beyond the range of applicability of the geometrically linear wave equation (12.1). The string is moving, and the sequential positions of a material point, which is originally in the middle of the string, are marked in Fig. 12.1 by a cross. We see two waves of unequal intensities, travelling with the velocities c + v and c − v respectively to the right and to the left and reflecting from the end points. From (12.4) we conclude that the type of the partial differential equation changes when v grows above c, which means loss of stiffness and dynamic instability.

Large Vibrations of an Axially Moving String
The mathematical model for large vibrations of an axially moving string or a beam was presented by the author in Vetyukov (2017). The key points of the analysis are discussed below, we restrict the presentation to the string case without bending stiffness. We consider the planar motion of a string across a given domain, see Fig. 12.2. The material enters the domain at the left end x = 0 and leaves it at the right one x = l with the same rate γ. This material length of the string per time unit is essentially different from the spatial velocity of the particles in the axial direction. Such kind of kinematically prescribed influx or outflux of the material may be implemented using a timing belt and a toothed pulley. It means that the total material length of the string, which remains in the control volume of the considered problem, remains constant.
The kinematically prescribed time-varying transverse deflection of the end point brings the system into a complicated in-plane motion with coupled axial and transverse dynamics whenẏ l is not small, see Fig. 12.3 for a plot of (12.6).
Although the equations of string dynamics at finite deformations are well established at Lagrangian description (Eliseev, 2006;Eliseev and Vetyukov, 2012), it is a numerically challenging task to solve them with the boundary conditions moving across the material domain. Thus, an attempt to impose kinematic constraints at  Kinematically prescribed transverse coordinate of the right end of the moving string moving points of a finite element mesh leads inevitably to numerically induced oscillations in the solution. Therefore, in the following we will use the principle of virtual work as a starting point and transform the equations to a more advantageous problem-specific kinematic description.
For the material description, we consider the displacement from an initially straight horizontal reference configuration being a function of the material coordinate s and time t, such that the actual position vector of a particle is r = xi + y j = si + u, x = s + u x , y = u y ; (12.8) the unit vectors in the directions of the coordinate axes are respectively called i and j.
The longitudinal strain measure reads The strain energy U per unit material length s is a quadratic function of the strain, (12.10) in which b is the tension stiffness.

Mixed Eulerian-Lagrangian Description of the Kinematics of Motion
Further exploiting the simple idea of the mixed Eulerian-Lagrangian description, we change the variable. Instead of the unknown fields u x (s, t), u y (s, t), we seek the material coordinate and the transverse displacement as functions of the axial spatial coordinate and time: (12.11) The presentation (12.11) is incapable of describing the formation of loops, at which several material points have the same coordinate x. However, we gain more efficiency in the analysis of practically relevant problems with large vibration amplitudes. Again introducing short notation for the derivatives with the pair of variables x, t similar to (12.3) and using we rewrite the strain measure (12.9) to the spatial form Now we address the velocity of a particle and compute it in the mixed description, in which a convective term comes into play: (12.14) The acceleration results into

Mixed Form of the Variational Equation of Virtual Work
The boundaries of the active material segment of the string, which is currently in the control domain, vary over time with the known rate: We restrict the presentation to the simple case of known rates of the material influx and outflux and begin with the material form of the variational equation of virtual work for the active segment of the string: With • δ we denote a material variation, i.e. a variation of a mechanical field in a given material point. Owing to the kinematic nature of the boundary conditions, we state that the variation • δu vanishes at the current boundaries, and there is no virtual work produced by the interaction forces with the parts of the continuum outside the active domain.
Using the mathematics similar to (12.14), we relate the material variation • δ at fixed material point and the spatial one δ at given x: ϕ is any mechanical field, which can be considered as a function of both, the material and the spatial coordinates. Transforming the elastic term in (12.17) to the spatial form is simply: The distributed strain energy U remains a quadratic form (12.10) of the strain measures (12.13). The inertial term in (12.17) reads (12.20) besides (12.12) we also used the transformation rule for the variations (12.18). Substituting (12.19) and (12.20) in (12.17), we obtain a non-material variational equation The latter can be used in a numerical procedure with a suitable Ritz-Galerkin approximation of the unknown displacement u(x, t). More elegance, simplicity in implementation and efficiency can be reached by transforming the equation to the form of Lagrangian equations of motion of the second kind for a given approximation of the unknown displacements.

Finite Element Approximation
The finite element approximation of the unknown field of displacements is set in the form Both cartesian components are independently represented by piecewise-cubic shape functions u k (x), which allows achieving C 1 inter-element continuity with the corresponding fine rate of mesh convergence. With equal finite element sizes and a local coordinate ξ on each finite element, which varies from −1 to 1, we have a linear mapping from ξ to x and use four cubic shape functions u el i (ξ) (Fig. 12.4) with the properties Fig. 12.4 Cubic shape functions, which provide C 1 inter-element continuity (12.23) the four conditions uniquely determine four coefficients of cubic polynomials, which often find use in the models of classical rods and shells (Vetyukov, 2014b(Vetyukov, , 2012. The global shape functions u k (x) are constructed from the local ones u el i (ξ), degrees of freedom q k,α determine either the value or the derivative of the corresponding component of u in a node. Each node has now four degrees of freedom, and n is twice the number of nodes in the finite element model. The boundary conditions in the first and the last nodes, kinematically prescribe some of the functions q k,α , which thus deleted from the set of active degrees of freedom in the model.

Lagrange's Equation of Motion of the Second Kind
The dynamic part of Eq. (12.21) can be represented in terms of the kinetic energy of the control domain In Vetyukov (2017) it was rigorously proved, that the Lagrange's equation of motion of the second kind d dt retains validity in the considered case of non-material control volume, provided that the shape functions of the active degrees of freedom vanish at both ends: (12.27) which is the case for kinematic boundary conditions at hand.

Example Solution
The finite element scheme including the time integration was implemented in the Wolfram Mathematica 1 environment. With n el being the number of elements in the model, we have n el + 1 nodes with 4(n el + 1) degrees of freedom q k,α , 4 of which are prescribed by the boundary conditions and excluded from the active set. The total strain energy (12.19) and kinetic energy (12.25) were computed using a Gaussian quadrature rule with 3 integration points per element. Substituting further the known time variations for the kinematically prescribed degrees of freedom, we construct the system of second order differential equations for the transient dynamics in the form (12.26). With the values of the numerical parameters, initial values and time derivatives of degrees of freedom we integrate the equations over time using standard NDSolve routine. Setting the option Method -> {"EquationSimplification" -> "Residual"} allows Mathematica to treat the problem as a system of differentialalgebraic equations and to use the corresponding solver, which is more efficient and robust in the present case. The obtained results are ready for post-processing. At first, we considered the linear problem of wave propagation in a moving string, discussed in the motivational example, Sect. 12.2. Setting the numerical values and the initial conditions in accordance with (12.5), we were able to reproduce the results of Fig. 12.1 with a high level of accuracy even using small n el .
The numerical experiments presented below feature the following parameters of the model with large deformations: (12.28) Scalability of the model allows testing the numerical scheme for unit tension stiffness and inertia. The structure is pre-stressed, and the initial tension strain determines the material length of the active part of the string s l − s 0 , whose length does not change during the simulation. In the beginning, the structure is moving steadily with The initial velocity of transverse waves is determined by the pre-tension and for the given numerical parameters can be estimated as √ 0.1 (square root of the tension divided by the linear density), which is comparable to the rate of material supply γ. The transverse motion of the right end y l (t) according to (12.6) increases the general tension level and the wave velocity over the time period 0 ≤ t ≤ t * .
Taking n el = 10, we sequentially solved the transient problem for both axial transport rates γ, given in (12.28); a simulation takes about 10-15 seconds of CPU time with a modern desktop computer with an Intel i7 CPU 2.3 GHz. The computed deformed configurations of the structure for the time instant t = 10 are presented in Fig. 12.5. The curves here are simple plots of u y in dependence on the axial coordinate x, the domains of the finite elements are indicated by gray vertical lines, and the nodes of the finite element model reside on this lines. More informative are the time histories of characteristic variables, which are shown for both solutions in Fig. 12.6. In the middle point of the domain x = l/2 we observe the time variation of the transverse displacement u y , of the strain ε and of the axial component of the velocity of particles v x ; we also present the time variation of the strain ε in the end point of the domain x = l. The starting values of the velocity v x are not equal to γ because of the pre-strain. We notice that the strain at the right end immediately responds to the growing y l (t), but it takes some time before the influence in the middle point is observed.
For the demonstration of the rapid mesh convergence, extension to the case of a beam with the bending stiffness and modelling dynamic instability at growing velocity of axial motion we refer to Vetyukov (2017). The formulation was also validated by successfully comparing the results in Fig. 12.6 against a lumped particles simulations at Lagrangian description, which featured thousands of particles in contrast to just a few number of degrees of freedom in the presented model, see the above reference for the details of comparison.

Finite Deformations of an Axially Moving Plate
We proceed to modelling finite deformations and bending of an elastic plate moving across a given domain and present the results, which were presented in Vetyukov et al (2016) under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/). Previous attempts to obtain a model for an axially moving deformable plate (Ghayesh et al, 2013;Banichuk et al, 2010) did not account for particular kinematic conditions at the boundaries of the domain in the form of time and space variations of the prescribed velocity profiles. These effects are particularly relevant in applications of hot rolling, strip coiling or belt drive simulations.
In this section we discuss a mathematical description for finite deformations of a thin plate, which is moving between two parallel lines, see Fig. 12.7. This gross axial motion takes place in the direction of the x-axis (called in the following axial direction) from left to right. Velocities of the plate at the entry to the domain x = 0 and at the exit from it x = L are kinematically prescribed and equal respectively v entry and v exit . In a transient process, a non-constant profile of the velocity v exit (y) results in the deformation of the plate in its plane xy and subsequent buckling in the out-of-plane direction of the z-axis. This kind of undesired behavior may happen during the metal rolling process, paper production, etc., which makes the respective methods of analysis practically relevant. The proposed Eulerian-Lagrangian formulation is motivated by computational difficulties that are inherent in conventional finite element modeling of the deformation of the plate. As a Lagrangian finite element mesh propagates across the domain in the course of a simulation, kinematic boundary conditions need to be imposed at the interior of the finite elements. This typical "variational crime" leads to numerically induced oscillatory behavior of the stressed state of the plate and this would not allow for accurate presentation of its bending, which is tightly coupled with the in-plane stresses. While extended finite element formulations (Moës et al, 1999) or a mortar approach (De Lorenzis et al, 2011) could theoretically be applied to solving the issue, an elegant and efficient alternative is to preserve the finite element mesh from moving in the axial direction. This guarantees that the boundaries of the active domain, at which the displacements need to be prescribed for each time step, remain aligned with the boundaries of the finite element mesh.
At present, we restrict the consideration to the constant profile of the velocity of material production v entry and aim at modeling the evolution of finite deformations of an elastic plate under quasistatic assumptions. Further extensions of the approach are discussed in the concluding remarks below. Further extensions of the approach are discussed below in Sect. 12.5.

Kinematic Description
Consider a plate as a two-dimensional material surface axially moving across a domain, bounded by two lines, see Fig. 12.8. Rolling of a sheet of metal, paper production or motion of a conveyor belt are typical examples of such a structure. At each time instant, the plate is clamped at the lines of contact: there are no out-of-plane displacements in the direction of the z axis at x = 0 and x = L nor does the plate rotate about these lines. The material is produced in a fixed segment −w/2 ≤ y ≤ w/2 at x = 0, such that the plate is here centered and has a fixed width w. The second line of contact, where the plate is leaving the domain, may travel in the transverse direction, see discussion in Sect. 12.4.5.
For the sake of mathematical description, we introduce an infinitely long reference configuration, see Fig. 12.9. Particles of the plate are identified by their material coordinates • x and • y, i.e., by the reference position vector with i, j and k being the unit base vectors of the global cartesian frame. In the actual state we consider just the particles, which are currently residing in the domain 0 ≤ x ≤ L. In the reference configuration these particles form a planar manifold, which will be called active material volume in the following (depicted grey in Fig. 12.9), and which is thus a pre-image of the actual state. Between the time instants t and t + dt, a material layer of the length v entry dt enters the domain from the left. As v entry = const, it means that the left boundary of the active material volume (which is a pre-image of the line x = 0, −w/2 ≤ y ≤ w/2 in the actual configuration) is moving across the stress free reference configuration to the left with the velocity v entry , the mathematical form of this statement being (12.71). The Lagrangian description, at which the actual configuration is defined by a mapping from • r to the actual place of a particle r = xi + y j + zk, (12.32) remains efficient as long as the active material volume is known (and suitable for being discretized in numerical simulations). This is not the case for the problem at hand, and it is efficient to decompose the mapping into two steps by introducing an =xi +ỹ j, 0 ≤x ≤ L, −w/2 ≤ỹ ≤ w/2 (12.33) differs from the reference state by the axial displacement u x only: At the same time, the transverse displacement u y and the out-of-plane displacement u z differ in the actual configuration from the intermediate one:

Deformation and Strain Energy
Treating displacements u = u x i + u y j + u z k (12.38) as primary variables, we seek them as functions of the place in the intermediate configuration: The quasistatic simulations below are based on the principle of stationarity of the total strain energy of the active volume of the plate for each time instance. For a classical (transverse shear-rigid) Kirchoff plate model (Eliseev, 2006;Eliseev and Vetyukov, 2010;Vetyukov, 2014b), the strain energy depends on the two symmetric strain tensors E = 1 2 (F T ·F − I 2 ), K = F T ·b·F. (12.40) The first strain measure E describes the membrane (in-plane) deformation of the plate, while the second one K is responsible for its bending. The gradient of deformation is the transposed gradient of the position vector of the actual state with respect to the reference one: (12.41) The in-plane identity tensor is which is the negative gradient of the unit normal vector n to the actual surface of the plate (Stoker, 1989;Ciarlet, 2005). We recall the relation between the differential operators of the reference and of the actual states, and rewrite the second equality in (12.40): Now, the displacement (12.39) is a field over the intermediate configuration, which does not allow us to directly compute the strain tensors using (12.40). We aim at transforming the above kinematic relations from the material description featuring derivatives with respect to the coordinates in the reference state to the differential operator of the intermediate configuratioñ (12.46) For a two-stage mapping from the reference configuration to the actual one, the total deformation gradient is a product which is mathematically equivalent to the chain rule of differentiation of a function of multiple arguments. We introduce the gradient of deformation from the reference configuration to the intermediate onẽ Indeed, using a relation between the differential operators analogous to (12.44) and recalling (12.34), we write • ∇ =F T ·∇, (12.49) The inverse of an in-plane tensor in (12.48) is again an in-plane one, and the total deformation gradient F can now be computed for any field of displacements (12.39), known in the intermediate configuration (12.33). While this is sufficient for the membrane strains E in (12.40), further mathematics is needed to compute the tensor of bending strains K. Using the first relation in (12.49), we rewrite (12.45) as in which a symmetric tensorK = −∇n·∇r T (12.51) is introduced. The vector of unit normal is easy to compute with the known relation r(r) = r(xi +ỹ j). Considering a surface, parametrized by coordinates x andỹ, we write n = ∂x r × ∂ỹr |∂x r × ∂ỹr| .
(12.52) From (12.35) we see, that n at a given point of the intermediate configuration does not depend on the field u x (r). Now we apply∇ to both sides of the identitỹ ∇r· n = i ∂x r + j ∂ỹ r · n = 0 (12.53) and arrive at an alternative form of (12.51), which shows the symmetry ofK more clear:∇∇ r· n = −∇n·∇r T =K. (12.54) The above general nonlinear kinematic relations may be simplified when particular deformations or coupling terms are negligible. Thus, in a geometrically linear model just small deformations are assumed to be superposed upon a regular axial motion of the plate, and linearizing we arrive at the conventional strains of the classical plate theory: In a linear theory, one does not need to differentiate between ∇ and∇, and conventional finite element models of a plate may be applied. However, in this linear formulation the in-plane stresses are decoupled from the bending of the plate, which makes it inapplicable for studying buckling behavior.
Returning to the fully nonlinear model, we write the strain energy of a plate per unit area of its reference configuration as a quadratic form of the strains (Eliseev and Vetyukov, 2010): Thus, we assume that the local strains remain small, which does not exclude large overall deformations of a thin structure. For a homogeneous isotropic plate with the thickness h, Young's modulus E and Poisson's ratio ν, the stiffness coefficients are (12.57) Now, we need to integrate the strain energy over the intermediate configuration. The change of variables is to be accounted for, In the absence of external loading, U Σ will be at minimum in a state of stable static equilibrium provided that the set of material particles is fixed. Although we are dealing with an open system, and particles are continuously entering and leaving the active material volume, this condition is fulfilled at each time step of a quasistatic simulation owing to the kinematic nature of the boundary conditions at both ends of the domain, see the discussion in Sect. 12.4.5.

Finite Element Scheme
In a numerical scheme, we need to compute the total strain energy U Σ and its derivatives for a given approximation of displacements u(r). Similar to the above onedimensional case (12.22), we have applied a four-node C 1 -continuous finite element approximation, which is based on the idea of a Bogner-Fox-Schmit rectangle with bi-cubic shape functions, see Vetyukov (2014a,b). A regular mesh was introduced in the intermediate configuration with n x × n y elements in the axial direction and in the transverse one, respectively. With twelve degrees of freedom at each node (vector of displacement, its two derivatives with respect to the local coordinates on the element and one mixed second-order derivative) we guarantee the continuity of both u and n across the boundaries of the element. Now, the first strain measure E is continuous and the second one K may undergo discontinuities of the first kind, which is sufficient for the regularity of the integral (12.60). At each time instant we fix the values of displacements at the boundaries x = 0 and x = L according to (12.75). Along with the condition of clamping n = k, this results in kinematic constraints for some of the degrees of freedom. Seeking for a stable equilibrium, we minimize U Σ with respect to remaining degrees of freedom of the model e with the help of a quasi-Newton scheme, which requires computing derivatives of the kind we see that either the first or the second term remain in the last expression depending on the kind of degree of freedom under consideration. Indeed, Again, either the first or the second term needs to be treated depending on the kind of e. Although programming the computation according to Eqs. (12.62)-(12.65) is feasible, it would be a challenging task to obtain the second-order derivatives in an algorithmic manner. In the present study, we used a system of computer algebra and first computed the first and the second order derivatives of the strain measures with respect to the local geometric characteristics ∂xu, ∂ỹu, ∂ 2 x u, ∂ 2 y u, ∂x∂ỹu. (12.66) The resulting expressions were exported as automatically generated C# code into the in-house finite element simulation software. Now, the derivatives of (12.66) with respect to the nodal degrees of freedom are just the shape functions of a finite element, and the element's stiffness matrix and force vector are easily computed using the chain rule.

Benchmark Problem
The kinematic formulation and the resulting finite element scheme were validated on a simple benchmark problem. Consider a plate of a trapezoidal shape. The width of the plate is w, while the other two edges have lengths L and L + u x0 as in Fig. 12.10. The right (inclined) edge is clamped and kinematically rotated such, that it becomes parallel to the left one, which is also clamped. Depending on the magnitude of u x0 , the resulting deformation of the plate in its own plane may become unstable and lead to the out-of-plane buckling as depicted in the figure by the surface with the grid lines on it.
As the deformed configuration is bounded between two parallel lines, the problem is perfectly suitable for the mixed Eulerian-Lagrangian description. At the same time, a solution with the conventional Lagrangian finite elements with the same approximation (Vetyukov, 2014a) is readily available, which allows comparing the solutions and thus testing the new numerical scheme.
Using SI system of units here and in the following, we summarize the parameters of the considered model: Seeking for a static equilibrium, we considered an intermediate configuration 0 ≤ x ≤ L, −w/2 ≤˜y ≤ w/2 and prescribed u x and u y at the edge x = L such, that its pre-image in the reference state • r would correspond to the given undeformed geometry and the length of the rotated edge is preserved. The out-of-plane buckling was promoted by imposing a small gravity force in the first iteration of the quasi-Newton minimization scheme, which was then released for the rest of iterations.
In accordance with Eq. (12.35), the transverse grid lines of the finite element mesh, shown in the deformed state in Fig. 12.10 remain parallel to the y-axis. The surface of the plate in the supercritical state after buckling may be better observed from an isometric viewpoint as shown in Fig. 12.11. We studied the mesh convergence and compared the solutions with the proposed mixed Eulerian-Lagrangian description and with the conventional Lagrangian finite elements by finding the maximal and the Fig. 12.10 Undeformed geometry of a trapezoidal plate and its postbuckling configuration after the inclined edge is kinematically rotated Fig. 12.11 Deformed configuration of a trapezoidal plate, kinematically loaded in its own plane minimal values of the displacement u z , which are observed at the edgesỹ = ±w/2 of the plate. The results of the comparison, which are summarized in Table 12.1, clearly demonstrate that both formulations converge to the same solution.

Time Stepping and Boundary Conditions
Time rates of mechanical entities need to be considered in a quasistatic simulation. We need to differentiate between the velocity of a particleu with the material time derivative (· · · )˙≡ ∂(· · · ) ∂t • r=const (12.68) and the time rate ∂ t u with the local derivative at a given point in the intermediate configuration defined as The general relation between the material and the local time derivatives with the known convective term readṡ Observing the time evolution of nodal variables, we need to deal with the local time derivatives as the nodes of the finite element mesh are fixed in the intermediate configuration. As discussed after Eq. (12.31), the left boundary of the active material volume moves across the reference state with a known velocity, which means that the local time rate of u is known at the left edge: (12.71) On the opposite, at the right boundary we know the actual velocities of the particles, with which they move across the line of contact: (12.72) This explains the fact that the right line of contact will move in the transverse direction as long as ∂xu has a non-zero component in the direction of the y-axis (the z component is always zero owing to the condition of clamping). Indeed, in Fig. 12.7 the plate is slightly inclined immediately before x = L and it is easy to imagine, that in the next time instant a new cross-section of the plate with greater values of y will arrive at the line of contact. Stationary motion of the plate with is possible only when the static deformation u s fulfills (12.74) These boundary conditions allow seeking the stationary solution by solving just a single static problem. We turn back to the problem of evolution of the deformation of the plate in time. The strategy of computing the solution u (k+1) in the end of a time step t (k+1) = (k + 1)τ for a known state of the plate u (k) in the beginning of the time step t (k) = kτ comprises the following two stages.
1. First we compute new displacements at the left and right boundaries of the domain, (12.75) in which the local time derivatives are determined by (12.71) and (12.72). This sort of explicit time integration requires a moderately small time step size and has a clear advantage: with known displacements u x at x = 0, L we have a fixed material volume for the end of the time step, which allows finding the equilibrium by seeking arg min U Σ . In practice, updating the nodal variables according to the incremental boundary conditions (12.75) requires a derivative of the velocity profile v exit (y) as the nodal unknowns comprise not only displacements, but also their derivatives in the transverse direction. 2. As soon as the kinematic constraints are applied for some of the degrees of freedom of the nodes at the left and right boundaries of the domain, we proceed to the quasi-Newton iterative scheme and seek the values of the rest of degrees of freedom in the model by minimizing the total strain energy (12.60). At each iteration we compute the force vector by the derivatives ∂U Σ /∂e, and the stiffness matrix is evaluated only when the rate of convergence at the previous iteration was unsatisfactory. Converged iterations provide us with the new equilibrium state u (k+1) .

Simulation of a Moving Plate
We considered a plate with the width w, length L and material properties as in the benchmark test (12.67). As the behavior of the system depends strongly on the ratio between the membrane and the bending stiffness of the plate, we considered two values of its thickness and started with a thick plate with h = 0.02. Velocities at the boundaries were chosen as follows: v entry = 5, v exit = 5 + 1y. grows in time. The mean value of v exit increases as well, and a tension force appears in the plate. However, the part of the plate nearỹ = −w/2 is compressed in the axial direction, which (along with the shear deformation owing to the motion of the line of contact) may result in the instability of the plane configuration. We integrated over the time span 0 ≤ t ≤ 1.25 with the step size τ = 0.00125, which allowed avoiding instability in the explicit time integration of the boundary conditions for finer meshes. Experiments showed that further reduction of the time step size results in nearly identical solutions. Varying the number of finite elements in the transverse direction n y , we kept the length to width aspect ratio of the elements (L/n x )/(w/n y ) ≈ 1.2.
The computed configuration of the plate at t = 0.4 with n y = 18 is shown in Fig. 12.12. As discussed after (12.76), the linear component in v exit leads to a growing in-plane deformation, which soon initiates the out-of-plane buckling of the plate. The computed time histories of the maximal and the minimal (negative) displacements u z as well as of the mean transverse displacement of the line of contact u * y are shown in Fig. 12.13 for three levels of mesh refinement. Besides rapid mesh convergence and transition to a stationary regime, we clearly observe the time instant of buckling are produced in the simulation. A thin plate with h = 0.01 has significantly lower bending stiffness in comparison to the membrane one, which makes the simulation more challenging from the numerical point of view. The deformed surface, presented in Fig. 12.14, is less smooth than the previous one in Fig. 12.12. Lines with high curvature ("folds") may now Fig. 12.14 Thin plate at t = 0.4 computed with 27 finite elements in the transverse direction, the "folds" are observable in the deformation pattern be observed, which require dense finite element meshes to be accurately resolved. In particular, the boundary layer near x = L tends to demonstrate a complicated deformation pattern. The effect is getting more pronounced for yet thinner plates (as the one, considered in the benchmark test in Sect. 12.4.4). Nevertheless, the solution with n y = 27 finite elements in the transverse direction is quite accurate, as it may be concluded from the time histories, shown in Fig. 12.15. Earlier buckling of the thin shell may also be observed.

Mixed Eulerian-Lagrangian Formulation in the Analysis of a Belt Drive
The belt drive is a very common technical solution. Mathematical modelling of its behavior is, however, often a challenging task at Lagrangian description, as the particles of the belt are permanently moving from the free spans to the contact zones and back. Our interest is commonly focused on the time history of mechanical fields in a particular point in space rather than on following them in a given material point of the belt. As it is typical for axially moving structures, Lagrangian (material) kinematic description is not an optimal choice for this sort of problems, in particular for numerical methods with spatial discretization. Purely Eulerian (spatial) (Eliseev and Vetyukov, 2012;Vetyukov et al, 2017b) or mixed Eulerian-Lagrangian (Vetyukov et al, 2016(Vetyukov et al, , 2017a formulations are advantageous as we discretize the problem in a domain-specific manner. Exact solutions of the problem of an extensible belt moving between the pulleys with dry friction law of contact were provided by Rubin (2000); Bechtel et al (2000); Morimoto and Iizuka (2012) for stationary regimes of motion. Leamy (2005) as well as Kim et al (2011) considered perturbations of stationary motion, which should allow studying transient regimes with sufficiently small and slow deviations from a given steady one. Hong and Ren (2011) suggested a mixed finite element formulation, which combines Lagrangian nodes with Arbitrary Lagrangian-Eulerian ones using constraint conditions, typical for Multibody Systems Dynamics, and which can theoretically be applied to the considered class of problems.
In the present section, we discuss, how the previosly used mixed Eulerian-Lagrangian kinematic description can be transferred to the looped motion of a

Problem Statement
We consider a flexible belt with the tension stiffness b, which is stretched on two identical rigid pulleys with radius R and is hanging in the field of gravity g, see Fig. 12.16. With the distance between the centers of the pulleys H, we denote by λ the ratio between the shortest possible contour of the belt L and its material length L mat : Fig. 12.16: Weakly pre-tensed rubber belt, hanging on two pulleys in the field of gravity; the curvilinear coordinate system is plotted in the background (12.78) In the following, we mainly focus on the case λ > 1, such that the belt is under tension just because of the gravity force. This situation is mostly challenging for the mathematical modelling.
As the pulley radius R is small compared to the distance H, it is certainly advantageous to keep the finite element discretization finer near the contact zones; larger elements can be used in the free spans to reduce the number of unknowns. This is feasible with the use of the mixed Eulerian-Lagrangian kinematics, which is demonstrated below on the simple example of seeking just a single equilibrium configuration of the hanging belt.

Problem-Specific Coordinate System
Splitting the spatial coordinates into a Eulerian and a Lagrangian one for a looped belt requires a looped coordinate system. The classical polar or elliptical coordinates would have too high distortions for longer belt drives with large H and small R. Therefore, we decided to develop a problem-specific coordinate system, which consists of 6 separately defined zones, see Fig. 12.17 for the geometrical parameters H = 2, R = 0.15. (12.79) The position vector of a point in space r is determined by two curvilinear coordinates σ and ν: r = g(σ) + νn(σ). (12.80)

Here
• g(σ) is a generatrix of the coordinate system; it is a C 2 continuous curve in the plane in the plane xy, which is defined piecewise analytically over 6 segments, see details below; n is C 1 continuous, which makes r also C 1 continuous everywhere and the metric of the resulting coordinate system is then C 0 continuous; • σ is a circumferential coordinate in the direction of gross axial motion of the belt, ν is a coordinate in the normal direction.
We choose polar coordinates in the segments 1 and 4 with x > H/2 and x < −H/2 coinciding with the outer half-circles of the pulleys. The analytical expressions for other four segments involve exponents, which provide the necessary continuity and which turn into ordinary Cartesian coordinates at a significant distance from the contact zones. The coordinate lines of the coordinate system are orthogonal, and in each segment and in each point we compute scalar products ∂ σ r·∂ ν r = 0, ∂ ν r·∂ ν r = 1, ∂ σ r·∂ σ r = Λ(σ, ν). (12.82) The metric coefficients of the coordinate system form a diagonal matrix, and the component in the circumferential direction reads Λ = Λ 0 (σ) + 2Λ 1 (σ)ν + Λ 2 (σ)ν 2 , Λ 0 = g · g , Λ 1 = g · n , Λ 2 = n · n . (12.83) The metric is continuous owing to the achieved C 2 smoothness of the generatrix, and the single coefficients Λ 0,1,2 may be pre-computed in the integration points of the finite elements, in which σ = const. In segments 1, 4 with polar coordinates x ≤ −H/2 and x ≥ H/2 we have Λ 0 = 1, Λ 1 = 0, Λ 2 = 1. More complicated expressions need to be used in −H/2 < x < H/2, which, however, result in Λ ≈ 1 far away from the pulleys, where the coordinates are almost Cartesian. It remains to notice, that the generatrix g(σ) is constructed such, that in each point σ is the arc coordinate of the projection of the point onto the undistorted contour of the belt, and at the free spans we have dσ = ±dx. Finally, the function g(σ) is periodic, and we have r(σ, ν) = r(σ + L, ν). (12.84)

Mixed Lagrangian-Eulerian Kinematic Description
Particles of the belt are identified by the material coordinate s. The belt is looped, which means that coordinates s and s + L mat correspond to the same material point. At material (Lagrangian) description, we would consider the curvilinear coordinates as functions of the material one: r = r(σ(s), ν(s)). (12.85) The proposed mixed Eulerian-Lagrangian description rests upon the change of variable: s = s(σ). (12.86) We observe the material particles, which cross spatial coordinate lines σ = const. The position vector reads r = r(σ, ν(σ)).

Finite Element Approximation and Energy
As discussed in Sect. 12.3.3, the cubic finite element approximation in the domain 0 ≤ σ ≤ L (12.89) provides us with the C 1 continuous approximation of s(σ), ν(σ). (12.90) As each node has fixed coordinate σ k , we mesh all six segments of the coordinate system individually. The periodicity conditions are ensured by looping the finite element mesh, i.e., by uniting the degrees of freedom at nodes at σ = 0 and σ = L. Instead of directly using the condition (12.88) on such a mesh, we account for the jump in s on the very last finite element when computing the derivatives. In contrast to (12.9), here we use the Green-Lagragian strain measure ε = 1 2 (∂ s r·∂ s r − 1) ; (12.91) see Eliseev and Vetyukov (2012) for a discussion of the difference between the possible strain measures in mechanics of strings. Transforming (12.91) to the new kinematic description (12.90), we use the rule of differentiating inverse and implicit functions along with the metric (12.82) and compute ε = Λ + ν 2 − s 2 2s 2 ; (12.92) the prime means a derivative with respect to σ. Using again a quadratic approximation for the strain energy of the belt per its material length, we write its total mechanical energy in the field of gravity as an integral over the spatial coordinate σ: (12.93) with j being the unit vector of the vertical cartesian axial y, ρ the density per unit material length and g the free fall acceleration. Using the simple penalty formulation for the normal contact of the belt and the pulleys, we find the equilibrium by minimizing U Σ + P → min, P = Here, K is the high penalty stiffness and γ is the penetration depth, which vanishes as long as the point of the belt does not touch a pulley.

Simulation Results
In the numerical example we consider a soft rubber belt with the square cross-section of the size 5 · 10 −3 , which is by 1% longer than the contour between the pulleys, and use the following parameters of the problem: b = 12.5, ρ = 0.03, g = 9.8, K = 2 · 10 5 , λ = 1.01. (12.95) We used a non-homogeneous mesh with smaller elements near the points, in which the belt touches the pulleys. The minimization (12.94) started with the initial approximation s = σ, ν = 0 and resulted in the configuration, shown in Fig. 12.16. In the equilibrium state we have the total strain energy of the belt approximately 0.0534, the potential of the field of gravity −0.252 and the penalty term resulted into a very small value P ≈ 3.89 · 10 −5 , which means that the penalty stiffness factor is sufficiently high. In Fig. 12.18 we demonstrate the computed distribution of strain along the belt. Locking (inability of the finite element approximation to represent the exact solution) results in irregularities, which appear mostly near the upper points of the pulleys. Dividing each element into two, we arrived at a much more regular solution for a refined mesh (red line in the figure).
We additionally validated the results by integrating the contact force, which acts on the pulleys from the side of the belt. The vertical component of the integrated force, which acts on the left pulley results approximately to −0.734268, while the half of the weight of the belt is ρgL mat /2 = 0.73381. Good correspondence was also achieved for the horizontal component of the contact force at the left pulley, which results into 0.9968 after integration. Computing the strains in the middle points of the upper and lower free spans of the belt, finding the respective tension forces according to Q = bε √ 1 + 2ε (for the particularity of computing the force with the Green-Lagrangian strain measure we refer to Eliseev and Vetyukov (2012)) and adding both we arrive at 0.9941.

Work in Progress and Outlook
While it is too early to present the preliminary results of modelling the motion of the belt with rotating pulleys and dry friction in detail, it can already be stated that the resulting time histories of the strains in different points of the belt are remarkably smooth, in particular in comparison with the attempts to achieve the same solutions using Lagrangian kinematic description.
Step-like discontinuities appear when the Fig. 12.18: Computed distribution of strain along the hanging belt for the original mesh and for the refined one contact state in a neighboring integration point switches from stick to slip as the zone of creep of the belt is growing (Rubin, 2000;Reynolds, 1874). The mesh convergence is very satisfactory.
The discussed above kinematics and simulation strategy may be applied to significantly more challenging problems, like e.g. modelling the motion of an endless steel belt between two rotating drums, Fig. 12.19. Being efficient in such technological processes as food industry, production of laminates, casting of optical films, etc., such belts often suffer from lateral run off during the motion because of intrinsic unsymmetries in the geometry, orientation of drums, temperature distribution, etc. Reliable mathematical modelling, which is necessary for the model-based design of a controller (Fig. 12.20) is particularly challenging for weakly tensed belts because of the high membrane and low bending stiffness along with the very hard contact conditions. Preliminary results of a non-material finite element simulation of a steel belt using a shell model and a three-dimensional curvilinear coordinate system, which is the extension of the one in Fig. 12.17, are shown in Fig. 12.21. Imperfect geometry of the belt is accounted for using the multiplicative decomposition from the undeformed state to the reference one (Vetyukov et al, 2017a), and then we apply the mixed Eulerian-Lagrangian kinematic description to prevent the mesh from motion in the circumferential direction. We observe the effects, known from the practical experience: skew hanging of the belt, partial loss of contact with the drums (indicated by missing integration points in the figure) and lateral run-off during the motion.  Non-material finite element simulation of a moving steel, which is running off in the lateral direction because of the imperfect geometry Currently, the research work is focused on the consistency of the analysis of the effects of the dry friction with stick and slip phenomena.