Stability of heterogeneous beams with three supports through Green functions

The present paper is devoted to the issue how the critical load of some heterogeneous beams with three supports can be determined by using Green functions. The stability problems of these beams are equivalent to three-point boundary value problems, paired with homogeneous boundary conditions. If the Green functions of these boundary value problems are known, the eigenvalue problems that provide the critical load can be transformed into eigenvalue problems governed by homogeneous Fredholm integral equations. The later eigenvalue problems can be reduced to algebraic eigenvalue problems which then can be solved numerically with effective algorithms.

propose a new standard as relevant differences were noted in comparison with the current ones. Work [13] by Vaz and Silva is based on the perturbation method. When there are terminal axial forces, the related set of first-order differential equations with boundary conditions determine a two-point boundary value problem that is solved with the mentioned technique. The finite element method is used by Cai et al. [14]. A mixed complementary energy functional is written and it is shown that the global maximum of the complementary energy is related to a stable equilibrium, while local maxima correspond to unstable configurations. A variational iterational method is used to find the buckling loads of uniform/variable cross-sections in [15] by Coskun and Atay. This method is proven to be effective in solving multiple classes of linear and non-linear differential equations, including boundary value problems as well. Functionally graded non-uniform beams are in the spotlight in [16] by Singh and Li. A new low dimensional mathematical model is presented. The solution method uses the transcendental functions of a uniform column, instead of polynomials. Non-homogeneous columns are replaced by (continuous) systems with constant properties leading to a lower dimensional transcendental eigenvalue problem.
Batista [17] uses a variational principle to find the equilibrium configuration of compressed elastic rods. The related bifurcation diagrams are created with the Jacobi test. In this way, the whole phase plane is assessed to reveal the typical stability regions. Bigoni et al. [18] focus on the stability analysis of equilibrium configurations. A novel phenomenon-namely the elastic tensile buckling of rods, when an internal slider constraint is present-is also addressed. Levyakov and Kuznetsov [19] consider inextensible rods. The second variation of the total potential energy is given using the eigensolutions of a Sturm-Liouville problem. The novelty of the proposed method is the incorporation of isoparametric constraints. Ruta et al. [20] tackle the flexural-torsional buckling of thinwalled hyperelastic beams with arbitrary cross-section. The bifurcation points are found using a static perturbation method. The subsequent study [21] is devoted to dynamic investigations of the issue.
Paper [22] by Harvey et al. is about how a geometrically imperfect beam can buckle. The ideal beam shapes are perturbed with their second-mode eigenshape by means of 3D printing and these specimen are experimentally investigated. Virgin [23] uses additive manufacturing to print beams with high-precision and test the behavior of these. Article [24] by Madah and Amir is about optimization against buckling when initial imperfections are present.
Furthermore, as for the concept of the Green function, the first publication of Green himself [25] dates back to 1828. This book presents the Green theorem, introduces the concept of the Green function and applies it to electrostatic problems, governed by partial differential equations. The Green function for two-point boundary value problems governed by ordinary differential equations (ODE) was established by Bocher [26]. The first book, written by Ince [27], that systematically contains the concept of the Green function came out in 1926. The books of Collatz [28,29] contain the definition of the Green function for two-point boundary value problems governed by ordinary linear differential equations and clarify their most important properties. In addition a number of Green functions are presented in closed-form. The concept of the Green function considering twopoint boundary value problems was generalized for a class of ordinary differential equation systems in [30,31] by Obadovics. The results of [31] were generalized for degenerated ordinary differential equation systems by Szeidl [32,33]. Existence proof for some three-point boundary value problems related to third-order nonlinear differential equations is given in paper [34] by Murty and Kumar using Green functions. For some three-point boundary value problems governed by linear ordinary differential equations of order two the corresponding Green functions are presented in the paper of Zhao [35]. The definition of the Green function for three-point boundary value problems with a solution algorithm is detailed in [36] by Szeidl and Kiss. In addition the free vibration issue of pinned-roller supported-pinned beams is solved numerically. Stability problems are, however, not investigated in that paper at all.
There are also recent works which present a solution method that incorporates the Green function for beam buckling. Wolfe [37] presents buckling investigations through multiple theories: global bifurcation theory, direct variational calculus, phase plane analysis and singular perturbation theory. These methods are compared and it is shown how these compliment each other. The bifurcation theory uses the implication of the Green function for a classical compressed column with two supports. Anghel and Mares [38] investigate the buckling of Euler-Bernoulli beams. One of the possible solution methods is rest on an integral equation with the kernel being the Green function. This function is found numerically. Köylüoglu et al. [39] dedicate their study to the stability of columns with material and geometrical uncertainties. The solution is based on integral equations and optimization techniques to set up conservative (practically, lower) bounds to the buckling load when there are two supports. Totry et al.
[40] focus on heterogeneous columns, using a direct functional perturbation method on the buckling equation. The analytical solution of eigenvalue functional differential equations is given. The method is compared with the Ritz approximation. Here, the functional derivatives act as Green functions. Khan and Al-Hayani [41] use the Adomian decomposition method together with the Green function technique to solve the non-linear buckling equations of a column. This method is proven to be efficient and, as an advantage, it does not require the use of the perturbation theory. Huang and Li [42] present an efficient method, by using Fredholm integral equations instead of differential equations, that is capable of dealing with the stability of axially graded material columns with variable cross-section and two supports. The proposed method can tackle various classes of non-constant flexural rigidities.
Based on the above literature research, in this article, the buckling problem of compressed beams are solved. The material is linearly elastic, isotropic and the material distribution can change over the crosssection-it is called cross-sectional heterogeneity. The stability equations are given for the selected, say, three-point boundary value problems. Utilizing the definition given in [36], these are then replaced by Fredholm integral equations, whose kernel is obtained from the Green function. Furthermore, a definition of this Green function for three-point boundary value problems with homogeneous boundary conditions is also provided. The Green function itself is also given in closed-form. Numerical solutions of the integral equations are given based on a boundary element technique-algebraic equations are introduced in this way. It is found that the position of the middle support and material distribution have significant effects on the linear buckling loads.

Governing equations
Let us consider three beams, each with three supports as shown in Fig. 1. The (E-weighted) centroidal axis of these beams coincide with axis x . The two further mutually perpendicular coordinates are ŷ,ẑ . Along the points of the (E-weighted) centroidal axis, the cross-sectional E-weighted first moment vanishes. The beams have uniform cross-section, symmetric to plane xẑ and are made of linearly elastic, isotropic materials. The modulus of elasticity E satisfies the condition E(ŷ,ẑ) = E(−ŷ,ẑ) over the points of the cross-section A, which means that (a) E is an even function of ŷ and (b) it is independent of the coordinate x . Thus, we speak about cross-sectional heterogeneity. The length of the beam is L, and the position of the middle support is given by b .
Equilibrium problems of beams with cross-sectional heterogeneity-the axial force N is zero-are governed by [43] the ordinary differential equation: where ŵ(x) is the vertical displacement of the centroidal axis, f z (x) is the intensity of the distributed load acting on the centroidal axis and I ey is If E is constant throughout, i.e., the beam is homogeneous, then in which I is the moment of inertia to ŷ.
It is practical to introduce some dimensionless variables defined by relations where ̂ is a coordinate on the axis x . Rewriting now equation (2.2) with dimensionless quantities yields This ordinary differential equation (ODE) is associated with the boundary-and continuity conditions presented in Table 1. A fixed support is abbreviated by F, a roller by r and a pin by P-in this way, e.g., FrF is a fixed-roller supported-fixed beam.
Here and now on, a simplified notation for the derivatives is used, that is The general solution of the homogeneous differential equation is the sum (2.7) d n w dx n = w (n) .
If the Green function G(x, ) for the boundary value problem determined by ODE (2.6) and the boundary and continuity conditions is known, the solution for the dimensionless deflection w can be given by the integral

Stability problem
Equilibrium problems of uniform heterogeneous beams subjected to an axial force are governed by the differential equation where N is an axial concentrated force and the sign of N is positive for compression and is negative for tension. When dealing with the stability issue, the force is compressive and the intensity f z is zero.
The three support arrangements actually make three eigenvalue problems-with N being the eigenvalue-determined by the differential equation and the related boundary-and continuity conditions. If −N w (2) is written for f z in (2.9) it is found that where since the Green function should satisfy the boundary conditions. Hence,

Deriving this equation with respect to x yields
After introducing new variables: a homogeneous Fredholm integral equation is found: In this way, the eigenvalue problems determined by differential equation (2.11) and the homogeneous boundary-and continuity conditions are replaced by eigenvalue problems governed by homogeneous Fredholm integral equations as of [29].

Remark 1 Differential equation (2.11) can be rewritten in the form
The ] are comparative functions if (a) they satisfy the boundary and continuity conditions given in Table 1 and The exact solution w(x) is obviously a comparative function.

Remark 2
The integrals taken on the set of the comparison functions u(x), v(x) are products defined on the differential operators K and M. The eigenvalue problems defined by differential equation (2.11) and the homogeneous boundary and continuity conditions presented in Table 1 are said to be self-adjoint if the products (2.19) are commutative, i.e., it holds that Conditions (2.20) are called conditions of selfadjointness. The differential operator K(w) is also comparative if the product (u, v) K under the corresponding boundary and continuity conditions is commutative. It can be checked by performing partial integrations and taking the corresponding boundary and continuity conditions into account that the differential operator K(w) and the eigenvalue problems mentioned are all self-adjoint.

Definition
Let us consider the inhomogeneous ODE where the differential operator of order 2k is defined by the following equation in which k is a natural number, the functions p n (x) and r(x) are continuous and It is assumed that the inhomogeneous differential equation (3.1) is paired with homogeneous boundary and continuity conditions given by the following equations: The notations ◻ I and ◻ II refer successively to the intervals [0, b] and [b, ] : y I and y II are the solutions to the differential equation in the interval I and II while nrI , nrI , nrII and nrII are arbitrary constants.
Solution of the three-point boundary value problem (3.1), (3.2) and (3.3) is sought as where the Green function is defined by the following formulas and properties [36]: In addition it is 2k times differentiable with respect to x and the derivatives are also continuous functions of x and in the triangular its derivatives should be continuous for x = : For a fixed ∈ [0, ] the product G(x, ) as a function of x ( x ≠ ) should satisfy the homogeneous differential equation (3.9) 5. The product G(x, ) as a function of x should satisfy both the boundary conditions and the continuity conditions The above continuity conditions should be satisfied by the function pairs Remark 3 It can be proven [36] by utilizing the previous definition that integral (3.4) satisfies Eq. (3.3).

Remark 4 If the boundary value problem defined by (3.1) and (3.3) is self-adjoint then the Green function is cross-symmetric [36]:
In what follows, the related Green functions are presented for the three selected support arrangements. The calculations are detailed for FrF beams. As regards FrP, PrP beams, solely the final formulae are given.

Calculation of the Green function if
The coefficients a mI ( ), b mI ( ) and c mI ( ) are unknown functions, w m (x) is given by (2.8b).  The previous linear equations in identical matrix form are After solving the system, it is found that As regards G 2I (x, ) the result is Then which is the Green function for fixed-fixed beamssee Table 8.1 in [33]. (3.22b)

Calculation of the Green function if
where the coefficients a mII ( ), b mII ( ) and c mII ( ) are again unknown functions. Continuity and discontinuity conditions (3.12) lead again to equation system (3.18) in which now the coefficients b mII ( ) , m = 1, 2, 3, 4 are the unknowns. Consequently b mII ( ) = b mI ( ) . Utilizing again the boundary and continuity conditions given in Table 1, equations can be obtained for the eight remaining unknown coefficients a mII ( ) and c mII ( ).  Since c 1II = c 2II = 0 , the final equation system has the following form: After solving the above equation system, the results for G 1II (x, ) and G 2II (x, ) are and (3.26e) where It is emphasised that the definition for the Green function of three-point boundary value problems is a constructive one since the calculations are based on those five properties which constitute the definition. According to Remark 2, the boundary value problem is self-adjoint. Consequently, the Green function, given by equations (3.5), (3.22) and (3.28), should satisfy symmetry condition (3.14). It can be proven by time consuming paper and pencil calculations that this condition is indeed fulfilled.

FrP beams
Repeating the calculation steps detailed in Subsect. 3.2 for FrP beams, the four corresponding elements for the Green function are and where The Green function given by Eqs. (3.5), (3.30) and (3.31) for FrP beams satisfies symmetry condition (3.14). Its graph is shown in Fig. 3 if L = 100 mm, b = 50 mm and ̂= 75 mm.

Remark 9
Substitute for b in (3.30a). We obtain For G 2II (x, ) | |b=0 we get which is the Green function for fixed-pinned beams [33].

PrP beams
The equations that constitute the Green function for PrP beams are taken from [36]: where The Green function (3.5), (3.35) and (3.36) valid for pinned-pinned-pinned beams also satisfies symmetry condition (3.14). Figure 4 shows the Green function if L = 100 mm, b = 50 mm and ̂= 75 mm.

Remark 10
Substituting for b in (3.36a) reults in (3.33) which is the Green function for pinned-fixed beams. Setting b to 0 in (3.36b) yields (3.34) which is the Green function for fixed-pinned beams.

Solution procedures
There are multiple ways to find the critical load. One can solve the eigenvalue problem determined by the homogeneous Fredholm integral equation (2.11) numerically by applying the boundary element technique or one can set up the characteristic equationsthese are nonlinear equations for the unknown critical load-which can also be solved numerically. The use of the boundary element technique is preferred primarily hereinafter and the characteristic equation will serve as a possibility to check the results. As regards the boundary element technique, the solution steps are detailed in book [33]-the reader is referred to Subsection 8.15.2. It is noted that the interval [0, 1] was divided into 40 elements and the resulted algebraic eigenvalue problem was solved by using subroutine DGVLRG of the International Mathematical Science Library [44].  The kernel in Eq. (2.11) assumes the following form where Accordingly, determination of the kernel K(x, ) requires the calculation of second derivatives.

Computational results for FrF beams
To calculate the elements of the kernel function, Eqs. (3.22), (3.28) are recalled, and thus, it is found that x , x .
(4.2a)  Based on the above, Table 2 is a summary of the computed results. The first column contains the dimensionless parameter b which identifies the location of the middle roller support. For symmetry reasons it is sufficient to consider its value in the interval [0, 0.5].
The second column contains the critical value for the dimensionless compressive force, more precisely, the quantity √ N crit ∕ against the discrete values of b. A polynomial of degree five was fitted onto the point pairs taken from the first two columns that is The third column contains the approximations computed using polynomial numbers in columns two and three are the same with the accuracy of four to five digits. Figure 6 represents the computational results. The discrete points are depicted by diamonds, while the polynomial itself is drawn by a continuous line.
Note that for b = 0 the beam behaves as if it were a fixed-fixed beam for which √ N(x) ∕ = 2.00 is the exact value. It is obvious that the critical force reaches its maximum if b = 0.5.    Table 3 Numerical results for the critical load when the beam is FrP With the kernel in hand, the second and fourth columns in Table 3 contain the values of √ N crit ∕ as a function of b. Figure 7 represents the computational results. The optimum location for the middle roller support is at b ≈ 0.635 . The fitting polynomials are given by If b = 0 , the beam behaves as if it were fixed-pinned for which √ N(x) ∕ = 1.4303 is the exact value. However, for b = 1.0 the beam behaves as fixed-fixed beam. Then √ N(x) ∕ = 2.000 . Columns 3 and 6 in Table 3 contain the values of the polynomials (4.6). The critical force reaches its maximum if b = 0.65.

Remark 11
(4.7a) With regard to the computational results, the second column in Table 4 shows the computational results, i.e., the critical values of √ N ∕ for twenty discrete points of x = b . For symmetry reasons, in the same manner as for FrF beams, it is sufficient to consider the results for the first half of the beam only.
A polynomial was fitted onto the point pairs in the first two columns. (4.7d) The third column contains the approximations computed using polynomial √ N(x) ∕ . Note that the numbers in columns two and three are the same with the accuracy of four to five digits. Figure 8 represents the computational results: the discrete points are depicted again by diamonds while the polynomial is drawn by a continuous line.
It is obvious that the critical force reaches its maximum if b = 0.5 . The PrP beam behaves as if it were a fixed-pinned beam if b tends to zero. Then √ N crit ∕ = 1.4303 is the exact result . If b = 0.5 the first half of the beam behaves as if it were a pinnedpinned beam for which √ N crit ∕ = 2.0.

Remark 13
According to Eq. (A.7) the characteristic equation which provides the eigenvalues, i.e., the critical force has the following form: This non-linear equation for p∕ . The solutions obtained are the same with five to six digits accuracy as those computed by solving the algebraic eigenvalue problem that belongs to the Fredholm integral equation (2.17).

Example for heterogeneous beam
It is now presented through a simple example how the material heterogeneity can influence the lowest buckling load of PrP beams. For the demonstration,  a simple two-layered rectangular section is selected, as shown in Fig. 9. The width is a, the total height is h, while the height of the lower part is h 1 = kh; k ∈ [0;1] . The typical Young modulus values are noted by E 1 and E 2 = dE 1 ; d ∈ ℝ.
To tackle the stability problem, first the relative position of the (E-weighted) centroid C e is neededit is identified by the coordinate z e (k, d) . It can be found based on the definition of Eq. (2.1). After that, Eq. (2.3) should be utilized. With Eq. (2.10) 2 and the results of Table 4 it is possible to evaluate the critical loads.
In Fig. 10, the N critical loads for a homogeneous ( d = 0 ) cross-section with two supports b = 10 −5 ≃ 0 and a heterogeneous cross-section are compared when d = 1∕3 . Otherwise, the length of the beams are identical. Computational results are given in terms of the support position b and the layer height parameter k. The plotted relationships are non-linear and the variation of the critical loads is absolutely relevant in this aspect.

Conclusions
The present paper determines the Green functions for three three-point boundary value problems, namely, for those boundary value problems which describe the mechanical behavior of FrF, FrP and PrP beams with cross-sectional heterogeneity. Making use of the Green functions, the linear stability problem of these beams is transformed into eigenvalues problem governed by the homogeneous Fredholm integral equation: Then the eigenvalue problems are replaced by an algebraic eigenvalue problem using the boundary element technique. Its numerical solution is compared to the solutions obtained by numerically solving the corresponding characteristic equations presented in the Appendix. Polynomial approximations for the critical force with a very good accuracy have also been derived.
Funding Open access funding provided by University of Miskolc.

Declarations
Conflict of interest Not applicable.