Bifurcations of an elastic disc coated with an elastic inextensible rod

An analytical solution is derived for the bifurcations of an elastic disc that is constrained on the boundary with an isoperimetric Cosserat coating. The latter is treated as an elastic circular rod, either perfectly or partially bonded (with a slip interface in the latter case) and is subjected to three different types of uniformly distributed radial loads (including hydrostatic pressure). The proposed solution technique employs complex potentials to treat the disc’s interior and incremental Lagrangian equations to describe the prestressed elastic rod modelling the coating. The bifurcations of the disc occur with modes characterized by different circumferential wavenumbers, ranging between ovalization and high-order waviness, as a function of the ratio between the elastic stiffness of the disc and the bending stiffness of its coating. The presented results find applications in various fields, such as coated fibres, mechanical rollers, and the growth and morphogenesis of plants and fruits.


Introduction
When an elastic cylindrical shell or a circular rod is subject to a radial external pressure of sufficient intensity, buckling occurs via ovalization, as shown with a teaching model in figure 1 on the left (where a segment of an acrylic polymer tube with a diameter of 42 mm and a wall thickness of 0.5 mm is used, maintaining a thickness-to-diameter ratio equivalent to that of a chicken eggshell).This buckling phenomenon is well known and has been analysed by various researchers, among others [1][2][3][4].In particular, a distinction has been introduced between three different mechanical models for the external uniform radial load [5,6]: (i.) 'Hydrostatic' or 'pressure' load, which always remains orthogonal to the structural element to which it is applied in any configuration (undeformed or deformed); (ii.) 'Centrally directed' load, which acts on the structural element remaining always directed towards the initial centre of the ring; (iii.) 'Constant directional' or, better, 'dead' load which remains aligned parallel to the unit normal to the structural element to which it is applied in its undeformed configuration.
All three above loads are conservative [5] and the difference between them emerges in the incremental equations, holding for departures from the trivial configuration, so that they lead to remarkably different bifurcation loads.Bifurcation also occurs in the case when an elastically deformable core is present inside of the shell (or the circular rod), as shown with a teaching model in figure 1 on the right (where Figure 1: A teaching model illustrates buckling in cylindrical shells.Left: a thin shell (undeformed in the upper part, 42 mm in diameter, 0.5 mm wall thickness made up of plastic) buckles under external pressure when air is extracted from it (lower part), resulting in ovalization.Right: buckling is also induced in a cylindrical shell with an inner core (a tubular rubber foam); in this experiment, the core's stiffness is insufficient to induce bifurcation in a small-wavelength mode, so that in both cases ovalization occurs.Note the position of the plunger inside the syringe showing that the pressure load is higher in the experiment with the internal core.The simple experimental set-up shown on the right highlights the idea beyond the mechanical modelling, in which a shell becomes a coating enclosing an elastic core.
the core is a rubber foam used for pipe insulation).The simple experimental set-up inspires the mechanical modelling that will be adopted in the present article, where the coating of the elastic core is provided by an axially inextensible elastic rod.This model is already well developed in mechanics [7][8][9], including cases involving circular and elliptic geometries [10][11][12][13][14][15][16][17][18][19], and has also been used to analyse bifurcation, but only with respect to a half space [20,21], while bifurcation for circular configurations has never been analysed.
The inner core inside the circular rod not only increases the buckling load (note the position of the plunger inside the syringe, which shows that the pressure for buckling is higher when the inner core is present), but also results in a complex bifurcation behaviour that allows for the emergence of short-wavelength wave modes, although this does not take place in the case shown in figure 1, due to insufficient stiffness of the core.The bifurcation problem of a coated elastic disc has been scarcely analysed and always only for hydrostatic pressure loading (i): in [22,23] (where a numerical, not analytical, solution is only found), under the assumption that the coating, modelled as an elastic rod, cannot transmit shear stress to the inner disc (the imperfect bonding condition also considered by us); in [24], where the coating is modelled as an elastic shell, either fully or partially bonded to the disc, and the obtained solution is analytical, though based on a number of mathematical simplifications.Therefore, the bifurcation of coated discs under radial forces is still an open and almost unexplored problem.
The present article employs the coating model for an elastic disc introduced in [25], which can be considered a specific case of the shell-coating model formulated in [26].This model assumes that the elastic disc is coated with an elastic rod that is axially inextensible and unshearable, thereby introducing a Cosserat and isoperimetric constraint for the disc.The analysis of buckling is carried out with all the three load variants (i)-(iii).Prior to bifurcation, the rod exhibits trivial equilibrium, characterized by a pure axial internal force.As a result, the elastic core, assumed to be isotropic, remains unloaded before bifurcation and follows incremental equations obeying linear isotropic elasticity, characterized by Lamé constants λ d and µ d .Two transmission conditions are analysed for the bonding between the elastic rod and the inner core, namely, perfect bonding, where full continuity of radial and tangential displacements is enforced, and tangential slip contact, in which the radial displacement is transmitted but the tangential is not, so that the shear stress at contact remains null.The latter condition can capture the behaviour of a partially detached coating.It can also model a coating attached at discrete points to the disc, as is the case of the piece (manufactured with a Stratasys J750 3D printer, following the multi-material Polyjet technique process, with a layer's printing resolution of 27 µm) documented in the photo reported in figure 2.  , 44 mm in diameter) is connected to a stiffer photopolymer (VeroYellow ™ , 1 mm thick) coating through radial ligaments.The latter can transmit axial but only negligible transverse forces.Therefore, a model of slip interface can be adopted to model the disc (manufactured with a Stratasys J750 3D printer) when radially loaded on its external surface.The piece is bioinspired by the joining through ligaments of the brain to the cranial vault.
The mechanical model adopted in the present article allows for an analytical solution to the bifurcation problem in a simple closed form, using Kolosov-Muskhelishvili complex potentials for the core and incremental Lagrangian equations for the coating.The superiority of the complex formalism becomes evident by comparing our closed-form and exact solution with the approximate solution provided in [24].The analysis demonstrates that imperfect bonding decreases the bifurcation threshold and that the hydrostatic-pressure model (i) results in the highest critical loads, while the centrally directed load model (ii) leads to the smallest.Significant differences in the bifurcation conditions discovered here highlight the importance of the role played by mechanical modelling of applied loads for accurate representations of load transfer mechanisms under structural deformation, a topic often underestimated and given here a new evidence.
Results of the bifurcation analysis show that, in cases where the inner core is sufficiently compliant, the critical bifurcation mode corresponds to ovalization.However, as the ratio between the stiffness of the elastic core and that of the coating rod increases, the modes display all possible waviness until they approach the vanishing-wavelength condition in the limit, where the stiffness ratio approaches infinity.
Coatings are widely used in various technologies, making the findings of this article applicable to several areas.For instance, the results may be useful in the design of mechanical rollers, as well as of coated fibres, at both the micro and nanoscales.Another exciting application is to the morphogenesis and growth of plants and fruits.In such cases, turgor pressure can reach up to 10 atmospheres, which is sufficient to trigger various types of bifurcation.For example, it can produce ruffle-like or dome-like shapes in leaves [27] and undulations in flowers, characterized by annular geometry [28].These findings suggest that coatings may play a significant role in shaping and enhancing the growth of plants and fruits, which can have implications for agriculture.In fact, the coated disc analysed here exhibits bifurcation modes that result in elegant shapes, resembling those observed during the maturation of some fruits or vegetables.In these cases, a soft pulp is enclosed in a stiff husk, so that pressure generated during drying may lead to the formation of gracious undulations.A prime example is shown in figure 3a, which illustrates a pumpkin (white swan acorn variety) and its cross section.At the initial growth stage, the fruit appears smooth on the outside.However, during maturation, the inner pulp dries up and generates a state of compression in the stiff husk, causing it to buckle and resulting in a wavy surface.Drying phenomena are highly complex as the case of colloidal drops shows [29,30], nevertheless our results can provide some insight in the maturation problem of the pumpkin.In particular, the bifurcation of a coated disc reveals that the presence of an inner core leads to complex undulated bifurcation modes, rather than the simple ovalization that occurs in the absence of any inner reinforcement.Figure 3b showcases the outcome of our bifurcation analysis for a coated elastic disc, illustrating the incremental displacement field (colours evidence growing values from the centre to the boundary) for a mode characterized by 10 wavelengths, generated by a coated elastic disc (with a ratio R/h = 7.75 evaluated for the pumpkin, where R is the radius of the disc and h is the thickness of the circular rod).By selecting for the pulp of pumpkin a Poisson's ratio equal to 0.5 and imposing the observed waviness of the skin, the bifurcation analysis permits the evaluation of the ratio between Young's moduli of the skin and of the pulp, which is found to fall between 0.543 (a) (b) Figure 3: (a) A pumpkin (white swan acorn variety) with a diameter of 20 cm, featuring circumferential waviness resulting from the drying of its inner pulp during maturation.(b) A bifurcation analysis of a coated elastic disc provides some insight in the complicated problem of pumpkin ripening.In particular, the incremental displacement field (where colours denote intensity increasing fromthe centre to the boundary) of a coated circular disc at bifurcation, involving amodewith 10 wavelengths (with a measured ratio between the disc's radius and the coating's thickness of 7.75, and a Poisson's ratio of 0.5), enables the evaluation of the ratio between the Young's moduli of the pulp and that of the skin.This ratio is found to be equal to 0.6, a value that is not easily determinable otherwise.and 0.736 (the figure is generated at 0.6).Therefore, while a direct measure may not be easy, our closed-form solution allows to determine the stiffness ratio between husk and pulp by simply counting the number of external undulations without cutting the vegetable.
2 Large deformation of a planar elastic rod

Kinematics of a curved rod
Figure 4a shows a rod that has undergone deformation in a plane defined by two unit vectors, e 1 and e 2 .The reference configuration of the rod is characterized by the arclength parameter s 0 , while its current configuration is characterized by the arclength parameter s.The points on the  rod in the reference and current configurations are parametrized, respectively, as x 0 (s 0 ) and x (s) and they are related to each other through the deformation g(x 0 ) and its inverse g −1 (x), The displacement of a point on the rod is defined as which can be expressed as a function of either s 0 or s.
The unit tangents, t 0 and t, principal normals, n 0 and n, and curvatures, κ 0 and κ, are defined in the reference and current configurations, respectively, by When axial deformation of the rod is negligible and thus axial inextensibility is enforced, the stretch λ, defined as the ratio between the strained and referential elements ds and ds 0 , becomes unity, ds = ds 0 .In this case, the geometrical elements (3) are related to each other through the equations It is instrumental to define unit normals to the rods in the two configurations, m 0 and m, which coincide with the principal normals, except possibly for the sign, where e 3 = e 1 × e 2 is the out-of-plane unit vector.The derivative of equations ( 5)with respect to the arclengths s 0 and s, respectively, leads to which can alternatively be expressed as The unit vector m in equation ( 5) and its derivative can be written in terms of referential quantities as 2.2 Statics of a curved rod

Equilibrium in the current configuration
An element of the rod in its current configuration, ideally 'excised' between the arclengths s 1 and s 2 , is subject to a distributed load q(s) and moment µ(s)e 3 .To mantain equilibrium, internal forces a(s) and bending moment M (s)e 3 must act at the ends s 1 and s 2 , as shown in figure 4b.The translational equilibrium of the rod is expressed as while the rotational equilibrium is where x − o represents the position vector of a generic point x on the rod.Equations ( 9) and ( 10) can be reduced to a unique integral, which can eventually be localized to yield The internal force a can be defined in terms of axial and shear components, N and T , both referred to the current configuration, as Substitution of equation ( 12) into equations ( 11) leads to the equilibrium equations, holding for any curved rod.

Equilibrium in the reference configuration
The equilibrium equations ( 13) are now re-derived in the referential description, by introducing the nominal, or 'Piola', internal force a 0 , defined in a way that and with referential, axial and shear force components Note that the bending moment M remains unchanged in the reference configuration, say, M 0 = M .The same procedure used in the spatial treatment of the equilibrium leads now to where the external loads in the reference configuration remain unchanged with respect to the deformed configuration because of the validity of the inextensibility constraint.Equation ( 5) yields so that substitution of the expression (15) into equations ( 16) leads to the equilibrium equations for a curved rod in the referential description It should be noted that equation (19) 3 can be rewritten as

Constitutive equations
Constitutive equations cannot determine the rod's axial and shear forces, which are to be understood as reactions to the inextensibility and unshearability constraints.However, the bending moment, which is independent of the referential and spatial distinctions, M = M 0 , is determined by the difference in the curvature where B = EJ is the bending stiffness, equal to the product between Young's modulus, E, of the rod and the second moment of inertia of its cross section, J, and ω and ω 0 are the angles between the tangent vectors in the current and reference configurations and the axis e 1 , so that t•e 1 = cos ω and t 0 • e 1 = cos ω 0 .

Incremental equations for a curved rod
The incremental form of the equilibrium equations for the elastic rod can directly be obtained from equations (13), holding for finite deformation, as where increments are denoted with a superimposed dot.Since the increments of referential quantities are null, ṫ0 = ṁ0 = 0, the incremental form of equations ( 4) and ( 5) become Hence, equations ( 22) can be rewritten as which are equivalent to equations (22), but expressed in terms of incremental displacement u, instead than ṫ and ṁ.
The incremental versions of the equilibrium equations (19) find their counterpart in the reference configuration as The incremental version of the rotational equilibrium equation (25) 3 can be rewritten as while the incremental version of the constitutive equation ( 21) is 3 The annular rod

Governing equations for a circular rod
The theory described above is applicable to rods of any shape, assuming that they are axially inextensible and unshearable, and is particularized now to the case of a circular rod of radius R and centered at point O, assumed as origin of a reference system defined by unit vectors e 1 and e 2 .The circumferential angle θ is measured positively in counter-clockwise direction, as depicted in figure 5. Due to the polar symmetry being s = Rθ, the tangent t 0 and the normal m 0 become while the position of a generic point x 0 is singled out by x 0 = Rm 0 .Hence, for the annular rod the following relations hold true: Equations ( 13) 1−2 govern the translational equilibrium of the rod in its spatial configuration where κ represents its deformed curvature.The derivative of the tangent vector t with respect to the arclength s becomes so that the curvature is to be used in equations (13) 1−2 , which are complemented by the rotational equilibrium, equation Equations ( 29) can be used to express the partial derivative of the internal 'Piola' force (15) as The displacement vector u can be expressed in polar components as so that, according to relations (29), the derivative with respect to s of the expression (33) leads to Using equations ( 32) and (34), the equilibrium equations for the rod in the material configuration are Kinematic considerations lead to the evaluation of the axial strain ϵ (to be set equal to zero because of the axial inextensibility), the rotation of cross-section Φ, and the change of curvature χ as Imposing rod's inextensibility, ϵ = 0, and using equations (36), the equation (35) 3 can be simplified to When an external radial load with uniform distribution, Π, is applied, the annular rod is only subject to a uniform internal compressive force thus, before bifurcation, the rod remains in its circular configuration without suffering any deformation.The incremental quantities ( 23) and ( 34) assume now the form so that the spatial equilibrium is governed by equations (24), where, from equation ( 31), the increment in the curvature κ is Taking the material time derivative of the equations (35) 1−2 , the incremental translational equilibrium equation in the reference configuration is obtained, so that equation (38) 1 leads to thus, from the constitutive equations ( 27) and (36) 3 , the left-hand side of equation (41) 3 can be rewritten as The use of relation (42), when combining together equations (41), yield the differential equations describing the kinematics of the annular rod, subject to an external uniform radial load Π where

The incremental applied load
Equation ( 43) requires the evaluation of the increment q in the applied external radial load, taking into account the direction assumed by the load in the deformed configuration [4,5,[31][32][33][34][35]. The modulus Π of the radial force is assumed to remain constant, so that Π = 0.In particular, as mentioned in the introduction, the three following different cases can be distinguished: • Hydrostatic pressure (i): the applied load remains aligned parallel to the normal m to the deformed rod element (figure5a), so that the load q h and its increment qh become • Centrally directed load (ii): the applied load remains directed toward the initial centre of the circular rod (figure 5b), so that the load q l and its increment ql become • Dead load (iii): the applied load is dead and does not change its original direction m 0 (figure 5c), so that the load q k and its increment qk become The cases of hydrostatic pressure (i), centrally directed (ii) and dead (iii) load, [2,3,31,33,36] are recovered by setting S = S Π in equation ( 43), being uθ for centrally directed load (ii.), 0 for dead load (iii.). ( Considering a circular elastic rod (without any inner part), the critical loads for bifurcation (i.e. the smaller buckling load) for the three different loading cases can be derived by imposing continuity of displacement, bending moment and shear force, via equation ( 43), solved for non-trivial solutions.These critical loads are given by for hydrostatic pressure (i), 9/2 for centrally directed load (ii), 4 for dead load (iii). (49) Note that the case of hydrostatic pressure does not require any external constraint for equilibrium in the undeformed and deformed configuration.This is different for the other two loadings (ii) and (iii), where although the undeformed configuration is of equilibrium, the latter is unstable, so that a constraint imposing null translations has to be enforced for case (ii) and a constraint imposing vanishing rotation about the centre has to be enforced for (iii).

Bifurcation of the coated elastic disc
The core of the coated disc is elastic, and for an elastic material under large strain, the constitutive equations, relating the Cauchy stress σ to the left Cauchy-Green deformation tensor B = FF T (where F is the deformation gradient), can be written as, [37] where the coefficients β j (j = 0, 1, 2) are functions of the invariants of B.
The elastic material inside the disc remains unstrained and unstressed up to bifurcation so that B = I and σ = 0.At bifurcation, the incremental relation between the Piola stress S and the Cauchy stress leads to The material time derivative of equation (50) eveals that the elastic response inside the disc follows the usual linear elastic relation, where the Lamé constants can be calculated as where I 1 and I 2 are the first and second invariants of the Cauchy-Green deformation tensor, respectively.Consider an elastic, homogeneous and isotropic disc characterized by a radius R, shear modulus µ d and Poisson's ratio ν d with reference to a Cartesian coordinate system (e 1 , e 2 , e 3 ) with origin O placed at the centre of the disc.A polar reference system (e r , e θ , e 3 ) is also introduced, so that the displacement for generalized plane conditions can be written as where u d r and u d θ are the radial and tangential displacement components.The disc is coated on its boundary by the previously introduced rod with a bending stiffness B. In its reference configuration and loaded with an external radial load Π, the annular rod is subject to an axial internal force N 0 = −ΠR, while the interior disc remains unstressed.At bifurcation, a non-trivial incremental Figure 6: A circular inextensible rod coating elastic disc.The bonding between the two is modelled as either perfect or allowing tangential frictionless slip.Loaded by a uniform external radial load, the rod is subject to only an axial force, but at bifurcation incremental internal forces develop so that tractions are transmitted from the external coating to the disc.
deformation occurs, causing the disc to experience incremental stress and strain.The resulting incremental traction at the disc's boundary, multiplied by its thickness b (to be set equal to the unity for plane strain), gives rise to an incremental force acting on the rod, denoted as qσ (figure 6).Hence, the incremental load on the coating is given by where qΠ represents the incremental contribution associated with the external radial load which is determined from the incremental radial and tangential components of the first Piola-Kirchhoff stress tensor S, evaluated on the disc's boundary r = R.The term M in equation ( 55) describes the shear transmission properties at the interface, so that M = 1 for perfect bonding between disc and coating or M = 0 for slip contact, when shear force is not transmissed.As a consequence, equation (43) becomes where the superscript 'c' stands for 'coating' and

Complex potential formulation for the elastic disc
In a region enclosed by a sufficiently smooth and non-intersecting curve L, each point can be identified with a complex number z = x 1 + ix 2 , where x 1 and x 2 represent the coordinates of the point and i denotes the imaginary unit.Each point can also be represented in terms of polar coordinates (r, θ), where r denotes the distance from the origin to the point, and θ is the angle between x 1 and the radius r (measured positively in the counter-clockwise direction), such that z = re iθ .
The complex displacement u d (z) in Cartesian coordinates at the point z inside the disc and the complex combination of the normal and shear traction components σ d (z) at that point can be obtained fromthe complex variables forms of the Somigliana's identities, which are the corollaries of Betti's reciprocal theorem [38].The corresponding expressions are where a bar over a symbol denotes complex conjugation, τ = Re iθ ∈ L, and u d (τ ) = u d 1 + iu d 2 , σ d (τ ) = σ rr + iσ rθ are the displacements and tractions at the boundary, respectively.The kernels K 1 (τ, z), K 2 (τ, z) and the Kolosov constant κ d in equations ( 58) are defined as Elastic displacement and stress fields can be determined everywhere in the disc via Kolosov-Muskhelishvili complex potentials φ(z) and ψ(z) as [39] where Re and Im denote real and the imaginary parts, respectively.Integral expressions for the complex potentials φ (z) and ψ (z) were obtained for a circular disc in [40], by evaluating the integrals in equation (58) and using equation (60) [41], to obtain where and z c denotes the centre of the disc, and A ±n are the complex coefficients in the Fourier series expansions for the displacements on the boundary of the disc, as explained below.

Complex combinations for elastic fields on the disc's boundary
The complex Fourier series representation for the displacement at every point τ = Re iθ on the boundary L of the disc is introduced as where A ±n are unknown complex coefficients and functions g ±n (τ ) are defined from equation (62 The relation between Cartesian and polar coordinates, allows for expressing the displacement components at every point τ in the polar coordinate system (r, θ) as so that the final representations for displacements are obtained from equation (63) as The complex Fourier series representation for the tractions at any point τ ∈ L are introduced as where σ d rr and σ d rθ are the radial and tangential components of the tractions at the point τ ∈ L, respectively, and B ±n are the unknown complex coefficients.The expressions for the tractions components can be obtained by separating the real and imaginary parts in equation (68) as The complex coefficients A ±n and B ±n are interrelated as [42] B −1 = 0, To satisfy the condition of inextensibility of the coating, additional relations for the complex coefficients A ±n can be obtained by using equation (99) 2 in [43] and the following relations [25] Re (A 1 ) = 0,

Complex variable formulation for bifurcation
Expression (55) can be derived with respect to the arclength s to yield so that the term S σ in equation ( 57) becomes (for j= σ) If the coating is perfectly connected to the disc or if sliding can occur, the following conditions (displacement continuity in the former case, partial continuity and vanishing of shear stress in the latter) have to be imposed, and equation ( 56) can be rearranged as Note that the term uc θ simplifies in equation ( 75) for all cases, except for the combination of slip contact and dead load, in which case, equation (75) has to be differentiated with respect to θ and the inextensibility condition has to be enforced.A governing sixth-order differential equation is obtained holding for slip contact and dead radial loading: In the following, only the cases pertinent to equation (75) will be explicitely derived, while for the sake of brevity analogous treatment of equation (76) will not be reported.Using complex variables formalism, the bifurcation equation ( 56) can be rewritten by adopting the Fourier series representation introduced in Section 4.2 for the incremental boundary displacement and stress components at a point τ ∈ L. The first three terms on the left-hand side of equation ( 75)have been already derived in [25] and can now be adapted as (where the superscripts 'c' and 'd' have been omitted) so that, using equations (43) reported in [25], equation (77) becomes The last three terms on the left-hand side of equation ( 75) can be rearranged by combining together equations (99) 4−5 and (99) 1 reported in [43], leading to Employing equations (42) reported in [25], equation (79) can be written as Using equation (99) 1 reported in [43] in equations (48), the complex form of S Π introduced by equation (75) becomes for hydrostatic pressure (i.),
• Determination of coefficients A 1 and A 0 • The bifurcation condition, holding for modes of order n ≥ 2: where When A 1−n = 0 the trivial solution is obtained, otherwise, the bifurcation radial load for the coated disc, corresponding to the n-th mode, is: where M = 1 (M = 0) for perfect bonding (for slip contact) at the rod/core interface and ξ = α = 1 for hydrostatic pressure, ξ = 1 and α = 0 for centrally directed load and ξ = α = 0 for dead load.
It is worth noting that equation (88) has been obtained for five load and interface combinations, excluded the sixth case of slip contact plus dead loading, which requires a separate treatment based on equation (76).This treatment is omitted for brevity but leads again to equation (88), which is found to hold true in all cases.Equation (88) shows that, when parameter µ d bR 3 /B tends to zero, the coated disc behaves as a rod subject to the radial load Π.For a given set of material and geometrical parameters (E c , E d , ν d , R, h, b) and varying the mode number n in equation ( 88), different values for the bifurcation load can be analyzed.The critical value corresponds to the integer number n that minimises equation (88), so that from the expressions (86) and (71), the only non-vanishing coefficients are A 1−ncr and A 1+ncr , define the bifurcation mode.For the elastic rod coating the disc, the latter corresponds to the displacement components where the non-vanishing coefficient A 1−n remains arbitrary, while A 0 and Im(A 1 ) can be computed by fixing the rigid-body displacement, as shown in [40,43].The stress components on the boundary of the disc are All displacement and stress fields at every point z within the boundary of the disc may be determined from equations (60), where the complex potentials and their derivatives can be obtained from equation (61) as The coefficients A 0 and A 1 are determined by fixing a rigid-body displacement, in particular, equation (60) 1 assumes the form so that the condition A 0 = 0 is obtained by imposing the displacement to be zero at point z = z c in equation (92), as in [25].Again, requiring the displacement component u θ to be zero at the point τ 0 = R, the following expression for A 1 is obtained: and hence, under the conditions A 0 = 0 and (93), fixing possible rigid-body roto-translations, the displacement field is determined by equations (89) leading to (94) The displacement amplitude is ruled by the non-vanishing coefficient A 1−n which remains arbitrary as in a standard Sturm-Liouville bifurcation problem.

Bifurcation results
The bifurcation radial load, equation (88), evaluated for the three different types of radial forces per unit length (i)-(iii), as listed in Section 1, and for the two models of bonding at the interface, has been normalized through division by the bifurcation load of the 'empty' coating, equation (49), and has been evaluated as a function of the contrast ratio E d /E c between Young's moduli of the coated disc.In addition to the latter parameter, the bifurcation load depends on κ d and on the dimensionless parameter bR 3 /[(1 + ν d )κ d J].The latter has been assumed equal to 1000, while the Kolosov constant has been selected as κ d = 2, corresponding to ν d = 1/4 in plane strain and ν d = 1/3 in plane stress.Bifurcation results are reported in figure 7, for perfect bonding at the disc/coating interface (a) and for frictionless slip (b).Both cases exhibit similar behaviour, in which the critical load increases with the increasing stiffness contrast between the disc and the coating.However, the critical loads are smaller under slip conditions than perfect bonding conditions.The increase in the critical load is accompanied by increase in the wavenumber n of the bifurcation mode, evidenced by the alternance of different colour marking the curves.
The shapes of the bifurcation modes are reported in figure 8 for the case of perfect bonding, from = 2 to n = 7.These have been obtained by choosing the non-vanishing coefficient A 1−n = 1 and fixing the rigid body motion accordingly to equation (93).The contour of the bifurcation modes are highlighted red (blue) where tensile (compressive) tractions in the radial direction prevail.The zones under incremental tension may be expected to detach as a consequence of adhesion failure.For the case of slip contact and hydrostatic pressure loading, the results have been compared with Herrmann and Forrestal [24] and included in the figure as a dashed black line.The solution by Herrmann and Forrestal was obtained under the plane strain assumption and by introducing several mathematical approximations, so that for ν d = 1/4 the comparison shows the correct trend, although results are not superimposed.Extensive analyses performed by us (and not reported for brevity) show that the approximate solution by Herrmann and Forrestal becomes tight to our solution ν d ≥ 0.455.In particular, figure 9 reports the plane strain analysis for ν d = 0.455,0.455,pertaining to the case of slip contact, and shows that the Herrmann and Forrestal approximation is superimposed to our bifurcation curve.pressure, centrally directed (ξ = 1, α = 0) and dead (ξ = 0) load, when slip contact holds at the interface, for ν d = 0.455 under plane strain assumption.The black dashed line corresponds to the approximate solution obtained in [24] for a hydrostatic load distribution.This approximate solution matches very well our bifurcation condition for ν d ≥ 0.455.

Conclusions
The Kolosov-Muskhelishvili complex potential technique has proved to enable the analytical solution of the bifurcation of an elastic disc coated with an inextensible elastic rod.The latter can be fully bonded or in slip contact with the inner disc and is subject to three different types of external radial loads, all uniformly distributed.This new solution reveals that the presence of the inner disc may lead to the predominance of high-wavenumber modes.It also emphasizes the importance of detachment at the interface between disc and coating, as well as, of proper modelling of how the external load acts on the structure during its incremental deformation.Applications of the results could be valuable in the development of various coating technologies and in the understanding of growth of plants and fruits.

Figure 2 :
Figure2: A circular rubber disc (Agilus30 ™ , 44 mm in diameter) is connected to a stiffer photopolymer (VeroYellow ™ , 1 mm thick) coating through radial ligaments.The latter can transmit axial but only negligible transverse forces.Therefore, a model of slip interface can be adopted to model the disc (manufactured with a Stratasys J750 3D printer) when radially loaded on its external surface.The piece is bioinspired by the joining through ligaments of the brain to the cranial vault.

Figure 4 :
Figure4: (a) The deformation g (x 0 ) maps a rod from its reference configuration (characterized by points x 0 , arclength parameter s 0 , unit tangent t 0 and normal n 0 ) to the deformed configuration (characterised by points x, arclength parameter s, unit tangent t and normal n) through the displacement u(x 0 ).(b): The external and internal forces acting on the current configuration and their counterparts in the reference configuration.

Figure 5 :
Figure5: The deformation of an annular rod, where points x 0 of the reference configuration are displaced to points x of the current configuration.The radial load Π, uniform in the reference configuration, is assumed to follow three different models under incremental deformation: (a) hydrostatic pressure, where the force resultant P over an elementary arc ds follows the normal m in the deformed configuration, (b) centrally directed load, where the resultant P is pointing towards the centre of the rod in the undeformed configuration, (c) dead load, where the original direction of the resultant P remains unaltered by deformation.

Figure 7 :
Figure 7: Dimensionless radial load for bifurcation Π(n) as a function of theratio E d /E c for hydrostatic (ξ = α = 1) pressure, centrally directed (ξ = 1, α = 0) and dead (ξ = 0) load.The different colours correspond to the critical bifurcation number ncr.(a) Perfect bonding at the disc/coating interface.(b) Slip contact, where the black dashed line corresponds to the approximate solution obtained in [24] for a hydrostatic load distribution.

Figure 8 :
Figure8: Bifurcation modes for the coated disc at increasing wavenumber n.Perfect bonding is assumed at the disc/coating interface.The parts highlighted red (blue) show zones at the interface where tensile (compressive) tractions in the radial direction prevail.

(Figure 9 :
Figure9: Dimensionless radial load for bifurcation Π(n) as a functionof the ratio E d /E c for hydrostatic (ξ = α = 1) pressure, centrally directed (ξ = 1, α = 0) and dead (ξ = 0) load, when slip contact holds at the interface, for ν d = 0.455 under plane strain assumption.The black dashed line corresponds to the approximate solution obtained in[24] for a hydrostatic load distribution.This approximate solution matches very well our bifurcation condition for ν d ≥ 0.455.