Geometric Numerical Integrators Based on the Magnus Expansion in Bifurcation Problems for Non-linear Elastic Solids

We illustrate a procedure based on the Magnus expansion for studying mechanical problems which lead to non-autonomous systems of linear ODE's. The effectiveness of the Magnus method is enlighten by the analysis of a bifurcation problem in the framework of three-dimensional non-linear elasticity. In particular, for an isotropic compressible elastic tube subject to an azimuthal shear primary deformation we study the possibility of axially periodic twist-like bifurcations. The approximate matricant of the resulting differential problem and the first singular value of the bifurcating load corresponding to a non-trivial bifurcation are determined by employing a simplified version of the Magnus method, characterized by a truncation of the Magnus series after the second term.


INTRODUCTION
n the framework of three-dimensional non-linear elasticity, a number of "small on large" bifurcation problems lead to the analysis of a system of linear second order ODE's with varying coefficients.It is a common practice to reduce the system of linear second order ODE's to a simpler non-autonomous first order linear ODE system.Nevertheless, the resulting differential system may be somewhat complex and only numerically tractable; thus, it is crucial to adopt an adequate strategy for obtaining an accurate numerical solution for determining the critical load and the corresponding bifurcation deformation field.For example, approximate expressions of the matricant of non-autonomous linear ODE systems, based either on the multiplicative integral of Volterra or on the truncated Peano expansion, have been recently proposed for bifurcation problems analyzed by means of the Stroh approach within the so-called sextic formalism [1][2][3]).In Section 2, we discuss an alternative numerical method -based on the Magnus expansion [4] -which furnishes an approximate exponential representations of the matricant of first order linear ODE systems.The Magnus expansion belongs to the so-called geometric numerical integrators, based on Lie-Group methods [5 -7] for appropriate references).To our knowledge, the applications of this method are not widespread in continuum mechanics, but it has been shown in [8] that for problems involving singularities, bifurcations and wave propagation the geometric numerical integrators may be useful and accurate.In particular, the Magnus expansion may be very efficient for a number of bifurcation analyses, since it leads to the determination of approximate solutions that preserve at any order of approximation the same qualitative properties of the exact (but unknown) solution; moreover, this method exhibits an improved accuracy with respect to other frequently used numerical schemes.As an application of the Magnus method, in An application of the Magnus method in solid mechanics: periodic twist-like deformations bifurcating from the azimuthal shear of a circular tube we investigate the possibility for a compressible hollow cylinder subject to I azimuthal shear to support axially periodic twist-like bifurcation modes similar to the Taylor-Couette axially periodic cellular patterns observed when the flow of a viscous fluid confined in the gap between two rotating concentric cylinders becomes unstable.Here, the hollow cylinder is made of a Levinson-Burgess isotropic compressible elastic material.By restricting our attention to a class of incremental displacements characterized by three unknown scalar functions of the radial coordinate and having an axial periodic structure, we then study the related incremental boundary-value problem.This leads to a somewhat complex non-autonomous homogeneous system of six first order linear ODE's with homogeneous boundary conditions, whose analysis is performed by a numerical approach.We show that a procedure based on the Magnus method outlined in An overview of geometric numerical integrators with special reference to the Magnus expansion reveals very effective either for approximating the matricant of the resulting differential problem or for determining the first singular value of the bifurcating load corresponding to a non-trivial twist-like solution.

EXPANSION
number of problems in scientific and engineering areas, ranging from Quantum Mechanics, Optics, Electromagnetism and Molecular Physics to Control Theory and Bifurcation, require the resolution of systems of differential equations with varying coefficients.Besides some very special cases, for which a fundamental matrix solution is explicitly available, the analysis of the problem usually requires the construction of an approximate representation of the solution of the non-autonomous differential system.In the last decade, this issue has attracted many scientists and has also represented a field of intense research activity; in particular, recent studies have shown that in many cases the so-called geometric numerical integration methods may offer better approximation schemes compared to perturbative procedures or to standard numerical integrators.The Magnus method belongs to the class of geometric numerical integrators, based on Lie-Group methods (see [6,7] for appropriate references), which, although not yet diffused in continuum mechanics, have recently found many applications in areas of molecular physics [9], electromagnetism [10] and celestial mechanics [11].In this Section, we first give a brief overview of geometric numerical integrators and then we deal with the Magnus method as a particular case.Geometric numerical integration concerns the development of numerical approximations to the solution of an ODE system with the main goal of preserving at any order of approximation the qualitative properties of the exact (but unknown) solution and thereby of exhibiting an improved accuracy with respect to other numerical methods.Although such methods apply also for nonlinear matrix differential equations, for the purposes of bifurcation problems here we consider geometric numerical integrators based on the Magnus expansion for the special case of homogeneous nonautonomous first order linear ODE systems of the form where the unknown 6x6 matrix   R Y , usually referred as the matricant of the Cauchy initial-value problem (1), is the 6dimensional identity matrix when the real independent variable R is equal to zero.Generally speaking, these new numerical approaches involve Lie-Group methods, which are based on Lie groups and their associated Lie algebras.We recall that a Lie group corresponds to a differential manifold which is endowed with a group structure; the corresponding Lie algebra, defined as the tangent space to the Lie group at the identity, is a linear space endowed with the Lie bracket commutation operation and a natural mapping.Since the unknown matricant of a first order ODE system like (1) must be an invertible 6x6 matrix with variable entries, the setting for finding solutions of (1) is a 6x6 matrix Lie group.In this case, the above general definition of a Lie algebra implies that the 6x6 matrix   R A governing (1) is an element of a matrix Lie algebra.Moreover, in this setting the associated matrix Lie algebra is endowed with the matrix commutation law   where   , A B is commonly called the matrix commutator.Finally, for these special matrix Lie groups the matrix exponential function plays the role of the natural mapping.We emphasize that in the theory of Lie groups the role of the exponential mapping is crucial.Indeed, within the context of differential geometry, the matrix exponential maps an element belonging to a matrix Lie algebra into an invertible matrix belonging to its corresponding Lie group; in other words, the exponential mapping allows to recapture the local A group structure from the corresponding Lie algebra.Thus, for a first order ODE system like (1), it is natural to attempt to construct the solution as a matrix exponential with exponent given by a matrix which belongs to the associated Lie algebra.This peculiarity of the exponential mapping probably inspired the idea developed by Magnus [4], who proposed to express the exponent by a matrix series expansion made of infinite terms, each of them belonging to the Lie algebra corresponding to the Lie group of the unknown matricant.In detail, with reference to the differential problem defined in (1) the Magnus method assumes the matricant to be of exponential form with the exponent given by the following series expansion, called the Magnus expansion: Here we do not report the analytical procedures for determining the terms of the series expansion in (3).We address the reader to the seminal paper [5] where, in view of ( 1) and ( 3), and according to the notions of the first derivative of a matrix exponential map and of its inverse (see formulae (33) and (38) in [5]), it is clearly showed how it is possible to construct all the terms of the Magnus series.For our purposes, we report here only the first two terms of the series: where [. , .] is the matrix commutator (2) and    A is the governing matrix of (1).Notice that each of the remaining terms (here not reported) contains nested matrix commutators of the form (2) involving the matrix operator   holds, then, in view of ( 3), ( 1) and ( 4), the matricant of (1) takes the exponential form which is analogous to the classical exponential solution for a scalar linear ODE.Clearly, (6) yields the solution of (1) if the entries of the matrix A governing (1) are independent of R. Obviously, (5) holds only in special cases because the matrices in general do not commute (this is closely related to the structure of the commutation law (2) which characterizes matrix Lie algebras).
A crucial aspect of the Magnus method related to its pertinence when seeking approximate solutions of differential problems like (1) emerges from the analysis of ( 3) and (4).By the above discussion, we recall that the matrix   R A in (1) belongs to a matrix Lie algebra for all R; thus, in view of (4) we see that any term of the Magnus expansion belongs to the same Lie algebra.This implies that any truncation of the Magnus series will also belong to the same Lie algebra and therefore the exponential map of any truncation will necessarily stay in the corresponding Lie group.This is basically the main reason why an approximated solution obtained by truncating the Magnus expansion at any order preserves the same qualitative features of the exact solution.
Another major issue concerns the conditions on the matrix   R A which guarantee the convergence of the Magnus series.This problem has been studied in [12], where it is shown that a sufficient condition for the local convergence of the Magnus series in a certain interval   2 0, R is given by : where the integrand function in (7) is the spectral norm of   R A , i.e., the square root of the maximum eigenvalue of the . Moreover, the fulfilment of the inequality (7) determines the range of the independent variable R for which it is possible to write the solution in the form (3).However, in many applicative cases it happens that the convergence condition (7) does not hold for the whole integration interval   2 0, R .In these cases, as usual for numerical integration methods, it is common dividing the interval   is expressed by the infinite Magnus series.Consequently, the matricant of (1) evaluated on the whole As it usually occurs when constructing approximate solutions, a crucial issue is that of choosing an appropriate order of which guarantees a suitable level of accuracy.To this aim, we recall the efficient procedure suggested in [7], which may be summarized in the following three steps: first, the series   i+1 R  is truncated at a suitable order in each subinterval wherein the convergence is guaranteed.Then, the integrals appearing in each truncated series are approximated by appropriate quadrature rules.Finally, the exponential of the matrix   i+1 R  is evaluated.The qualified literature on this issue contains a number of proposals which correlate the level of accuracy of the method to the order of truncation of the Magnus series.We recall the paper [6] for an overview on this problem.Finally, the accurate calculation of the exponential of a matrix is still an open problem, as one may infer from the paper [13].In the spirit of geometric numerical integration methods, the consequence of an inaccurate computation of the exponential matrix due, for example, by the replacement of the exponential map with one of its analytic (for example rational) approximants may lead to a solution not belonging to the desired Lie group, thus rendering the main advantages of geometric numerical integration ineffective.Here, we report two references on this issue, namely the splitting method proposed in [14] and the scaling and squaring with a Padè approximation method summarized in [13].Notice that the command MatrixExp of the widely used software Mathematica employs the Padè approximation method for the calculation of the exponential matrix.We conclude by observing that the non-trivial evaluation of the exponential matrix may be circumvented when the unknown matrix of a differential problem belongs to special Lie groups.One of these cases involves a group commonly occurring in mechanical problems, namely the orthogonal group, that is the group of 3x3 matrices Q such that T 3x3 . QQ I Unlike many other Lie-group solvers, they do not require the evaluation of matrix exponentials because, as outlined in [15], the back translation from the Lie algebra to the corresponding Lie group may be obtained by using the Cayley transform.

AN APPLICATION OF THE MAGNUS METHOD IN SOLID MECHANICS: PERIODIC TWIST-LIKE DEFORMATIONS BIFURCATING FROM THE AZIMUTHAL SHEAR OF A CIRCULAR TUBE
n the present section, we introduce a bifurcation problem set within the context of the nonlinear theory of elastic solids.In particular, the last part of the following analysis contains an explicit application of the Magnus method based on the procedure outlined in the previous section.

The fundamental azimuthal shear equilibrium deformation
We consider a homogeneous, isotropic, compressible, elastic tube which occupies the region in its natural reference configuration.  R, , Z  denote the cylindrical coordinate of a point X in a cylindrical coordinate system with orthonormal basis   R Z , ,  e e e .The boundary of  is divided into two disjoint complementary parts:  of  is assumed to be a smooth function with gradient     :  F X f X .Since the tube is elastic and isotropic, the general form of the strain energy function is given by B , where are the orthogonal principal invariants of T :  B F F ; thus, the Piola stress takes the form with , W : , W : We assume that the inner cylinder is kept fixed, whereas the outer is subject to a uniform angular displacement 0   around its axis.On the bases of  only tangential displacements are admitted.This leads to the following mixed boundary-value problem: For the equilibrium problem ( 15)-( 17), we consider the possibility of an azimuthal shear deformation f  defined by where  , the angular displacement field, is assumed to be a smooth function satisfying Consequently, in view of (18) we have where the prime denotes differentiation with respect to R and    .Notice that the azimuthal shear is an isochoric deformation whose principal invariants (12) are given by Moreover, (18) trivially satisfies the displacement boundary conditions ( 16)-( 17) 1 and, in view of (13) and (20), it is easily seen that also the traction boundary condition (17)2 holds.It remains to study the equilibrium field Eq. ( 15).The explicitly evaluation of ( 15) by means of (13) and ( 20) yields an over-determined system of two ODE for the single unknown function   R  . In particular, here we consider the class C 2 Levinson-Burgess strain energy function where and are the referential shear modulus and referential Poisson's ratio, respectively, and is a dimensionless material parameter.For the constitutive class (22), the Piola stress takes the form By following a result contained in [16], it is possible to generate classes of strain energy functions such that the azimuthal shear satisfies both the two equilibrium ODE; in particular, for the present case of a Levinson-Burgess material, deformations of the form ( 18) are allowed only if the following condition on the material parameters holds: Finally, in view of ( 23), ( 19)-( 21) and ( 24), after a non-trivial calculation we obtain from (15) a single ODE, whose solution is

Periodic twist-like bifurcations
We now consider the possibility of deformations which bifurcate from the primary equilibrium pure circular shear (18) as the loading parameter (angle)  increases from 0 (natural state).A condition (only necessary) for the occurrence of a local branch of bifurcating solutions at a fixed  is the existence of a non-trivial solution of the homogeneous linearized equilibrium problem which corresponds to an adjacent equilibrium state.Such solutions superposed to the primary azimuthal shear must satisfy the adjacent equilibrium Eq. ( 26) and the incremental boundary conditions ( 27)-( 28), obtained by linearizing ( 15)-( 17) around the fundamental deformation    x f X  , which is conveniently assumed as the independent variable: In particular,  

:  u x 
 represents an incremental displacement field, "div" and "grad" are the divergence and gradient operators with respect to x, respectively, and  is the fourth-order instantaneous elasticity tensor, given by: For the constitutive class ( 22), it follows from ( 29), ( 20)-( 21) and ( 24) that where B  is given by (20)2.We now restrict our analysis to the class of periodic bifurcating displacements which are reminiscent of the instability patterns observed in the classical Taylor-Couette shear flow of viscous fluids.A more detailed discussions on such a choice, together with an exhaustive analysis of problems closely related to the one here presented, are developed in [17], [18]; in particular, in [17] one may also find stability issues based on [19] and partly on [20].This allows us to skip here the description of certain analytical developments, for which we refer to the works above mentioned.Notice that (31) models the occurrence of an axially periodic cellular pattern in the gap between the inner and outer cylinders (n represents the number of possibly forming cells in the axial direction); in particular, inside each of the n possible forming cells, (31) describes a twist-like displacement perpendicular to the  e -direction of primary annular shear, so that the vector lines of the incremental displacement u are similar to the streamlines of the twisting Taylor-like effects for fluids.We easily check from (31) that the displacement incremental boundary condition (28) 1 trivially holds, whereas (27) holds provided the smooth scalar functions   For what concerns the traction boundary condition (28)2 and the field Eq. ( 26), we first observe by (31) that ; thus, in view of (30)-(31), condition (28)2 immediately holds.Finally, by evaluating the divergence of (30) with the aid of (20) 2 and (33), after some non-trivial calculation and also by employing separation of variables allowed by the choice (31) of the displacement field, we reduce the partial differential problem ( 26)-(28) to the following system of three homogeneous second-order ODE (see also (32)): where   are non-constant matrices determined by the components  of the symmetric fourth- order tensor  in the coordinate system (r, , z).Here, we do not report the explicit expression of such components, which may be obtained by (30), (20) 2 and (26).We only observe that   The latter property of   r,  P is crucial.Indeed, a common practice is that of transforming the set of three linear second order ordinary differential equations like (34) into a system of six linear first order ordinary differential equations, for which well-known existence theorems and procedures for constructing the solutions are available in the literature.With this aim, once introduced the 3x3 matrices we rewrite (34) in the form , r r r r r Then, by setting 6x1 6x6 6x6 6x6 : , : : , : we easily see that (37) is equivalent to the linear homogeneous first order ODE boundary-value problem Observe that in the notation above we have dropped the dependence of the basic parameters (the inner and outer radii, the height of the tube, the material moduli and the angular displacement  ), usually assumed as prescribed.We only maintain explicit the dependence on r, in order to emphasize that the linear ODE system is non-autonomous.Now, by following elementary results on systems of linear ODE's, we may write a solution of (39) in the form (see (38)1 and (37)2) where   1 r y is an unknown constant vector and , with det r 0, for r < r < r and r r r is a particular fundamental matrix solution for (39)1, usually called matricant.Then, substitution of (40) into (39)2 yields which, in view of (38) 3-4 , (41) and (40), is equivalent to the homogeneous system of three algebraic equations for which non-trivial solution   Since (46) has the same structure of the Cauchy initial-value problem (1), in the following we illustrate a strategy based on the Magnus method for determining the solution of (46).

A strategy for determining the bifurcating load through the Magnus method
The theoretical preliminaries in An overview of geometric numerical integrators with special reference to the Magnus expansion combined with the bifurcation analysis of Subsections The fundamental azimuthal shear equilibrium deformation -Periodic twist-like bifurcations allow us to describe a strategy for determining the critical value cr  of the primary angular displacement, that is the lowest value of the angle  for which the tube may support (during a loading process) a non-trivial twist-like bifurcation of the form (31) from the fundamental azimuthal shear deformation.The Magnus method is now the key ingredient for approximately determining the matricant of (46), whereas (44) plays the role of bifurcation condition.Following the steps outlined in the Section 2, we first fix the radii R1 and R2, the height H, the referential shear modulus  and the referential Poisson's ratio  (recall that 3 / 4   by ( 24)).Then, we: -fix the number n of twisting cells in the axial direction (by (31) this is equivalent to fixing  ) starting from n =1; -fix the value of the loading parameter  starting from zero; -compute the matrices (35), which now depend only on r, and then compute   after the second term in each subinterval (see ( 4)) and by approximating the integrals using a Gauss quadrature rule.This allow us to write where are two Gauss points in correspondence of which   r A and its employed linear approximation agree.
compute the matrix exponential of   i+1 r  using the function MatrixExp in the software Mathematica; -evaluate the matricant of (46) at the external radius r2 = R2 by ( 8) and (9); -compute the "north-east" block of the matricant at r2 = R2 , that is the 3x3 matrix   Following this procedure, we define n cr  as the smallest value of  which satisfies the bifurcation condition (44) corresponding to the fixed value of n.Of course, when one repeats the above procedure by varying the number of possible axial cells, the corresponding critical load n cr  varies.We have performed a large number of computations which show that as n increases, the corresponding n cr  determined by (44) decreases until a critical value of n above which (44) definitively does not hold for any  .In correspondence, a load cr  corresponding to cr n is determined, and this is the critical value of the circular angle of shear which allows for the occurrence of a non-trivial axial periodic deformation of the form (31). R /R between the outer and the inner radius of the tube.We may easily infer, as expected, that twist-like bifurcation may be more easily supported in thin tubes.

CONCLUSIONS
e have studied a bifurcation problem concerning the possibility of periodic twits-like deformations superposed to a primary azimuthal shear state for a Levinson-Burgess isotropic elastic tube.The procedure proposed for determining approximate solutions of the resulting non-autonomous systems of linear ODE's is based on the Magnus method, which represents an innovative numerical scheme belonging to the more general class of geometric numerical integrators.We have the feeling that these new methods, already employed in many scientific areas, may find interesting applications also in the field of continuum mechanics.In the future we are going to further explore the applicability and the advantages of the Magnus method by studying other bifurcation problems yielding non-autonomous ODE systems like, for example, a Taylor-Couette type bifurcation for a sheared incompressible annular tube modeled by an extended version of the Gent constitutive equation.

2 0
, R into N subintervals   i+1 iR , R such that the Magnus series converges in each of them.Given such a subinterval   i+1 i R , R , in view of (3) we may write: since  is symmetric; moreover, for the constitutive class (22)   r,  P is always invertible.
solution of the primary ODE problem (34) can be written in the following form the boundary conditions (34)2 (see (41) and (43)).The above analysis shows that the possibility of determining non-trivial solutions for the bifurcation differential problem (34) is strictly related to the knowledge of the 3x3 matrix   2 r U , which represents the "north-east" block of the 6x6 matricant (41) of the Cauchy initial-value problem such that the Magnus series converges in each of them (see the convergence condition (7)); -employ for simplicity a fourth-order method by truncating the Magnus series   i+1 r  the bifurcation condition (44) holds.

Fig. 1
reports the bifurcation curve for a referential shear modulus 1   MPa, referential Poisson's ratio 0H/R 1 =1.It shows the critical load cr  as a function of the ratio 2 1