Optimization of the dimensions of a tapered thin-walled cantilever, aimed at reaching a maximum critical buckling load

Analyses of the thin-walled elements of bar structures have been conducted by many authors, as evidenced by the extensive literature on this subject. Andrade and Camotim [1] derived differential equations for the analysis of the lateral torsional buckling (LTB) of monosymmetric thin-walled beams with a tapered cross section. The derived equations were solved using the Rayleigh–Ritz method. The LTB of simply supported I-beams and cantilever I-beams, both with a tapered web, was analysed. The results were compared with those reported by other authors. The LTB problem was also the subject of study by Benyamina et al. [2], who derived nonlinear equilibrium equations and used the Ritz method to solve them. Bisymmetric thin-walled beams with tapered webs were analysed. The results were compared with those obtained using the finite element method (FEM). Andrade et al. [3] presented a one-dimensional model enabling the analysis of the LTB of nonprismatic thin-walled beams with a monosymmetric open cross section. I-shaped thin-walled beams and cantilever beams, both with linearly tapered webs, were analysed. In order to verify the model, the numerical results obtained using this model were compared with those obtained by FEM. Asgarian et al. [4] presented a nonprismatic thin-walled beam model and used it to analyse the LTB stability of nonprismatic thin-walled beams with any cross sections and various support schemes. The differential equations with variable coefficients, derived using Hamilton’s principle, were solved using the power series method. In order to verify the efficacy of the method used, the results were compared with those reported by other authors and with those obtained using FEM. Branford [5] used the latter method to analyse the LTB of monosymmetric nonprismatic I-beams under any load and any support conditions. The obtained results were compared to those reported in the literature on the subject. Challamel et al. [6] analysed the LTB of cantilever beams with a linearly tapered web. An analytical solution to the problem was obtained in the form of hypergeometric functions. Also, a solution to the cantilever shape optimization problem was presented, assuming the maximum warping resistance of the cantilever as the optimization criterion. The effect of shear strains on the buckling of nonprismatic thin-walled beams with a bisymmetric cross section was studied by Osmani and Meftah [7]. Beams with an open and closed cross-section, loaded with a combination of forces causing bending and axial deformations, were analysed. The Ritz method was used to solve the nonlinear equilibrium equations. The critical load determined using the proposed method was shown to agree with the results yielded by FEM (ABAQUS). Raftoyiannis and Adamakos [8] developed a numerical procedure for determining the critical load causing the LTB of I-beams with nonprismatic webs. The obtained results were compared with those reported in the literature on the subject. Rezaiee-Pajand et al. [9] studied the LTB of a nonprismatic thin-walled beam at variable material parameters (the elasticity modulus and the rigidity modulus) in both the longitudinal and transverse directions. Beams with a symmetric double-tee cross-section and a linearly variable depth of the web were analysed. The sought total displacement functions were approximated with sine trigonometric series. Soltani [10] presented a method of analysing the LTB of nonprismatic composite beams with a symmetric double-tee cross section. The problem was solved using the Ritz method. In order to demonstrate the accuracy of the proposed method, the results were compared with those obtained by FEM (ANSYS) in which shell elements were used for the approximation. Soltani and Asgarian [11,12] used the FEM and the differential quadrature method to analyse the LTB of bisymmetric nonprismatic beams with variable (along the length) material parameters. Beams with various load diagrams were analysed. Energy methods were used to derive equilibrium equations. The effect of factors such as material property variation along the bar length, bending moment gradient variation, or the transverse load application level on LTB was studied. The LTB of beams with variable material parameters along the length was investigated by Soltani et al. [13]. Pinned nonprismatic I-beams with material parameters (continuously) varying along the length were analysed using FEM. Soltani et al. [14] used the finite difference method to solve the LTB problem for thin-walled beams with any load diagram. In order to verify the applied method, the results were compared with those obtained by FEM, in which shell elements were used for approximation. Zhang and Tong [15] analysed the LTB of I-beams with a nonprismatic web. Using Vlasov’s assumptions and variational methods, they derived equations describing the behaviour of this type of beam. It was shown, among other things, that the FEM using prismatic beam elements for approximation may yield results significantly different from the results obtained using other methods. El-Mahdy and El-Saadawy [16] created a three-dimensional thin-walled beam model in ANSYS. They used this model to analyse the effect of the variable ratio of the compression flange dimensions to the tension flange dimensions (at a fixed quantity of material used) on the optimal performance characteristics of the system, i.e. on the beam’s load-bearing capacity considering the LTB.

In the present study, displacement differential equations describing the global buckling problem for a thin-walled nonprismatic bar with a bisymmetric cross section, subjected to a uniformly distributed load, are derived using the minimal elastic strain energy principle. The recurrence algorithm described in earlier papers by the author(s) [17,18,19,20,21,22,23,24,25] is used to solve the system of displacement differential equations with variable coefficients. In this algorithm, Chebyshev series of the first kind are employed to approximate the generalized displacement functions. A tapered cantilever with a bisymmetric cross-section, subjected to a uniformly distributed load, is analysed. The load can be applied to the I-beam’s upper flange, its lower flange, or it can act in the web’s centre. The determined critical load values are compared with those obtained using FEM and the commercial Sofistik [26] software. The obtained results show that among the linearly tapered cantilevers with a bisymmetric double-tee cross section, cantilevers with the highest web and flanges taper coefficients (the free end having the smallest possible dimensions, while the fixed end has the largest possible dimensions) are the optimal ones (in their case, at a constant mass, the critical load causing the LTB reaches its maximum value. This finding can be surprising as it indicates that the cantilever’s “stronger” part (at the fixed end) contributes considerably more to an increase in critical load than its “weaker” part (at the free end).

A model presented previously [23,25] was used in this study to describe a thin-walled beam with an open cross-section. The model was based on the momentless theory of shells and Vlasov’s assumptions, and it is a generalization of Wilde’s model [27]. Here, only the essential equations describing the model are quoted, specifically for the case when the system’s cross sections are bisymmetric. The detailed derivation of the general equations for any cross-section can be found from the previous study [23]. Figure 1 shows the coordinate systems used in the description of the considered model.

Figure 1

Coordinate system and local basis vectors for a nonprismatic thin-walled beam.

Under the adopted assumptions, the displacement of any point of the thin-walled bar is given by the following formula: (1) u = U x e x + U y e y + U z e z = u x ( x , s ) e x + [ u y ( x ) − z ( x , s ) θ ( x ) ] e y + [ u z ( x ) + y ( x , s ) θ ( x ) ] e z , \begin{array}{c}{\boldsymbol{u}}={U}_{x}{{\boldsymbol{e}}}_{x}+{U}_{y}{{\boldsymbol{e}}}_{y}+{U}_{z}{{\boldsymbol{e}}}_{z}\\ \hspace{1em}={u}_{x}(x,s){{\boldsymbol{e}}}_{x}+{[}{u}_{y}(x)-z(x,s)\theta (x)]{{\boldsymbol{e}}}_{y}+{[}{u}_{z}(x)+y(x,s)\theta (x)]\hspace{.25em}{{\boldsymbol{e}}}_{z},\end{array} where θ ( x ) \theta (x) is the angle of rotation relative to the axis defined by the points with coordinates ( C y ( x ) = 0 , C z ( x ) = 0 ) ({C}_{y}(x)=0,\hspace{.25em}{C}_{z}(x)=0) specifying the shear centre; u y ( x ) , u z ( x ) {u}_{y}(x),{u}_{z}(x) are the cross-section displacement components in the shear centre, measured along the axis of the global coordinate system, u x ( x , s ) {u}_{x}(x,s) is a measure of cross-sectional displacement along the x-axis, which takes into account warpings relative to the plane x = const x=\text{const} (Figure 2).

Figure 2

Description of displacements for the cross-section.

As shown in previous studies [23,27], the use of Vlasov’s assumptions ( γ x s = 0 {\gamma }_{xs}=0 ) leads to the following formula for displacements u x {u}_{x} : (2) u x = u o x − y u y , x − z u z , x − ω θ , x , {u}_{x}={u}_{ox}-y{u}_{y,x}-z{u}_{z,x}-\omega {\theta }_{,x}, where u o x = u o x ( x ) {u}_{ox}={u}_{ox}(x) is an integration constant independent of s s , and ω = ω ( x , s ) \omega =\omega (x,s) is a sectorial coordinate described by the relationship (3) ω ( x , s ) = ∫ [ − z y , s + y z , s ] d s . \omega (x,s)=\int {[}-z{y}_{,s}+y{z}_{,s}]\text{d}s.

The linear and nonlinear components of the strain tensor are defined by the respective formulas (4) γ α β = 1 2 ( a α ∘ u , β + a β ∘ u , α ) , {\gamma }_{\alpha \beta }=\frac{1}{2}({{\bf{a}}}_{\alpha }\circ {{\bf{u}}}_{,\beta }+{{\bf{a}}}_{\beta }\circ {{\bf{u}}}_{,\alpha }), (5) γ x x n l = 1 2 ( u , x ∘ u , x ) , γ x s n l = 1 2 ( u , x ∘ u , s ) . {\gamma }_{xx}^{nl}=\frac{1}{2}({{\bf{u}}}_{,x}\circ {{\bf{u}}}_{,x}),\hspace{1em}{\gamma }_{xs}^{nl}=\frac{1}{2}({{\bf{u}}}_{,x}\circ {{\bf{u}}}_{,s}).

Substituting the displacement function u ( x , s ) {\bf{u}}(x,\hspace{.25em}s) , defined in formula (1), into formulas (4) and (5), the appropriate derivatives were calculated and laborious transformations were carried out and the following formulas were obtained for:

• Linear strains: (6) γ x x = u o x , x − y u y , x x − z u z , x x − ω θ , x x − ψ θ , x , {\gamma }_{xx}={u}_{ox,x}-y{u}_{y,xx}-z{u}_{z,xx}-\omega {\theta }_{,xx}-\psi {\theta }_{,x}, where (7) ψ = ω , x + y , x z − z , x y . \psi ={\omega }_{,x}+{y}_{,x}z-{z}_{,x}y.

• Nonlinear strains:

(8) γ x x n l = 1 2 ( α x , 1 x , 1 u o x , x 2 + α y , 1 y , 1 u y , x 2 + α y , 2 y , 2 u y , x x 2 + α z , 1 z , 1 u z , x 2 + α z , 2 z , 2 u z , x x 2 + α θ θ θ 2 + α θ , 1 θ , 1 θ , x 2 + α θ , 2 θ , 2 θ , x x 2 ) + α x , 1 y , 1 u o x , x u y , x + α x , 1 y , 2 u o x , x u y , x x + α x , 1 z , 1 u o x , x u z , x + α x , 1 z , 2 u o x , x u z , x x + α x , 1 θ , 1 u o x , x θ , x + α x , 1 θ , 2 u o x , x θ , x x + α y , 1 z , 1 u y , x u z , x + α y , 2 z , 2 u y , x x u z , x x + α y , 1 y , 2 u y , x u y , x x + α y , 1 z , 2 u y , x u z , x x + α z , 1 y , 2 u z , x u y , x x + α z , 1 z , 2 u z , x u z , x x + α y , 1 θ u y , x θ + α y , 1 θ , 1 u y , x θ , x + α y , 1 θ , 2 u y , x θ , x x + α y , 2 θ , 1 u y , x x θ , x + α y , 2 θ , 2 u y , x x θ , x x + α z , 1 θ u z , x θ + α z , 1 θ , 1 u z , x θ , x + α z , 1 θ , 2 u z , x θ , x x + α z , 2 θ , 1 u z , x x θ , x + α z , 2 θ , 2 u z , x x θ , x x + α θ θ , 1 θ θ , x + α θ , 1 θ , 2 θ , x θ , x x , \begin{array}{c}{\gamma }_{xx}^{nl}=\frac{1}{2}({\alpha }_{x,1}^{x,1}{u}_{ox,x}^{2}+{\alpha }_{y,1}^{y,1}{u}_{y,x}^{2}+{\alpha }_{y,2}^{y,2}{u}_{y,xx}^{2}+{\alpha }_{z,1}^{z,1}{u}_{z,x}^{2}+{\alpha }_{z,2}^{z,2}{u}_{z,xx}^{2}+{\alpha }_{\theta }^{\theta }{\theta }^{2}+{\alpha }_{\theta ,1}^{\theta ,1}{\theta }_{,x}^{2}+{\alpha }_{\theta ,2}^{\theta ,2}{\theta }_{,xx}^{2})+{\alpha }_{x,1}^{y,1}{u}_{ox,x}{u}_{y,x}+{\alpha }_{x,1}^{y,2}{u}_{ox,x}{u}_{y,xx}+{\alpha }_{x,1}^{z,1}{u}_{ox,x}{u}_{z,x}+{\alpha }_{x,1}^{z,2}{u}_{ox,x}{u}_{z,xx}+{\alpha }_{x,1}^{\theta ,1}{u}_{ox,x}{\theta }_{,x}+{\alpha }_{x,1}^{\theta ,2}{u}_{ox,x}{\theta }_{,xx}+{\alpha }_{y,1}^{z,1}{u}_{y,x}{u}_{z,x}+{\alpha }_{y,2}^{z,2}{u}_{y,xx}{u}_{z,xx}+{\alpha }_{y,1}^{y,2}{u}_{y,x}{u}_{y,xx}+{\alpha }_{y,1}^{z,2}{u}_{y,x}{u}_{z,xx}+{\alpha }_{z,1}^{y,2}{u}_{z,x}{u}_{y,xx}+{\alpha }_{z,1}^{z,2}{u}_{z,x}{u}_{z,xx}+{\alpha }_{y,1}^{\theta }{u}_{y,x}\theta +{\alpha }_{y,1}^{\theta ,1}{u}_{y,x}{\theta }_{,x}+{\alpha }_{y,1}^{\theta ,2}{u}_{y,x}{\theta }_{,xx}+{\alpha }_{y,2}^{\theta ,1}{u}_{y,xx}{\theta }_{,x}+{\alpha }_{y,2}^{\theta ,2}{u}_{y,xx}{\theta }_{,xx}+{\alpha }_{z,1}^{\theta }{u}_{z,x}\theta +{\alpha }_{z,1}^{\theta ,1}{u}_{z,x}{\theta }_{,x}+{\alpha }_{z,1}^{\theta ,2}{u}_{z,x}{\theta }_{,xx}+{\alpha }_{z,2}^{\theta ,1}{u}_{z,xx}{\theta }_{,x}+{\alpha }_{z,2}^{\theta ,2}{u}_{z,xx}{\theta }_{,xx}+{\alpha }_{\theta }^{\theta ,1}\theta {\theta }_{,x}+{\alpha }_{\theta ,1}^{\theta ,2}{\theta }_{,x}{\theta }_{,xx},\end{array}

(9) γ x s n l = 1 2 ( ζ y , 1 y , 1 u y , x 2 + ζ z , 1 z , 1 u z , x 2 + ζ θ θ θ 2 + ζ θ , 1 θ , 1 θ , x 2 ) + ζ x , 1 y , 1 u o x , x u y , x + ζ x , 1 z , 1 u o x , x u z , x + ζ x , 1 θ , 1 u o x , x θ , x + ζ y , 1 y , 2 u y , x u y , x x + ζ y , 1 z , 2 u y , x u z , x x + ζ y , 1 z , 1 u y , x u z , x + ζ y , 2 z , 1 u y , x x u z , x + ζ z , 1 z , 2 u z , x u z , x x + ζ y , 1 θ u y , x θ + ζ y , 1 θ , 1 u y , x θ , x + ζ y , 1 θ , 2 u y , x θ , x x + ζ y , 2 θ , 1 u y , x x θ , x + ζ z , 1 θ u z , x θ + ζ z , 1 θ , 1 u z , x θ , x + ζ z , 1 θ , 2 u z , x θ , x x + ζ z , 2 θ , 1 u z , x x θ , x + ζ θ θ , 1 θ θ , x + ζ θ , 1 θ , 2 θ , x θ , x x . \begin{array}{c}{\gamma }_{xs}^{nl}=\frac{1}{2}({\zeta }_{y,1}^{y,1}{u}_{y,x}^{2}+{\zeta }_{z,1}^{z,1}{u}_{z,x}^{2}+{\zeta }_{\theta }^{\theta }{\theta }^{2}+{\zeta }_{\theta ,1}^{\theta ,1}{\theta }_{,x}^{2})+{\zeta }_{x,1}^{y,1}{u}_{ox,x}{u}_{y,x}+{\zeta }_{x,1}^{z,1}{u}_{ox,x}{u}_{z,x}+{\zeta }_{x,1}^{\theta ,1}{u}_{ox,x}{\theta }_{,x}+{\zeta }_{y,1}^{y,2}{u}_{y,x}{u}_{y,xx}+{\zeta }_{y,1}^{z,2}{u}_{y,x}{u}_{z,xx}+{\zeta }_{y,1}^{z,1}{u}_{y,x}{u}_{z,x}+{\zeta }_{y,2}^{z,1}{u}_{y,xx}{u}_{z,x}+{\zeta }_{z,1}^{z,2}{u}_{z,x}{u}_{z,xx}+{\zeta }_{y,1}^{\theta }{u}_{y,x}\theta +{\zeta }_{y,1}^{\theta ,1}{u}_{y,x}{\theta }_{,x}+{\zeta }_{y,1}^{\theta ,2}{u}_{y,x}{\theta }_{,xx}+{\zeta }_{y,2}^{\theta ,1}{u}_{y,xx}{\theta }_{,x}+{\zeta }_{z,1}^{\theta }{u}_{z,x}\theta +{\zeta }_{z,1}^{\theta ,1}{u}_{z,x}{\theta }_{,x}+{\zeta }_{z,1}^{\theta ,2}{u}_{z,x}{\theta }_{,xx}+{\zeta }_{z,2}^{\theta ,1}{u}_{z,xx}{\theta }_{,x}+{\zeta }_{\theta }^{\theta ,1}\theta {\theta }_{,x}+{\zeta }_{\theta ,1}^{\theta ,2}{\theta }_{,x}{\theta }_{,xx}.\end{array}

The coefficients α p , i q , j {\alpha }_{p,i}^{q,\hspace{.25em}j} , ζ p , i q , j {\zeta }_{p,i}^{q,\hspace{.25em}j} occurring in (8) and (9) are defined by the formulas: (10) α x , 1 x , 1 = 1 α y , 1 y , 1 = 1 + y , x 2 α z , 1 z , 1 = 1 + z , x 2 α y , 2 y , 2 = y 2 α z , 2 z , 2 = z 2 α θ θ = ( y , x ) 2 + ( z , x ) 2 α θ , 1 θ , 1 = ( ω , x ) 2 + y 2 + z 2 α θ , 2 θ , 2 = ω 2 α x , 1 y , 1 = − y , x α x , 1 y , 2 = − y α x , 1 z , 1 = − z , x α x , 1 z , 2 = − z α x , 1 θ , 1 = − ω , x α x , 1 θ , 2 = − ω α y , 1 z , 1 = y , x z , x α y , 1 z , 2 = y , x z α z , 1 y , 2 = y z , x α y , 2 z , 2 = y z α y , 1 y , 2 = y y , x α z , 1 z , 2 = z z , x α y , 1 θ , 2 = y , x ω α z , 1 θ , 2 = z , x ω α y , 1 θ = − z , x α z , 1 θ = y , x α y , 1 θ , 1 = y , x ω , x − z α z , 1 θ , 1 = z , x ω , x + y α y , 2 θ , 1 = y ω , x α z , 2 θ , 1 = z ω , x α y , 2 θ , 2 = y ω α z , 2 θ , 2 = z ω α θ θ , 1 = z , x z + y , x y α θ , 1 θ , 2 = ω , x ω , \begin{array}{c}{\alpha }_{x,1}^{x,1}=1\hspace{1em}{\alpha }_{y,1}^{y,1}=1+{y}_{,x}^{2}\hspace{1em}{\alpha }_{z,1}^{z,1}=1+{z}_{,x}^{2}\hspace{1em}{\alpha }_{y,2}^{y,2}={y}^{2}{\alpha }_{z,2}^{z,2}={z}^{2}\\ {\alpha }_{\theta }^{\theta }={({y}_{,x})}^{2}+{({z}_{,x})}^{2}\hspace{1em}{\alpha }_{\theta ,1}^{\theta ,1}={({\omega }_{,x})}^{2}+{y}^{2}+{z}^{2}\hspace{1em}{\alpha }_{\theta ,2}^{\theta ,2}={\omega }^{2}\\ {\alpha }_{x,1}^{y,1}=-{y}_{,x}\hspace{1em}{\alpha }_{x,1}^{y,2}=-y\hspace{1em}{\alpha }_{x,1}^{z,1}=-{z}_{,x}\hspace{1em}{\alpha }_{x,1}^{z,2}=-z{\alpha }_{x,1}^{\theta ,1}=-{\omega }_{,x}\hspace{1em}{\alpha }_{x,1}^{\theta ,2}=-\omega \\ {\alpha }_{y,1}^{z,1}={y}_{,x}{z}_{,x}\hspace{1em}{\alpha }_{y,1}^{z,2}={y}_{,x}z\hspace{1em}{\alpha }_{z,1}^{y,2}=y{z}_{,x}\hspace{1em}{\alpha }_{y,2}^{z,2}=yz{\alpha }_{y,1}^{y,2}=y{y}_{,x}\hspace{1em}{\alpha }_{z,1}^{z,2}=z{z}_{,x}\\ {\alpha }_{y,1}^{\theta ,2}={y}_{,x}\omega \hspace{1em}{\alpha }_{z,1}^{\theta ,2}={z}_{,x}\omega \hspace{1em}{\alpha }_{y,1}^{\theta }=-{z}_{,x}\hspace{1em}{\alpha }_{z,1}^{\theta }={y}_{,x}{\alpha }_{y,1}^{\theta ,1}={y}_{,x}{\omega }_{,x}-z\hspace{1em}{\alpha }_{z,1}^{\theta ,1}={z}_{,x}{\omega }_{,x}+y\\ \hspace{.25em}{\alpha }_{y,2}^{\theta ,1}=y{\omega }_{,x}\hspace{1em}{\alpha }_{z,2}^{\theta ,1}=z{\omega }_{,x}\hspace{1em}{\alpha }_{y,2}^{\theta ,2}=y\omega \hspace{1em}{\alpha }_{z,2}^{\theta ,2}=z\omega {\alpha }_{\theta }^{\theta ,1}={z}_{,x}z+{y}_{,x}y\hspace{1em}{\alpha }_{\theta ,1}^{\theta ,2}={\omega }_{,x}\omega ,\end{array} (11) ζ y , 1 y , 1 = y , x y , s ζ z , 1 z , 1 = z , x z , s ζ θ θ = y , x y , s + z , x z , s ζ θ , 1 θ , 1 = ω , x ω , s ζ x , 1 y , 1 = − 1 2 y , s ζ x , 1 z , 1 = − 1 2 z , s ζ x , 1 θ , 1 = − 1 2 ω , s ζ y , 1 θ = − 1 2 z , s ζ z , 1 θ = 1 2 y , s ζ y , 1 y , 2 = 1 2 y y , s ζ z , 2 y , 1 = 1 2 z y , s ζ y , 1 z , 1 = 1 2 ( y , x z , s + z , x y , s ) ζ y , 2 z , 1 = 1 2 y z , s ζ z , 1 z , 2 = 1 2 z z , s ζ y , 1 θ , 1 = 1 2 ( y , x ω , s + ω , x y , s ) ζ z , 1 θ , 1 = 1 2 ( z , x ω , s + ω , x z , s ) ζ y , 1 θ , 2 = 1 2 ω y , s ζ z , 1 θ , 2 = 1 2 ω z , s ζ y , 2 θ , 1 = 1 2 y ω , s ζ z , 2 θ , 1 = 1 2 z ω , s ζ θ θ , 1 = 1 2 ( y y , s + z z , s ) ζ θ , 1 θ , 2 = 1 2 ω ω , s . \begin{array}{c}{\zeta }_{y,1}^{y,1}={y}_{,x}{y}_{,s}\hspace{1em}{\zeta }_{z,1}^{z,1}={z}_{,x}{z}_{,s}\\ {\zeta }_{\theta }^{\theta }={y}_{,x}{y}_{,s}+{z}_{,x}{z}_{,s}\hspace{1em}{\zeta }_{\theta ,1}^{\theta ,1}={\omega }_{,x}{\omega }_{,s}\\ {\zeta }_{x,1}^{y,1}=-\frac{1}{2}{y}_{,s}\hspace{1em}{\zeta }_{x,1}^{z,1}=-\frac{1}{2}{z}_{,s}\hspace{1em}{\zeta }_{x,1}^{\theta ,1}=-\frac{1}{2}{\omega }_{,s}{\zeta }_{y,1}^{\theta }=-\frac{1}{2}{z}_{,s}\hspace{1em}{\zeta }_{z,1}^{\theta }=\frac{1}{2}{y}_{,s}\\ {\zeta }_{y,1}^{y,2}=\frac{1}{2}y{y}_{,s}\hspace{1em}{\zeta }_{z,2}^{y,1}=\frac{1}{2}z{y}_{,s}\hspace{1em}{\zeta }_{y,1}^{z,1}=\frac{1}{2}({y}_{,x}{z}_{,s}+{z}_{,x}{y}_{,s}){\zeta }_{y,2}^{z,1}=\frac{1}{2}y{z}_{,s}\hspace{1em}{\zeta }_{z,1}^{z,2}=\frac{1}{2}z{z}_{,s}\\ {\zeta }_{y,1}^{\theta ,1}=\frac{1}{2}({y}_{,x}{\omega }_{,s}+{\omega }_{,x}{y}_{,s})\hspace{1em}{\zeta }_{z,1}^{\theta ,1}=\frac{1}{2}({z}_{,x}{\omega }_{,s}+{\omega }_{,x}{z}_{,s}){\zeta }_{y,1}^{\theta ,2}=\frac{1}{2}\omega {y}_{,s}\hspace{1em}{\zeta }_{z,1}^{\theta ,2}=\frac{1}{2}\omega {z}_{,s}\\ {\zeta }_{y,2}^{\theta ,1}=\frac{1}{2}y{\omega }_{,s}\hspace{1em}{\zeta }_{z,2}^{\theta ,1}=\frac{1}{2}z{\omega }_{,s}\hspace{1em}{\zeta }_{\theta }^{\theta ,1}=\frac{1}{2}(y{y}_{,s}+z{z}_{,s}){\zeta }_{\theta ,1}^{\theta ,2}=\frac{1}{2}\omega {\omega }_{,s}.\end{array}

The coefficients α p , i q , j {\alpha }_{p,i}^{q,\hspace{.25em}j} , ζ p , i q , j {\zeta }_{p,i}^{q,\hspace{.25em}j}\hspace{.25em} in formulas (8) and (9) are the multipliers of products ∂ x i u p ⋅ ∂ x j u q {\partial }_{x}^{i}{u}_{p}\cdot \hspace{.25em}{\partial }_{x}^{j}{u}_{q} , respectively, where u r = u o x , u y , u z , θ {u}_{r}=\hspace{.25em}{u}_{ox}\hspace{.25em},\hspace{.25em}{u}_{y}\hspace{.25em},\hspace{.25em}{u}_{z}\hspace{.25em},\hspace{.25em}\theta .

Using the constitutive equations [23,27] and taking into account the fact that in the considered case, γ x s = γ s s = 0 , {\gamma }_{xs}={\gamma }_{ss}=0, the stress n x x {n}^{xx} can be written in the form (12) n x x = E ′ g ( 1 + β 2 ) − 2 γ x x , {n}^{xx}=E^{\prime} g{(1+{\beta }^{2})}^{-2}{\gamma }_{xx}, where E ′ = E / 1 − ν 2 E^{\prime} =E/1-{\nu }^{2} , ν \nu is the Poisson’s ratio, E E is the Young’s modulus, and g g is the wall thickness. As it was assumed that γ x s = γ s s = 0 {\gamma }_{xs}={\gamma }_{ss}=0 , the actual forces n x s {n}^{xs} cannot be determined from the constitutive relations, i.e. from the strain–stress relations given above. Forces n x s {n}^{xs} should be determined from the equilibrium equation n x α ∣ α = 0 {{n}^{x\alpha }| }_{\alpha }=0 . Since the metric tensor for the considered coordinate system is known, one can determine the Christoffel symbols and the covariant derivatives. In the special case when γ = 0 \gamma =0 , i.e. the coordinate system is orthogonal (as in the case of I-beams), the equation reduces to the form [27] (13) ( 1 + β 2 ) 3 2 n x s , s + ( 1 + β 2 ) 1 2 ( ( 1 + β 2 ) n x x ) , x = 0 . {\left(,{(1+{\beta }^{2})}^{\frac{3}{2}}{n}^{xs}\right)}_{,s}+{(1+{\beta }^{2})}^{\frac{1}{2}}{((1+{\beta }^{2}){n}^{xx})}_{,x}=0.

As n x x {n}^{xx} is known, the above equality can be easily integrated, where one obtains (14) n x s = − ( 1 + β 2 ) − 3 2 ∫ s − s ( 1 + β 2 ) 1 2 ( ( 1 + β 2 ) n x x ) , x . {n}^{xs}=-{(1+{\beta }^{2})}^{-\frac{3}{2}}\underset{{s}_{-}}{\overset{s}{\int }}{(1+{\beta }^{2})}^{\frac{1}{2}}{((1+{\beta }^{2}){n}^{xx})}_{,x}.

It appears from this equation that when the change in the cross section is so small that one can assume 1 + β 2 ≈ 1 1+{\beta }^{2}\approx 1 , the equation reduces to an equation describing the relationship for a prismatic bar.

A system of displacement equations describing the buckling problem (a second-order theory problem) was derived from the minimal potential energy principle by calculating the first-order variation of this energy: (15) δ ( U S − U 0 − W ) = δ U S − δ U 0 − δ W = 0 , \delta ({U}_{S}-{U}_{0}-W)=\delta {U}_{S}-\delta {U}_{0}-\delta W=0, where U S {U}_{S} is the linear part of the potential energy of the elastic strain, U 0 {U}_{0} is the potential energy generated by the initial stress, and W W is the work of the external load.

The potential energy of the elastic strain is expressed by the formula: (16) U S = 1 2 ∬ E ′ γ x x 2 d A ⁎ d x + 1 2 ∫ − a a GI s ( θ , x ) 2 d x , {U}_{S}=\frac{1}{2}\iint E^{\prime} {\gamma }_{xx}^{2}\text{d}{A}^{\ast }\text{d}x+\frac{1}{2}\underset{-a}{\overset{a}{\int }}{\text{GI}}_{s}{({\theta }_{,x})}^{2}\text{d}x, where GI s \hspace{.25em}{\text{GI}}_{s} is stiffness in Saint-Venant torsion and d A ⁎ = g ⁎ d s , g ⁎ = ( 1 + β 2 ) − 3 2 g \text{d}{A}^{\ast }={g}^{\ast }\text{d}s,\hspace{.5em}{g}^{\ast }={(1+{\beta }^{2})}^{-\frac{3}{2}}g .

The formula for the potential energy generated by the initial stress is as follows: (17) U 0 = ∬ ( n 0 x x γ x x n l + n 0 x s γ x s n l + n 0 s x γ s x n l ) d S = ∬ ( n 0 x x γ x x n l + 2 n 0 x s γ x s n l ) d S , {U}_{0}=\iint ({n}_{0}^{xx}{\gamma }_{xx}^{nl}+{n}_{0}^{xs}{\gamma }_{xs}^{nl}+{n}_{0}^{sx}{\gamma }_{sx}^{nl})\text{d}S=\iint ({n}_{0}^{xx}{\gamma }_{xx}^{nl}+2{n}_{0}^{xs}{\gamma }_{xs}^{nl})\text{d}S, where n 0 x x = E ' g ( 1 + β 2 ) − 2 γ x x 0 , {n}_{0}^{xx}=E\text{'}g{(1+{\beta }^{2})}^{-2}{\gamma }_{xx}^{0}\hspace{.25em}, n 0 x s {n}_{0}^{xs} is defined by formula (14), γ x x 0 {\gamma }_{xx}^{0} is the strain generated by the external loads, determined for the linear problem, γ x x n l {\gamma }_{xx}^{nl} is the nonlinear part of the strains, defined by formula (8), and γ x s n l {\gamma }_{xs}^{nl} is the nonlinear part of the strains, defined by formula (9).

The work of the external load is expressed as follows: (18) W = ∫ a b ( p ∘ u p ) d x = ∫ a b ( p x u x p + p y u y p + p z u z p ) d x , W=\underset{a}{\overset{b}{\int }}({\bf{p}}\circ {{\bf{u}}}^{p})\text{d}x=\underset{a}{\overset{b}{\int }}({p}_{x}{u}_{x}^{p}+{p}_{y}{u}_{y}^{p}+{p}_{z}{u}_{z}^{p})\text{d}x, where p = [ p x , p y , p z ] {\bf{p}}={[}{p}_{x},{p}_{y},{p}_{z}] is the external load vector, u p = [ u x p , u y p , u z p ] T {{\bf{u}}}^{p}={{[}{u}_{x}^{p},{u}_{y}^{p},{u}_{z}^{p}]}^{T} and is the vector of the displacements of the external load application points.

The components of the displacement vector u p {{\bf{u}}}^{p} are defined by the following formulas: (19) u x p = u o x − y p u y , x − z p u z , x − ω p θ , x , u y p = u y − z p θ − 1 2 y p θ 2 , u z p = u z + y p θ − 1 2 z p θ 2 , \begin{array}{c}{u}_{x}^{p}={u}_{ox}-{y}^{p}{u}_{y,x}-{z}^{p}{u}_{z,x}-{\omega }^{p}{\theta }_{,x},\\ {u}_{y}^{p}={u}_{y}-{z}^{p}\theta -\frac{1}{2}{y}^{p}{\theta }^{2},\\ {u}_{z}^{p}={u}_{z}+{y}^{p}\theta -\frac{1}{2}{z}^{p}{\theta }^{2},\end{array} where [ x p , y p , z p ] {[}{x}^{p},{y}^{p},{z}^{p}] is the vector of the external load application point coordinates, and ω p {\omega }^{p} is the sectorial coordinate of the load application point.

As a result of the transformations performed in formula (15) and taking into account the fact that the load vector in the considered problem has the form p = [ 0 , 0 , p z ] {\bf{p}}={[}0,0,{p}_{z}] , one obtains a system of displacement equations. The system comprises four coupled differential equations, three of which are fourth-order equations. The equations are expressed by the following formulas: (20) E ( − A ⁎ u o x , x ) , x + χ − I y , 1 x , 1 0 u y , x − I z , 1 x , 1 0 u z , x − I z , 2 x , 1 0 u z , x x , x = 0 , E{(-{A}^{\ast }{u}_{ox,x})}_{,x}+\chi {\left(-{I}_{\begin{array}{c}y,1\\ x,1\end{array}}^{0}{u}_{y,x}-{I}_{\begin{array}{c}z,1\\ x,1\end{array}}^{0}{u}_{z,x}-{I}_{\begin{array}{c}z,2\\ x,1\end{array}}^{0}{u}_{z,xx}\right)}_{,x}=0, (21) E ( I z ⁎ u y , x x ) , x x + χ I θ , 1 y , 2 0 θ , x + I θ , 2 y , 2 0 θ , x x , x x − I y , 1 x , 1 0 u o x , x + I θ y , 1 0 θ + I θ , 1 y , 1 0 θ , x , x = 0 , E{({I}_{z}^{\ast }{u}_{y,xx})}_{,xx}+\chi \left({\left({I}_{\begin{array}{c}\theta ,1\\ y,2\end{array}}^{0}{\theta }_{,x}+{I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{\theta }_{,xx}\right)}_{,xx}-{\left({I}_{\begin{array}{c}y,1\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta \\ y,1\end{array}}^{0}\theta +{I}_{\begin{array}{c}\theta ,1\\ y,1\end{array}}^{0}{\theta }_{,x}\right)}_{,x}\right)=0, (22) E ( I y ⁎ u z , x x ) , x x + χ I z , 2 x , 1 0 u o x , x + I θ , 1 z , 2 0 θ , x , x x − I z , 1 x , 1 0 u o x , x + I θ z , 1 0 θ + I θ , 1 z , 1 0 θ , x , x − χ p z = 0 , E{({I}_{y}^{\ast }{u}_{z,xx})}_{,xx}+\chi \left({\left({I}_{\begin{array}{c}z,2\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta ,1\\ z,2\end{array}}^{0}{\theta }_{,x}\right)}_{,xx}-{\left({I}_{\begin{array}{c}z,1\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta \\ z,1\end{array}}^{0}\theta +{I}_{\begin{array}{c}\theta ,1\\ z,1\end{array}}^{0}{\theta }_{,x}\right)}_{,x}\right)-\chi {p}_{z}=0, (23) E ( ( I ω ⁎ θ , x x + I ψ ω ⁎ θ , x ) , x x − ( I ψ ⁎ θ , x + I ψ ω ⁎ θ , x x ) , x − G ( I s θ , x ) , x ) + χ I θ , 2 y , 2 0 u y , x x , x x − I θ , 1 y , 1 0 u y , x + I θ , 1 y , 2 0 u y , x x + I θ , 1 z , 1 0 u z , x + I θ , 1 z , 2 0 u z , x x , x + I θ y , 1 0 u y , x + I θ z , 1 0 u z , x − χ ( p z y p − p z z p θ ) = 0 . \begin{array}{c}E({({I}_{\omega }^{\ast }{\theta }_{,xx}+{I}_{\psi \omega }^{\ast }{\theta }_{,x})}_{,xx}-{({I}_{\psi }^{\ast }{\theta }_{,x}+{I}_{\psi \omega }^{\ast }{\theta }_{,xx})}_{,x}-G{({I}_{s}{\theta }_{,x})}_{,x})+\chi \left({\left({I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{u}_{y,xx}\right)}_{,xx}-{\left({I}_{\begin{array}{c}\theta ,1\\ y,1\end{array}}^{0}{u}_{y,x}+{I}_{\begin{array}{c}\theta ,1\\ y,2\end{array}}^{0}{u}_{y,xx}+{I}_{\begin{array}{c}\theta ,1\\ z,1\end{array}}^{0}{u}_{z,x}+{I}_{\begin{array}{c}\theta ,1\\ z,2\end{array}}^{0}{u}_{z,xx}\right)}_{,x}\right.\\ \hspace{1em}\left.+\left({I}_{\begin{array}{c}\theta \\ y,1\end{array}}^{0}{u}_{y,x}+{I}_{\begin{array}{c}\theta \\ z,1\end{array}}^{0}{u}_{z,x}\right)\right)-\chi \hspace{.25em}({p}_{z}\hspace{.25em}{y}^{p}-{p}_{z}\hspace{.25em}{z}^{p}\theta )=0.\end{array}

The associated equations describing the boundary conditions are as follows: (24) [ N δ u o x ] a b = 0 , {{[}N\delta {u}_{ox}]}_{a}^{b}=0,

where N = E ′ A ⁎ u o x , x + χ I y , 1 x , 1 0 u y , x + I z , 1 x , 1 0 u z , x + I z , 2 x , 1 0 u z , x x , N=E^{\prime} \left({A}^{\ast }{u}_{ox,x}+\chi \left({I}_{\begin{array}{c}y,1\\ x,1\end{array}}^{0}{u}_{y,x}+{I}_{\begin{array}{c}z,1\\ x,1\end{array}}^{0}{u}_{z,x}+{I}_{\begin{array}{c}z,2\\ x,1\end{array}}^{0}{u}_{z,xx}\right)\right), (25) [ Q y δ u y ] a b = 0 , {{[}{Q}_{y}\delta {u}_{y}]}_{a}^{b}=0, where Q y = E ′ ( − I z ⁎ u y , x x ) , x − χ I y , 1 x , 1 0 u o x , x + I θ y , 1 0 θ + I θ , 1 y , 1 0 θ , x − χ I θ , 1 y , 2 0 θ , x + I θ , 2 y , 2 0 θ , x x , x , {Q}_{y}=E^{\prime} \left({(-{I}_{z}^{\ast }{u}_{y,xx})}_{,x}-\chi \left({I}_{\begin{array}{c}y,1\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta \\ y,1\end{array}}^{0}\theta +{I}_{\begin{array}{c}\theta ,1\\ y,1\end{array}}^{0}{\theta }_{,x}\right)-\chi {\left({I}_{\begin{array}{c}\theta ,1\\ y,2\end{array}}^{0}{\theta }_{,x}+{I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{\theta }_{,xx}\right)}_{,x}\right), (26) [ M z δ u y , x ] a b = 0 , {{[}{M}_{z}\delta {u}_{y,x}]}_{a}^{b}=0, where M z = − E ′ − I z ⁎ u y , x x + χ I θ , 1 y , 2 0 θ , x + I θ , 2 y , 2 0 θ , x x , {M}_{z}=-E^{\prime} \left(-{I}_{z}^{\ast }{u}_{y,xx}+\chi \left({I}_{\begin{array}{c}\theta ,1\\ y,2\end{array}}^{0}{\theta }_{,x}+{I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{\theta }_{,xx}\right)\right), (27) [ Q z δ u z ] a b = 0 , {{[}{Q}_{z}\delta {u}_{z}]}_{a}^{b}=0, where Q z = E ′ − ( I y ⁎ u z , x x ) , x + χ I z , 1 x , 1 0 u o x , x + I θ z , 1 0 θ + I θ , 1 z , 1 0 θ , x − χ I z , 2 x , 1 0 u o x , x + I θ , 1 z , 2 0 θ , x , x , {Q}_{z}=E^{\prime} \left(-{({I}_{y}^{\ast }{u}_{z,xx})}_{,x}+\chi \left({I}_{\begin{array}{c}z,1\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta \\ z,1\end{array}}^{0}\theta +{I}_{\begin{array}{c}\theta ,1\\ z,1\end{array}}^{0}{\theta }_{,x}\right)-\chi {\left({I}_{\begin{array}{c}z,2\\ x,1\end{array}}^{0}{u}_{ox,x}+{I}_{\begin{array}{c}\theta ,1\\ z,2\end{array}}^{0}{\theta }_{,x}\right)}_{,x}\right), (28) [ M y δ u z , x ] a b = 0 , {{[}{M}_{y}\delta {u}_{z,x}]}_{a}^{b}=0, where M y = E ′ − I y ⁎ u z , x x + χ − I z , 2 x , 1 0 u o x , x − I z , 2 y , 1 0 u y , x − I z , 2 y , 2 0 u y , x x − I z , 2 z , 1 0 u z , x − I z , 2 z , 2 0 u z , x x − I θ , 1 z , 2 0 θ , x − I θ , 2 z , 2 0 θ , x x , {M}_{y}=E^{\prime} \left(-{I}_{y}^{\ast }{u}_{z,xx}+\chi \left(-{I}_{\begin{array}{c}z,2\\ x,1\end{array}}^{0}{u}_{ox,x}-{I}_{\begin{array}{c}z,2\\ y,1\end{array}}^{0}{u}_{y,x}-{I}_{\begin{array}{c}z,2\\ y,2\end{array}}^{0}{u}_{y,xx}\hspace{2em}-{I}_{\begin{array}{c}z,2\\ z,1\end{array}}^{0}{u}_{z,x}-{I}_{\begin{array}{c}z,2\\ z,2\end{array}}^{0}{u}_{z,xx}-{I}_{\begin{array}{c}\theta ,1\\ z,2\end{array}}^{0}{\theta }_{,x}-{I}_{\begin{array}{c}\theta ,2\\ z,2\end{array}}^{0}{\theta }_{,xx}\right)\right), (29) [ M x δ θ ] a b = 0 , {{[}{M}_{x}\delta \theta ]}_{a}^{b}=0, where M x = E ′ ( ( − I ω ⁎ θ , x x − I ψ ω ⁎ θ , x ) , x + ( I ψ ⁎ θ , x + I ψ ω ⁎ θ , x x ) − χ I θ , 1 y , 1 0 u y , x + I θ , 1 y , 2 0 u y , x x + I θ , 1 z , 1 0 u z , x + I θ , 1 z , 2 0 u z , x x − χ I θ , 2 y , 2 0 u y , x x , x + GI s θ , x , \begin{array}{c}{M}_{x}=E^{\prime} ({(-{I}_{\omega }^{\ast }{\theta }_{,xx}-{I}_{\psi \omega }^{\ast }{\theta }_{,x})}_{,x}+({I}_{\psi }^{\ast }{\theta }_{,x}+{I}_{\psi \omega }^{\ast }{\theta }_{,xx})\\ \hspace{1em}\left.-\chi \left({I}_{\begin{array}{c}\theta ,1\\ y,1\end{array}}^{0}{u}_{y,x}+{I}_{\begin{array}{c}\theta ,1\\ y,2\end{array}}^{0}{u}_{y,xx}+{I}_{\begin{array}{c}\theta ,1\\ z,1\end{array}}^{0}{u}_{z,x}+{I}_{\begin{array}{c}\theta ,1\\ z,2\end{array}}^{0}{u}_{z,xx}\right)-\chi {\left({I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{u}_{y,xx}\right)}_{,x}\right)+{\text{GI}}_{s}{\theta }_{,x},\end{array} (30) [ B δ θ , x ] a b = 0 , {{[}B\delta {\theta }_{,x}]}_{a}^{b}=0, where B = E ′ I ω ⁎ θ , x x + I ψ ω ⁎ θ , x + χ I θ , 2 y , 2 0 u y , x x . B=E^{\prime} \left({I}_{\omega }^{\ast }{\theta }_{,xx}+{I}_{\psi \omega }^{\ast }{\theta }_{,x}+\chi \left({I}_{\begin{array}{c}\theta ,2\\ y,2\end{array}}^{0}{u}_{y,xx}\right)\right).

The coefficients in formulas (20)–(30) are defined by the following formulas: (31) A ⁎ = ∫ S d A ⁎ , I y ⁎ = ∫ S z 2 d A ⁎ , I z ⁎ = ∫ S y 2 d A ⁎ I ω ⁎ = ∫ S ω 2 d A ⁎ , I ψ ω ⁎ = ∫ S ω ψ d A ⁎ , I ψ ⁎ = ∫ S ψ 2 d A ⁎ , \begin{array}{c}{A}^{\ast }=\mathop{\int }\limits_{S}\text{d}{A}^{\ast },\hspace{1em}{I}_{y}^{\ast }=\mathop{\int }\limits_{S}{z}^{2}\text{d}{A}^{\ast },\hspace{1em}{I}_{z}^{\ast }=\mathop{\int }\limits_{S}{y}^{2}\text{d}{A}^{\ast }\\ {I}_{\omega }^{\ast }=\mathop{\int }\limits_{S}{\omega }^{2}\text{d}{A}^{\ast },\hspace{1em}{I}_{\psi \omega }^{\ast }=\mathop{\int }\limits_{S}\omega \psi \text{d}{A}^{\ast },\hspace{1em}{I}_{\psi }^{\ast }=\mathop{\int }\limits_{S}{\psi }^{2}\text{d}{A}^{\ast },\end{array} (32) I q , 2 p , 2 0 = ∫ S n 0 x x α q , 2 p , 2 d s ⁎ , when p , q = y , z , θ ; I q , 2 x , 1 0 = ∫ S n 0 x x α q , 2 x , 1 d s ⁎ , when q = y , z , θ ; I x , 1 x , 1 0 = ∫ S n 0 x x α x , 1 x , 1 d s ⁎ I q , j p , i 0 = ∫ S n 0 x x α q , j p , i + 2 n 0 x s ζ q , j x , i d s ⁎ − in other cases, d s ⁎ = 1 + β 2 d s . \begin{array}{c}{I}_{\begin{array}{c}q,2\\ p,2\end{array}}^{0}=\mathop{\int }\limits_{S}{n}_{0}^{xx}{\alpha }_{\begin{array}{c}q,2\\ p,2\end{array}}\text{d}{s}^{\ast },\hspace{1em}\text{when}\hspace{.5em}p,q=y,\hspace{.5em}z,\hspace{.5em}\theta ;\\ {I}_{\begin{array}{c}q,2\\ x,1\end{array}}^{0}=\mathop{\int }\limits_{S}{n}_{0}^{xx}{\alpha }_{\begin{array}{c}q,2\\ x,1\end{array}}\text{d}{s}^{\ast },\hspace{1em}\text{when}\hspace{.5em}q=y,\hspace{.5em}z,\hspace{.5em}\theta ;\\ {I}_{\begin{array}{c}x,1\\ x,1\end{array}}^{0}=\mathop{\int }\limits_{S}{n}_{0}^{xx}{\alpha }_{\begin{array}{c}x,1\\ x,1\end{array}}\text{d}{s}^{\ast }\\ {I}_{\begin{array}{c}q,\hspace{.25em}j\\ p,\hspace{.25em}i\end{array}}^{0}=\mathop{\int }\limits_{S}\left({n}_{0}^{xx}{\alpha }_{\begin{array}{c}q,j\\ p,i\end{array}}+2{n}_{0}^{xs}{\zeta }_{\begin{array}{c}q,j\\ x,i\end{array}}\right)\text{d}{s}^{\ast }-\text{in}\hspace{.5em}\text{other}\hspace{.5em}\text{cases,}\text{d}{s}^{\ast }=\sqrt{1+{\beta }^{2}}\text{d}s.\end{array}

The theorem given by Lewanowicz [28] was used to solve the considered problem. The theorem describes a method of solving ordinary differential equations using Gegenbauer polynomials (a detailed description of this method can be found in Lewanowicz [28] and in the present authors’ publications [23,24,25]). Here, zero-order Gegenbauer polynomials ( λ = 0 ) (\lambda =0) , i.e. Chebyshev polynomials of the first order, were used to solve the problem, (33) T k ( x ) = ∑ m = 0 [ k / 2 ] ( − 1 ) m Γ ( k − m ) m ! ( k − 2 m ) ! ( 2 x ) k − 2 m . {T}_{k}(x)=\mathop{\sum }\limits_{m=0}^{{[}k/2]}{(-1)}^{m}\frac{\text{Γ}(k-m)}{m\!(k-2m)\!}{(2x)}^{k-2m}.

The system of differential equations being solved has the form (34) ∑ i = 0 n P i ( x ) f ( i ) ( x ) = P ( x ) . \mathop{\sum }\limits_{i=0}^{n}{{\bf{P}}}_{i}(x)\hspace{.25em}{{\bf{f}}}^{(i)}(x)={\bf{P}}(x).

In the method presented by Paszkowski [29] and Lewanowicz [28], solutions are obtained in the form of the series (35) f ( x ) = 1 2 b 0 [ f ] T 0 ( x ) + ∑ k = 1 ∞ b k [ f ] T k ( x ) , {\bf{f}}(x)=\frac{1}{2}{b}_{0}{[}{\bf{f}}]{T}_{0}(x)+\mathop{\sum }\limits_{k=1}^{\infty }{b}_{k}{[}{\bf{f}}]{T}_{k}(x), where f = f 1 f 2 ⋯ f m T , b k [ f ] = b k [ f 1 ] b k [ f 2 ] ⋯ b k [ f m ] T {\bf{f}}={\left[\begin{array}{cc}{f}_{1}& {f}_{2}\cdots {f}_{m}\end{array}\right]}^{{\bf{T}}},\hspace{1em}{b}_{k}{[}{\bf{f}}]={\left[\begin{array}{cc}{b}_{k}{[}{f}_{1}]& {b}_{k}{[}{f}_{2}]\cdots {b}_{k}{[}{f}_{m}]\end{array}\right]}^{{\bf{T}}} , and b k [ f i ] {b}_{k}{[}{f}_{i}] is the kth term of the Chebyshev series expansion of the function f i {f}_{i} . It was shown in Lewanowicz [28] that the initial system of equations (34) is equivalent to the following system: (36) ∑ i = 0 n ( Q i ( x ) f ( x ) ) ( i ) = P ( x ) , \mathop{\sum }\limits_{i=0}^{n}{({{\bf{Q}}}_{i}(x){\bf{f}}(x))}^{(i)}={\bf{P}}(x), where the matrix Q i ( x ) {{\bf{Q}}}_{i}(x) is defined by the formula (37) Q i ( x ) = ∑ j = i n ( − 1 ) j − i j j − i ( P j ( x ) ) ( j − i ) , i = 0 , 1 , … , n . {{\bf{Q}}}_{i}(x)=\mathop{\sum }\limits_{j=i}^{n}{(-1)}^{j-i}\left(\begin{array}{c}j\\ j-i\end{array}\right){({{\bf{P}}}_{j}(x))}^{(j-i)},\hspace{1em}i=0,\hspace{.25em}1,\hspace{.25em}\ldots ,\hspace{.25em}n.

The system of equations (36) can be reduced (see Lewanowicz [28]) to the following infinite system of equations: (38) ∑ i = 0 n 2 i ∑ m = 0 n − i ϱ n i m ( k ) b k − n + i + 2 m [ Q i ( x ) f ( x ) ] = ∑ m = 0 n ϱ n 0 m ( k ) b k − n + 2 m [ P ( x ) ] for k ≥ n , \mathop{\sum }\limits_{i=0}^{n}{2}^{i}\mathop{\sum }\limits_{m=0}^{n-i}{{\varrho }}_{nim}(k)\hspace{.25em}{b}_{k-n+i+2m}{[}{{\bf{Q}}}_{i}(x)\hspace{.25em}{\bf{f}}(x)]=\mathop{\sum }\limits_{m=0}^{n}{{\varrho }}_{n0m}(k){b}_{k-n+2m}{[}{\bf{P}}(x)]\hspace{1em}\text{for}\hspace{.5em}k\ge n, where (39) ϱ i j m ( k ) = ( − 1 ) m i − j m ( k − i ) j + m ( k − i + j + 2 m ) ( k + m + 1 ) i − m ( k 2 − i 2 ) − 1 , ( α ) 0 = 1 , ( α ) k = α ( α + 1 ) … ( α + k − 1 ) for k ≥ 1 . \begin{array}{c}{{\varrho }}_{ijm}(k)\hspace{.25em}={(-1)}^{m}\left(\begin{array}{c}i-j\\ m\end{array}\right){(k-i)}_{j+m}(k-i+j+2m)\hspace{.25em}{(k+m+1)}_{i-m}\hspace{.25em}{({k}^{2}-{i}^{2})}^{-1},\\ {(\alpha )}_{0}=1,\hspace{1em}{(\alpha )}_{k}=\alpha (\alpha +1)\ldots (\alpha +k-1)\hspace{.5em}\text{for}\hspace{.5em}k\ge 1.\end{array}