**2. Work Methods**

#### *2.1. Analytical Model for Buckling of a Column with Pin Connections at Ends with Stepwise Variable Section*

In Figure 1, a geometrical model is shown of a bar having a stepwise variable circular cross-section whose bottom end is pin-connected, while the upper end is simply supported. The undeformed and deformed shapes of the bar under compression are shown in Figure 1a. The bar consists in three portions. The first and the third portions have the same values for length *l* and for the second moment of inertia *I* of the cross-section. The second bar portion has a cross-section whose moment of inertia *I*<sup>1</sup> is equal to *k*<sup>1</sup> *I*, with *k*<sup>1</sup> > 1, and its length is equal to 2*k*2*l*. In other words, the parameters *k*<sup>1</sup> and *k*<sup>2</sup> represent the ratios of the lengths and of the second moment of inertia, respectively, corresponding to the second and first portions of the bar (Figure 1a). The total length of the column is denoted with *lc* (Figure 1a). All portions of the bar are made of the same isotropic material having a modulus of elasticity *E*. The column is symmetric with respect to its midpoint (Figure 1a), and consequently, the analysis model may be reduced to half of the bar, considering the symmetry conditions (Figure 1b). The shape of the bending-moment diagram for half of the bar is shown in Figure 1c.

**Figure 1.** Bar analyzed having a stepwise variable cross-section whose bottom end is pin-connected, while the upper end is simply supported: (**a**) geometrical model of the entire bar and deformed shape at stability loss; (**b**) geometrical model considering the symmetry condition; (**c**) bendingmoment diagram.

The main purpose of this subsection is to compare the critical buckling force corresponding to a column having a stepwise variable circular cross-section with the critical buckling force corresponding to a column having the same length and a constant crosssection whose moment of inertia is equal to *I*. For this purpose, the main objectives are to compute the critical buckling force corresponding to the column having a stepwise variable cross-section and to compare it with the one corresponding to the column with a constant cross-section, as well as an analysis of the rational shapes for buckling in the case of columns having stepwise variable circular cross-sections.

Due to symmetry, the buckling analysis was made by considering a model of half of a column (Figure 1b) with boundary conditions corresponding to the midpoint of the column located on the horizontal symmetry axis.

It was assumed that the column buckled in the elastic domain. This meant that Bernoulli's hypothesis and Euler's relation were valid at buckling.

The bending moment *Mb*1(*x*) developed at buckling and caused by the compressive force *P* at the level of the arbitrary cross-section located on the first portion of the column was computed using Equation (1):

$$M\_{b1}(\mathbf{x}) = -Pv\_1(\mathbf{x}),\tag{1}$$

where *v*1(*x*) represents the deflection function of the arbitrary cross-section of the first column portion located at distance *x* with respect to the end of the column (Figure 1b).

In the same manner, the bending moment *Mb*2(*x*) developed at buckling at the level of the arbitrary cross-section located on the second portion of the column was computed using Equation (2):

$$M\_{b2}(\mathbf{x}) = -P\upsilon\_2(\mathbf{x}),\tag{2}$$

where *v*2(*x*) represents the deflection function of the arbitrary cross-section of the second column portion located at distance *x* with respect to the end of the column (Figure 1b).

To compute the critical buckling force, we used an analytical method. The differential equation of the approximate deformed median fiber corresponding to the first portion of the bar is given in Equation (3):

$$\frac{d^2\upsilon\_1(\mathbf{x})}{d\mathbf{x}^2} = -\frac{M\_{b1}(\mathbf{x})}{EI}.\tag{3}$$

Replacing Equation (1), Equation (3) became the following [2]:

$$\frac{d^2\upsilon\_1(\mathbf{x})}{dx^2} + \frac{P\upsilon\_1(\mathbf{x})}{EI} = 0.\tag{4}$$

In the same manner, the differential equation of the approximate deformed median fiber corresponding to the second portion of the bar is given in Equation (5) [2]:

$$\frac{d^2v\_2(\mathbf{x})}{d\mathbf{x}^2} + \frac{Pv\_2(\mathbf{x})}{k\_1EI} = 0.\tag{5}$$

The notation *α* was introduced for the ratio given in Equation (6) [2]:

$$\alpha = \sqrt{\frac{P}{k\_1 EI}} \quad \text{or} \quad \alpha^2 = \frac{P}{k\_1 EI}. \tag{6}$$

By using Equation (6), differential Equations (4) and (5) of the approximate deformed median fibers for the bar portions became:

1. Equation (7) for the first portion of the column [2]:

$$\frac{d^2v\_1(\mathbf{x})}{d\mathbf{x}^2} + \left(\mathfrak{a}\sqrt{k\_1}\right)^2 \cdot v\_1(\mathbf{x}) = 0,\tag{7}$$

2. Equation (8) for the second portion of the column [2]:

$$\frac{d^2v\_2(\mathbf{x})}{d\mathbf{x}^2} + \mathfrak{a}^2 \cdot v\_2(\mathbf{x}) = 0. \tag{8}$$

Solutions of inhomogeneous second-order differential Equations (7) and (8) were given by Equations (9) and (10) for the first portion and for the second portion of the column, respectively:

$$w\_1(\mathbf{x}) = \mathbb{C}\_1 \sin \left( \mathfrak{a} \sqrt{k\_1} \mathbf{x} \right) + \mathbb{C}\_2 \cos \left( \mathfrak{a} \sqrt{k\_1} \mathbf{x} \right), \tag{9}$$

$$w\_2(\mathbf{x}) = \mathbb{C}\_3 \sin(a\mathbf{x}) + \mathbb{C}\_4 \cos(a\mathbf{x}).\tag{10}$$

Integration constants *Ci i* = 1, 4 were computed using the boundary conditions given in Equation (11):

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{0}, & \upsilon\_1(\mathbf{0}) = \mathbf{0};\\ \text{for } \mathbf{x} = l + k\_2 l, & \upsilon\_2(l + k\_2 l) = f;\\ \text{for } \mathbf{x} = l + k\_2 l, & \frac{\text{d}\upsilon\_2}{\text{d}\mathbf{x}}(l + k\_2 l) = \mathbf{0}, \end{cases} \tag{11}$$

where *f* represents the deflection of the midpoint of the bar at buckling. The continuity conditions for the deformed shape of the median fiber of the column are given in Equation (12) at the level of the bar cross-section located at distance *l* with respect to the upper simply supported end of the column:

$$\begin{cases} \text{for } \mathbf{x} = l, \quad v\_1(l) = v\_2(l);\\ \text{for } \mathbf{x} = l, \quad \frac{\mathrm{d}v\_1}{\mathrm{d}\mathbf{x}}(l) = \frac{\mathrm{d}v\_2}{\mathrm{d}\mathbf{x}}(l). \end{cases} \tag{12}$$

The first derivatives of the functions *v*1(*x*) and *v*2(*x*) of the deflections of the arbitrary cross-section corresponding to each column portion were computed using Equations (9) and (10), respectively. In fact, these derivatives represented the functions of the rotations of the arbitrary cross-section and were expressed by Equations (13) and (14):

$$\frac{d\mathbf{v}\_1(\mathbf{x})}{d\mathbf{x}} = a\sqrt{k\_1}\mathbb{C}\_1\left[\cos\left(a\sqrt{k\_1}\mathbf{x}\right) - \mathbb{C}\_2\sin\left(a\sqrt{k\_1}\mathbf{x}\right)\right],\tag{13}$$

$$\frac{\mathrm{d}v\_{2}(x)}{\mathrm{d}x} = a[\mathrm{C\_{3}}\cos(kx) - \mathrm{C\_{4}}\sin(ax)].\tag{14}$$

By replacing Equations (9), (10), (13), and (14) in the boundary conditions given by Equation (11) and in the continuity conditions given in Equation (12), the system of Equation (15) was obtained:

$$\begin{cases} \mathsf{C}\_{2} = 0; \\ \begin{aligned} \mathsf{C}\_{3}\sin[al(1+k\_{2})] + \mathsf{C}\_{4}\cos[al(1+k\_{2})] &= f; \\ \varkappa\{\mathsf{C}\_{3}\cos[al(1+k\_{2})] - \mathsf{C}\_{4}\sin[al(1+k\_{2})] \} &= 0; \\ \mathsf{C}\_{1}\sin[al\sqrt{k\_{1}}] + \mathsf{C}\_{2}\cos[al\sqrt{k\_{1}}] &= \mathsf{C}\_{3}\sin(al) + \mathsf{C}\_{4}\cos(al); \\ \varkappa\sqrt{k\_{1}}[\mathsf{C}\_{1}\cos(al\sqrt{k\_{1}}) - \mathsf{C}\_{2}\sin(al\sqrt{k\_{1}})] &= \varkappa[\mathsf{C}\_{3}\cos(al) - \mathsf{C}\_{4}\sin(al)], \end{aligned} \end{cases} \tag{15}$$

whose unknown quantities are the integration constants *Ci i* = 1, 4 and the deflection *f* of the midpoint of the column.

Equation system (15) was reduced practically to the following system of four equations:

$$\begin{cases} \mathcal{C}\_3 \sin[al(1+k\_2)] + \mathcal{C}\_4 \cos[al(1+k\_2)] - f = 0; \\ \mathcal{C}\_3 \cos[al(1+k\_2)] - \mathcal{C}\_4 \sin[al(1+k\_2)] = 0; \\ \mathcal{C}\_1 \sin(al\sqrt{k\_1}) - \mathcal{C}\_3 \sin(al) - \mathcal{C}\_4 \cos(al) = 0; \\ \mathcal{C}\_1 \sqrt{k\_1} \cos(al\sqrt{k\_1}) - \mathcal{C}\_3 \cos(al) + \mathcal{C}\_4 \sin(al) = 0. \end{cases} \tag{16}$$

The homogeneous system of Equation (16) had nonzero solutions for *Ci i* = 1, 4 and *f* if the determinant of the coefficients was zero. This condition led to Equation (17):

$$\text{Det}\begin{bmatrix} 0 & \sin[al(1+k\_2)] & \cos[al(1+k\_2)] & -1\\ 0 & \cos[al(1+k\_2)] & -\sin[al(1+k\_2)] & 0\\ \sin(al\sqrt{k\_1}) & -\sin(al) & -\cos(al) & 0\\ \sqrt{k\_1}\cos(al\sqrt{k\_1}) & -\cos(al) & \sin(al) & 0 \end{bmatrix} = 0. \tag{17}$$

The notation *ξ* was introduced and computed with Equation (18):

$$
\xi = \mathfrak{a}l\_\prime \tag{18}
$$

and Equation (17) became:

$$\text{Det}\begin{bmatrix} 0 & \sin[\xi(1+k\_2)] & \cos[\xi(1+k\_2)] & -1\\ 0 & \cos[\xi(1+k\_2)] & -\sin[\xi(1+k\_2)] & 0\\ \sin(\xi\sqrt{k\_1}) & -\sin\xi & -\cos\xi & 0\\ \sqrt{k\_1}\cos(\xi\sqrt{k\_1}) & -\cos\xi & \sin\xi & 0 \end{bmatrix} = 0,\tag{19}$$

whose unknown is *ξ*. To compute the critical buckling force, the minimum value of the absolute values of the solutions *ξ* must be used because, for other solutions *ξ*, the value of the critical force would be greater.

*ξmin* denoted the minimum value of the absolute values of the solutions *ξ* of Equation (19).

From Equation (18), the minimum value of *α* could be computed using Equation (20):

$$
\alpha\_{\rm min} = \mathfrak{J}\_{\rm min} / l. \tag{20}
$$

Equation (20) was replaced in relation (6) in order to compute the critical buckling force for a column having a stepwise variable cross-section:

$$P\_{cr} = k\_1 a\_{\text{min}}^2 EI = k\_1 \zeta\_{\text{min}}^2 \frac{EI}{l^2}.\tag{21}$$

The length *l* of the first portion was expressed in the function of the total length *lc* of the column using Equation (22) according to Figure 1a:

$$l = l\_c / [2(1 + k\_2)].\tag{22}$$

Equation (22) was replaced in Equation (21) in order to compute the critical buckling force *Pcr* for a column having a stepwise variable cross-section:

$$P\_{cr} = 4k\_1(1+k\_2)^2 \xi\_{\text{min}}^2 \frac{EI}{l\_c^2} = \left(\frac{2\sqrt{k\_1}(1+k\_2)\xi\_{\text{min}}}{\pi}\right)^2 \frac{\pi^2 EI}{l\_c^2}.\tag{23}$$

On the other hand, the critical buckling force *Pcr*<sup>0</sup> for a column with a constant crosssection whose moment of inertia of the section was *I* having the same total length *lc* with pin connections at its ends was computed as follows:

$$P\_{\rm cr0} = \pi^2 EI/l\_\odot^2. \tag{24}$$

Using Equations (23) and (24), the normalized buckling force (denoted with *c*) was computed as being equal with the ratio between the critical buckling force *Pcr* for a column having a stepwise variable cross-section and the critical buckling force *Pcr*<sup>0</sup> for a column having a constant cross-section whose moment of inertia of the section was *I*:

$$c = \frac{P\_{cr}}{P\_{cr0}} = \left(\frac{2\sqrt{k\_1}(1+k\_2)\xi\_{\text{min}}}{\pi}\right)^2. \tag{25}$$

It is said that a column is rationally designed if the critical buckling force is increased while the mass of the material of the column is optimal. In structure design, both the mass of the material, which influences the material costs, and the weight of the structure are also important. The manufacturing costs increase due to the additional manufacturing operations required for a column having a stepwise variable cross-section.

In this research, it was said that a column having a stepwise variable cross-section was rationally designed if its ratio between the critical buckling force *Pcr* and the total volume *V* of the column was greater than the similar ratio computed for a column with a constant cross-section whose moment of inertia of the section was *I* with the same total length *lc*.

The ratio between the critical buckling force *Pcr* and the total volume *V* of the column was computed as follows, taking into account the geometry of the column shown in Figure 1a:

$$\frac{P\_{cr}}{V} = \frac{P\_{cr}}{2l(A + A\_1k\_2)},\tag{26}$$

where *A* represents the cross-sectional areas for the first and the third portions of the column having the length *l*, while *A*<sup>1</sup> represents the area of the cross-section for the second column portion whose length is equal to 2*k*2*l*.

It was assumed that each portion of the column had a circular cross-section. It may be remarked in Equation (27), which gives the ratio between the moment of inertia *I* and the second power of the area of the cross-section *A* for the first portion, whose diameter is denoted with *d*:

$$\frac{1}{A^2} = \frac{\pi d^4 / 64}{\left(\pi d^2 / 4\right)^2} = \frac{1}{4\pi} \tag{27}$$

or

$$I = A^2 / 4\pi.$$

In the same manner, Equation (29) was written for the second portion of the column:

$$I\_1 = A\_1^2 / 4\pi. \tag{29}$$

On the other hand, the relation between the moments of inertia *I*<sup>1</sup> and *I* corresponding to the first two portions of the column, respectively, is given in Equation (30):

$$I\_1 = k\_1 I. \tag{30}$$

Equations (28) and (29) were replaced in Equation (30), and it obtained Equation (31):

$$A\_1 = A\sqrt{k\_1}.\tag{31}$$

Then, relations (22) and (31) were replaced in relation (26), which became:

$$\frac{P\_{cr}}{V} = \frac{(1+k\_2)P\_{cr}}{l\_c A \left(1+k\_2\sqrt{k\_1}\right)}.\tag{32}$$

From Equation (25), the critical buckling force *Pcr* for a column having a stepwise variable cross-section could be computed in the function of the critical buckling force *Pcr*<sup>0</sup> of a column with a constant cross-section:

$$P\_{cr} = \mathcal{c}P\_{cr0}.\tag{33}$$

The volume of the column with a constant cross-section was:

$$V\_0 = l\_\text{c} A.\tag{34}$$

Replacing Equations (33) and (34) in Equation (32) obtained Equation (35):

$$\frac{P\_{cr}}{V} = \frac{(1+k\_2)c}{1+k\_2\sqrt{k\_1}} \frac{P\_{cr0}}{V\_0} \,. \tag{35}$$

which led to Equation (36), which computed the rationality factor *krat*:

$$k\_{\rm rat} = \frac{P\_{cr}/V}{P\_{cr0}/V\_0} = \frac{(1+k\_2)c}{1+k\_2\sqrt{k\_1}}\,. \tag{36}$$

Replacing the ratio *c* given by Equation (25) in Equation (36) led to Equation (37):

$$k\_{\rm rdl} = \frac{P\_{\rm cr}/V}{P\_{\rm cr0}/V\_0} = \frac{(1+k\_2)}{1+k\_2\sqrt{k\_1}} \left(\frac{2\sqrt{k\_1}(1+k\_2)\xi\_{\rm min}}{\pi}\right)^2 = \frac{4\xi\_{\rm min}^2 k\_1 (1+k\_2)^3}{\pi^2 \left(1+k\_2\sqrt{k\_1}\right)},\tag{37}$$

which was used to compute the ratio between the rationality factor *Pcr*/*V* corresponding to a column having a stepwise variable cross-section and the rationality factor *Pcr*0/*V*<sup>0</sup> corresponding to a column with a constant cross-section.

*2.2. Numerical Modeling and Simulation for Loss in Stability of a Column with Pin Connections at Ends with Stepwise Variable Cross-Section*

If the compressive stress, which acts on the slenderness of a column, is variable, then the relation between force and displacement may be written with Equation (38):

$$\{dP\} = [K\_T] \{\Delta U\},\tag{38}$$

where infinitesimal variation in the force is denoted with d*P*. A similar relation may be written with Equation (39):

$$\{\Delta P\} = [\mathbb{K}\_{\mathbb{S}}]\{\Delta lI\},\tag{39}$$

where finite variation in the force is denoted with Δ*P*. In Equations (38) and (39), [*KT*] = [*KT*{*U*}] and [*KT*] = [*KT*{*U*}] represent the tangent stiffness matrix and the secant stiffness matrix, respectively. The finite variation in the displacements was computed with Equation (40):

$$\{\Delta \mathcal{U}\} = \left[K\tau(\mathcal{U})\right]^{-1} \{\Delta P\} = \frac{\left[K\_T(\mathcal{U})\right]^\*}{\det[K\_T(\mathcal{U})]} \{\Delta P\},\tag{40}$$

where [*KT*(*U*)] ∗ is the adjunct matrix of the stiffness matrix.

The phenomenon of loss in stability (transition from one equilibrium shape to another equilibrium shape) takes place when the displacements tend toward infinity for a variation Δ*P* in the compressive force. From a mathematical point of view, this condition is fulfilled if the determinant of the tangent stiffness matrix [*KT*] is equal to zero, which was expressed by Equation (41):

$$\det[\mathbb{K}\_T(\mathcal{U})] = 0.\tag{41}$$

Equation (41) could be written with Equation (42):

$$\det\left[\left[\mathbb{K}\right]-\lambda\left[\mathbb{K}\_{\mathbb{S}}\right]\right]=\mathbf{0},\tag{42}$$

where [*K*] represents the elastic stiffness matrix of the structure (e.g., the column) obtained by assembling of the stiffness matrices [*Ke*] corresponding to the finite elements that form the structure; *Kg* is the geometric stiffness matrix of the structure obtained by assembling of the geometric stiffness matrices *Kge* corresponding to the finite elements that form the structure; and *λ* is the common multiplier of the axial forces *N* acting in the slender bar.

For the finite element of the double-embedded bar type (Figure 2), the elastic stiffness matrix [*Ke*] and the geometric stiffness matrix *Kge* were expressed with Equations (43) and (44), respectively:

 *Kge* = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *EA <sup>l</sup>* 0 0 0 <sup>12</sup>*EI l*3 6*EI l*2 0 <sup>6</sup>*EI l*2 4*EI l* −*EA <sup>l</sup>* 0 0 <sup>0</sup> <sup>−</sup>12*EI l*3 6*EI l*2 <sup>0</sup> <sup>−</sup>6*EI l*2 2*EI l* −*EA <sup>l</sup>* 0 0 <sup>0</sup> <sup>−</sup>12*EI <sup>l</sup>*<sup>3</sup> <sup>−</sup>6*EI l*2 0 <sup>6</sup>*EI l*2 2*EI l EA <sup>l</sup>* 0 0 0 <sup>12</sup>*EI <sup>l</sup>*<sup>3</sup> <sup>−</sup>6*EI l*2 <sup>0</sup> <sup>−</sup>6*EI l*2 4*EI l* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (43) *Kge* <sup>=</sup> *<sup>N</sup> l* ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 00 0 0 <sup>6</sup> <sup>5</sup> *<sup>l</sup>* 10 0 *<sup>l</sup>* <sup>10</sup> <sup>2</sup>*<sup>l</sup>* 2 15 00 0 <sup>0</sup> <sup>−</sup><sup>6</sup> <sup>5</sup> *<sup>l</sup>* 10 <sup>0</sup> <sup>−</sup> *<sup>l</sup>* <sup>10</sup> <sup>−</sup>*<sup>l</sup>* 2 30 00 0 <sup>0</sup> <sup>−</sup><sup>6</sup> <sup>5</sup> <sup>−</sup> *<sup>l</sup>* 10 0 ll <sup>10</sup> <sup>−</sup>*<sup>l</sup>* 2 30 00 0 0 <sup>6</sup> <sup>5</sup> <sup>−</sup> ll 10 <sup>0</sup> <sup>−</sup> ll <sup>10</sup> <sup>2</sup>*<sup>l</sup>* 2 15 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (44)

where *l* , *A* , and *I* are the length, area of the cross-section, and the second moment of inertia of the cross-section, respectively, corresponding to the finite elements double-embedded at both ends.

**Figure 2.** Characteristics of the finite elements of the bar double-embedded at both ends in terms of both the internal forces and flections (displacement and rotation) developed at the nodes.

In this context, the solving of the stability equation involved solving a problem of eigenvectors and eigenvalues. The solutions of the stability equation were the eigenvalues *λk k* = 1, *n* corresponding to the multiplier of the axial forces. For the eigenvalues *λk k* = 1, *n* , the corresponding eigenvectors {*Uk*}  *k* = 1, *n* were determined, which represented the geometric shapes (equilibrium shapes) of the loss in stability. From a practical point of view, only the lowest eigenvalue *λ*min was of interest, the other values being of interest just from a theoretical point of view.

The numerical model for the calculation of the critical buckling force was validated for a bar having a pin connection at one end and a simple support at the other end. Considering the numerical model previously described, a computer calculation program was written with MATLAB R2014a software for the calculation of the eigenvalues and eigenvectors for the loss in stability of a slender bar subjected to compression. Because just the first value of the critical buckling force and the corresponding deformed shape were of interest, the calculation program reported just the first three values of the critical buckling force and plotted the corresponding eigenvectors.

The assumptions considered in the numerical analysis of the finite elements were the following: (i) the material of the bar was isotropic, homogeneous, and linearly elastic; (ii) the hypothesis of small strains was valid; and (iii) the normal stress *σ<sup>p</sup>* at the proportionality limit was approximately equal to the normal stress at yielding, denoted with *fy*, for the material of the bar.

In order to obtain the numerical solution with the finite element method, the main steps covered by the MATLAB program were the following: (i) meshing of the bar in a certain number of finite elements; (ii) computing both the elastic stiffness matrix [*Ke*] and the geometric stiffness matrix *Kge* corresponding to each finite element according to Equations (43) and (44), respectively; (iii) assembling all the stiffness matrices [*Ke*] and geometric stiffness matrices *Kge* in order to obtain the elastic stiffness matrix [*K*] and geometric stiffness matrix *Kg* of the column analyzed; (iv) computing the eigenvalues *λ<sup>k</sup> k* = 1, *n* by solving Equation (42); (v) computing the eigenvectors {*Uk*} *k* = 1, *n* using Equation (40); and (vi) by using the eigenvalues *λ*k, computing the critical buckling forces *Pcr*. The smallest critical buckling force corresponded to the minimum eigenvalue *λ*min.

In order to validate the numerical model, the loss in stability was analyzed using the numerical model with 18 finite elements for a bar having a pin connection at one end and a simple support at the other end, for which the geometrical characteristics and the material properties are given in Table 1. Considering the geometrical characteristics of the bar given in Table 1, the following quantities were computed: the area *A* of the cross-section of 2826 mm2; the second moment of inertia *I* of the cross-section, whose value was 2,896,650 mm4; and the radius *i* of inertia, having a value of 32.01562 mm.


**Table 1.** Geometrical characteristics and material properties for the bar analyzed by numerical modeling.

For the bar involved, the slenderness ratio *λ* was computed with Equation (45):

$$
\lambda = l/i = 8000/32.01562 = 249.878.\tag{45}
$$

For steel of type S355, whose properties are shown in Table 1, the slenderness ratio *λ*<sup>0</sup> that limited buckling in the elastic field was computed with Equation (46):

$$
\lambda\_0 = \pi \sqrt{\mathcal{E}/\sigma\_p} = \pi \sqrt{\mathcal{E}/f\_y} = 76.37041,\tag{46}
$$

where it is assumed that the normal stress *σ<sup>p</sup>* at the proportionality limit is approximately equal to the normal stress at yielding, denoted with *fy*, for the material of the bar (Table 1).

The first three eigenshapes and the corresponding eigenvalues for the bar analyzed are shown in Figure 3.

**Figure 3.** The first three shapes of stability loss and corresponding critical buckling forces for the bar, whose geometrical and material characteristics are given in Table 1.

In Figure 4, the convergence of the critical buckling load *Pcr* obtained with the algorithm of the numerical model is shown related to the number of the finite elements of the

numerical model, and it is analyzed with respect to the critical buckling load of 93.807 kN computed with Equation (24) using the analytical model. By analyzing Figure 4, it can be remarked that the solutions obtained through numerical modeling with the FEM tended asymptotically to the value of 93.86 kN for the numerical model, which had at least six elements. It can be concluded that the numerical model consisting of 18 finite elements provided results that were sufficiently accurate concerning the critical buckling load.

**Figure 4.** Analysis concerning the convergence of the solution obtained for the critical buckling load by the FEM compared with respect to the value computed with the analytical model.

Six cases of bars with pin connections at one end and simple supports at the other end, shown in Table 2, were analyzed using the numerical model with finite elements in order to show the effects of a stepwise variable cross-section on the critical force of stability loss.

**Table 2.** Values of the ratios *k*<sup>1</sup> and *k*<sup>2</sup> considered for the analyzed cases with the FEM for bars with pin connections at one end and simple supports at the other end.


\* Ratios *k*<sup>1</sup> and *k*<sup>2</sup> are shown in Figure 1a.

In Table 2, two extreme cases were considered for bars whose the second moments of inertia were constant along the bar length *l*: (i) a bar having a second moment of inertia equal to *I* (code CONST\_1I); and (ii) a bar having a second moment of inertia equal to 4*I* (code CONST\_1I). For both cases, the length *l* of the bar was that given in Table 1, and it was assumed that the bar had an annular cross-section. The value *I* for the second moment of inertia corresponded to an annular cross-section having an inner diameter *d* and an outer diameter *D*, which are also given in Table 1. For the second bar, for which the second moment of inertia was equal to 4*I*, the inner diameter *d* and outer diameter *D* were computed by considering the area *A* of the annular cross-section in order for the normal stress in compression to be equal to the normal stress *σ<sup>p</sup>* at the proportionality limit (which was approximately equal to the normal stress at yielding *fy*).

The results obtained using the numerical model with 32 finite elements were comparatively analyzed for all the types of columns shown in Table 2.

#### *2.3. Numerical Modeling and Simulation for Loss in Stability of a Column with Pin Connections at Ends with Continuous Variable Cross-Section*

2.3.1. Cases Approached and Assumptions concerning the Variation in the Second Moment of Inertia

Four cases of variation in the second moment of inertia *I*(*x*) were investigated: (i) parabolic variation; (ii) sinusoidal variation; (iii) triangular variation; (iv) trapezoidal variation; and (v) constant variation.

It was assumed that the arbitrary cross-section of a bar whose area was *A*(*x*) was an annular cross-section having an inner diameter *d*(*x*) and an outer diameter *D*(*x*). The length of the bar was denoted with *l*. Assuming that the bar must lose its stability in the elastic field (for which Euler's relation is valid), the normal stress *σ* developed at the arbitrary point of the cross-section must be smaller or equal to the normal stress at the proportionality limit denoted with *σp*, which was considered to be equal to the normal stress at yielding *fy* (given in Table 1 for S355 steel). In this context, the area *A*(*x*) of the cross-section was computed with Equation (47), considering that the normal stress *σ* developed at that cross-section was equal to the normal stress *σ<sup>p</sup>* at the proportionality limit for the material of the bar at stability loss:

$$A(\mathbf{x}) = \mathbf{P}\_{cr} / \sigma\_{p\prime} \tag{47}$$

where *Pcr* is the critical buckling force (or the critical force of stability loss).

For the arbitrary cross-section of the bar, knowing the area *A*(*x*) computed with Equation (47) and the variation in the second moment of inertia *I*(*x*), the inner and outer diameters of the cross-section were computed with Equations (48) and (49), respectively [20]:

$$d(\mathbf{x}) = \sqrt{\frac{8I(\mathbf{x})}{A(\mathbf{x})} - \frac{2A(\mathbf{x})}{\pi}},\tag{48}$$

$$D(\mathbf{x}) = \sqrt{\frac{8I(\mathbf{x})}{A(\mathbf{x})} + \frac{2A(\mathbf{x})}{\pi}}.\tag{49}$$

Equations (48) and (49) show that the wall thickness of the bar cross-section changed along the axis of the bar as long as the ratio between *d*(*x*) and *D*(*x*) was not constant.

To analyze the critical buckling force *Pcr* for a column with pin connections at both ends whose second moment of inertia *I* had continuous variation along the axis of the bar, numerical models were used that were similar to the one corresponding with the bar with stepwise variation in the cross-section. In order to analyze the solution convergence for the critical buckling force *Pcr*, the numerical analysis was repeated for FEMs consisting of 20, 80, 100, 200, 400, and 800 elements, respectively. The convergence analysis made for FEMs corresponding to bars with continuous variation in the second moment of inertia along the bar axis led to the conclusion that the solutions for the critical buckling force converged for FEMs consisting of a minimum of 100 elements.

2.3.2. Parabolic Variation in the Second Moment of Inertia of the Cross-Section along the Bar

The consideration of the second moment of inertia *I*(*x*) of an arbitrary cross-section with parabolic variation along the bar axis is given by Equation (50):

$$I(\mathbf{x}) = A\mathbf{x}^2 + B\mathbf{x} + \mathbf{C},\tag{50}$$

where *x* represents the position of the arbitrary cross-section with respect to the pinconnected end of the bar.

The constants *A*, *B*, *C* were computed considering that the second moment of inertia was equal with *I* for the cross-sections of the bar ends, and it was equal to 4*I* for the middle of the bar. These conditions are given by Equation (51):

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{0}, & \mathcal{I}(\mathbf{0}) = \mathbf{C} = \mathbf{I};\\ \text{for } \mathbf{x} = \mathbf{l}/2, & \mathcal{I}(\mathbf{l}/2) = A(\mathbf{l}/2)^2 + B\mathbf{l}/2 + \mathbf{C} = 4\mathbf{I};\\ \text{for } \mathbf{x} = \mathbf{l}, & \mathcal{I}(\mathbf{l}) = A\mathbf{l}^2 + B\mathbf{l} + \mathbf{C} = \mathbf{I}; \end{cases} \tag{51}$$

and led to the constants given by Equation (52):

$$A = -12I/l^2; \quad B = 12I/l; \quad C = I. \tag{52}$$

By replacing the above constants *A*, *B*, and *C* in Equation (50), the parabolic variation in the second moment *I*(*x*) along the axis of the bar is given by Equation (53):

$$I(\mathbf{x}) = -\left(12I/l^2\right)\mathbf{x}^2 + (12I/l)\mathbf{x} + I,\tag{53}$$

which is plotted in Figure 5 for a bar whose length *l* is given in Table 1 and whose second moment of inertia *I* of the cross-sections located at the bar ends was computed for an annular cross-section having an inner diameter *d* and an outer diameter *D*, which are also given in Table 1.

**Figure 5.** Parabolic variation in the second moment of inertia *I* along the bar axis.

In Figure 6, a geometrical model is shown of a bar whose parabolic variation in the second moment of inertia *I*(*x*) is given by Equation (53).

**Figure 6.** Geometrical model for a bar with parabolic variation in the second moment of inertia along the bar axis: (**a**) isometric view and (**b**) longitudinal section.

2.3.3. Sinusoidal Variation in the Second Moment of Inertia of the Cross-Section along the Bar

The variation in the second moment of inertia *I*(*x*) of an arbitrary cross-section along the bar axis as a sinusoidal function is given by Equation (54):

$$I(x) = 3I\sin(\pi x/l) + I.\tag{54}$$

Considering the conditions given by Equation (55), it could be checked that the second moment of inertia was equal to *I* for the cross-sections located at the bar ends, and it was 4*I* for the cross-section located at the middle of the bar:

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{0}; \qquad I(0) = 3I\sin(\pi \cdot 0/l) + I = I;\\ \text{for } \mathbf{x} = l/2; \quad I(l/2) = 3I\sin(\pi(l/2)/l) + I = 3I\sin(\pi/2) + I = 4I;\\ \text{for } \mathbf{x} = l; \qquad I(l) = 3I\sin(\pi l/l) + I = 3I\sin\pi + I = I. \end{cases} \tag{55}$$

The sinusoidal variation in the second moment *I*(*x*) along the axis of the bar given by Equation (54) is plotted in Figure 7 for a bar whose length *l* is given in Table 1 and whose second moment of inertia *I* of the cross-sections located at bar ends was computed for an annular cross-section having an inner diameter *d* and an outer diameter *D*, which are also given in Table 1.

**Figure 7.** Sinusoidal variation in the second moment of inertia along the bar axis.

In Figure 8, a geometrical model is shown of a bar whose sinusoidal variation in the second moment of inertia *I*(*x*) is given by Equation (54).

**Figure 8.** Geometrical model for s bar with sinusoidal variation in the second moment of inertia along the bar axis: (**a**) isometric view and (**b**) longitudinal section.

2.3.4. Triangular Variation in the Second Moment of Inertia of the Cross-Section along the Bar

It was considered that the second moment of inertia linearly increased from the value *I* to the value 4*I* from the pin-connected end of the bar to the middle of the bar, which meant position *x* of the cross-section was in the range of [0, *l*/2]. Then, the second moment of inertia linearly decreased from the value 4*I* to the value *I* from the middle of the bar to the other end of the bar, which meant a variation of *x* in the range of [*l*/2, *l*].

A linear function for the second moment of inertia for the first half of the bar, which meant *x* ∈ (0, *l*/2), and it was given by Equation (56):

$$I\_1(\mathbf{x}) = A\mathbf{x} + B\_\prime \tag{56}$$

whose constants *A* and *B* are computed using the conditions given by Equation (57) regarding the values of the second moments of inertia for the cross-sections located at the pin-end connection and at the middle of the bar:

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{0}; \qquad I\_1(\mathbf{0}) = A \cdot \mathbf{0} + B = I;\\ \text{for } \mathbf{x} = \mathbf{l}/2; \quad I\_1(\mathbf{l}/2) = A \mathbf{l}/2 + B = 4I. \end{cases} \tag{57}$$

Using the conditions given by Equation (57), the constants *A* and *B* were computed, and the results are given in Equation (58):

$$A = \text{6I}/\text{l}; \quad B = I. \tag{58}$$

By replacing the constants *A* and *B* in Equation (56), a linear function was obtained for the second moment of inertia for the first half of the bar, which is given in Equation (59):

$$I\_1(\mathbf{x}) = \mathbf{6}I\mathbf{x}/l + l, \quad \text{for } \mathbf{x} \in (0, \, 0.5l). \tag{59}$$

A linear function was assumed for the second moment of inertia for the second half of the bar, which meant *x* ∈ (0.5, *l*), and it was given by Equation (60):

$$I\_2(\mathbf{x}) = \mathbb{C}\mathbf{x} + D\_\prime \tag{60}$$

whose constants *C* and *D* are computed by using the conditions given by Equation (61) regarding the values of the second moments of inertia for the cross-sections located at the middle of the bar and at bar end, which is simple-supported:

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{l}/2; & \mathbf{l}\_2(\mathbf{l}/2) = \mathbf{C}l/2 + D = 4\mathbf{l};\\ \text{for } \mathbf{x} = \mathbf{l}; & \mathbf{l}\_2(\mathbf{l}) = \mathbf{C}l + D = I. \end{cases} \tag{61}$$

Using the conditions given by Equation (61), the constants *C* and *D* were computed, and the results are given in Equation (62):

$$\mathbf{C} = -6\mathbf{I}/\mathbf{l}; \quad \mathbf{D} = 7\mathbf{I}.\tag{62}$$

By replacing the constants *C* and *D* in Equation (60), a linear function was obtained for the second moment of inertia for the first half of the bar, which is given in Equation (63):

$$I\_2(\mathbf{x}) = -6I\mathbf{x}/l + 7I,\quad \text{for } \mathbf{x} \in (0.5l, l). \tag{63}$$

The triangular variation in the second moment *I*(*x*) along the axis of the bar given by Equations (59) and (63) is plotted in Figure 9 for a bar whose length *l* is given in Table 1 and whose second moment of inertia *I* of the cross-sections located at the bar ends was computed for an annular cross-section having an inner diameter *d* and an outer diameter *D*, which are also given in Table 1.

**Figure 9.** Triangular variation in the second moment of inertia along the bar axis.

In Figure 10, a geometrical model is shown of a bar whose triangular variation in the second moment of inertia *I*(*x*) is given by Equations (59) and (63).

**Figure 10.** Geometrical model for a bar with triangular variation in the second moment of inertia along the bar axis: (**a**) isometric view and (**b**) longitudinal section.

2.3.5. Trapezoidal Variation in the Second Moment of Inertia of the Cross-Section along the Bar

For trapezoidal variation in the second moment of inertia *I*(*x*) along the axis of a bar, three portions of the bar length *l* were considered in order to obtain the functions for such a variation: (i) for the first third of the bar, which meant *x* (0, *l*/3), the second moment of inertia linearly increased from the value *I* to the value 4*I*; (ii) for the second third of the bar, which meant *x* (*l*/3, 2*l*/3), the cross-section remained constant, and the second moment of inertia was equal with 4*I*; and (iii) for the last portion, which meant *x* (2*l*/3, *l*), the second moment of inertia linearly decreased from the value 4*I* to the value *I*.

The linear function corresponding to the first portion is given in Equation (64):

$$I\_1(\mathbf{x}) = A\mathbf{x} + B\_\prime \tag{64}$$

for which the constants *A* and *B* are computed using the conditions written in Equation (65):

$$\begin{cases} \text{for } \mathbf{x} = \mathbf{0}; \qquad I\_1(\mathbf{0}) = A \cdot \mathbf{0} + B = I;\\ \text{for } \mathbf{x} = \mathbf{l}/3; \quad I\_1(\mathbf{l}/3) = A \mathbf{l}/3 + B = 4I. \end{cases} \tag{65}$$

Solving the above system of two equations led to the following values for the constants *A* and *B*:

$$A = \Theta I/l; \quad B = l. \tag{66}$$

By replacing the constants *A* and *B* in Equation (64), a linear function was obtained of the second moment of inertia corresponding to the first portion of the bar, as shown in Equation (67):

$$I\_1(\mathbf{x}) = \theta I \mathbf{x} / l + I, \quad \text{for } \mathbf{x} \in [0, l/\\$]. \tag{67}$$

The function for the second moment of inertia corresponding to the second portion of the bar is given by Equation (68):

$$I\_2(\mathbf{x}) = 4I, \quad \text{for } \mathbf{x} \in \left[l/3, 2l/3\right]. \tag{68}$$

For the third portion of the bar, the function of the second moment of inertia was assumed, as given by Equation (69):

$$I\_3(\mathbf{x}) = \mathbb{C}\mathbf{x} + D\_\prime \tag{69}$$

for which constants *C* and *D* are computed using the following conditions:

$$\begin{cases} \text{for } \mathbf{x} = 2l/3; & l\_3(l/3) = 2\mathbf{C}l/3 + D = 4l; \\ \text{for } \mathbf{x} = l; & l\_3(l) = \mathbf{C}l + D = I. \end{cases} \tag{70}$$

By solving the system of two equations given by Equation (70) and replacing constants *C* and *D* in Equation (69), the following function of the second moment of inertia was obtained:

$$I\_3(\mathbf{x}) = -\Re \mathbf{l} \mathbf{x} / l + 10I\_\prime \quad \text{for } \mathbf{x} \in [2l / \Im \mathbf{\hat{z}}, l]. \tag{71}$$

Considering the functions given by Equations (67), (68), and (71), trapezoidal variation in the second moment of inertia along the bar axis is graphically shown in Figure 11 for a bar whose length *l* is given in Table 1 and whose second moment of inertia *I* of the cross-sections located at the bar ends was computed for an annular cross-section having an inner diameter *d* and an outer diameter *D*, which are also given in Table 1. The geometrical model of such a bar is shown in Figure 12.

**Figure 11.** Trapezoidal variation in the second moment of inertia along the bar axis.

**Figure 12.** Geometrical model for a bar with trapezoidal variation in the second moment of inertia along the bar axis: (**a**) isometric view and (**b**) longitudinal section.

#### **3. Results**

#### *3.1. Results Obtained by Analytical Model*

Figure 13 shows the variation in the ratio *Pcr*/*Pcr*0. computed with Equation (25) for the case of *k*<sup>1</sup> = 4. related to the value *k*<sup>2</sup> ∈ [0; 3]. It may be observed that, for this design of a column with a variable cross-section, the maximum value of the ratio *Pcr*/*Pcr*0. was equal to 3.69 in the case when *k*2. was equal to 3.

**Figure 13.** Variation in the normalized critical buckling force *Pcr*/*Pcr*0. computed with Equation (25) considering *k*<sup>1</sup> = 4 related to the ratio *k*2.

Figure 14 shows the variation in the rationality factor *krat* computed with Equation (37) for the case of *k*<sup>1</sup> = 4. related to the value *k*<sup>2</sup> ∈ [0; 3]. It was observed that the column with a stepwise variable cross-section was more rationally designed because the rationality factor *krat*. was always greater than 1, and the greatest value was 2.1086 in the case when *k*2. was equal to 3. In fact, for the ratio *k*2. in the range of [2, 3], the rationality factor *krat*. varied between 1.9973 and 2.1086.

**Figure 14.** Variation in the rationality factor *krat*. computed with Equation (37) considering *k*<sup>1</sup> = 4. related to the ratio *k*2.

The least squares method was used for the approximation of the data in both Figures 13 and 14 considering the second-degree polynomial functions. The approximation functions both for the normalized critical buckling load *Pcr*/*Pcr*<sup>0</sup> and for the rationality factor *krat*. are given in Figures 13 and 14, where the value *R*<sup>2</sup> close to 1 shows that the data were accurately approximated.

#### *3.2. Results by Numerical Modeling for Loss in Stability of a Bar with Pin Connections at Ends with Stepwise Variable Cross-Section*

Considering the column with pin connections at its ends and a stepwise variable crosssection shown in Figure 1, the results obtained using the MATLAB calculation program for the numerical model with 32 finite elements are plotted in Figure 15 for all six cases given in Table 2.

**Figure 15.** The first three shapes of stability loss and corresponding critical forces *Pcr*. obtained by the FEM for the six cases of bars analyzed: (**a**) CONST\_1I; (**b**) STEPWISE405; (**c**) STEPWISE410; (**d**) STEPWISE420; (**e**) STEPWISE430; and (**f**) CONST\_4I (details about each case are given in Table 2).

#### *3.3. Validation of the Numerical Model by Theoretical Results for a Column with Stepwise Variable Cross-Section*

Table 3 shows the results obtained for the critical buckling force *Pcr*. using both the FEM and the analytical model in the cases of columns with stepwise variable crosssections and with constant cross-sections involved in this research. It was observed that the numerical model was validated by the results obtained with the analytical model because the maximum error was equal to 3.84%, as shown in Table 3.


**Table 3.** Comparison between the critical buckling force *Pcr* obtained by the FEM with the one obtained with the analytical model for different cases of bars involved in this research.

\* Details about bars of each case are given in Table 2.

The values of the normalized critical buckling forces obtained by the FEM and by the analytical model are comparatively plotted in Figures 16 and 17, respectively, for the columns with stepwise variable cross-sections and with constant cross-sections. These graphs also show a good correlation between the numerical model and the analytical model.

**Figure 16.** Comparison of the normalized critical buckling forces *Pcr*/*Pcr*0. computed with the finite element models for the bars involved (details about each case are given in Table 2).

**Figure 17.** Comparison of the normalized critical buckling forces *Pcr*/*Pcr*<sup>0</sup> computed with the analytical model for the bars involved in this research (details about each case are given in Table 2).

*3.4. Results of Numerical Modeling for Loss in Stability of a Column with Pin Connections at Ends with Continue Variable Cross-Section*

3.4.1. Critical Buckling Forces for the Parabolic Variation in the Second Moment of Inertia of the Cross-Section along the Bar

Using the numerical model with 100 finite elements for the bar with parabolic variation in the second moment of inertia, the first three eigenshapes of stability loss were obtained that corresponded to the first three values for the critical buckling force shown in Figure 18. It was observed that the smallest value of the critical force for stability loss *Fcr*. was equal to 329.7439 kN.

**Figure 18.** The first three shapes of stability loss and corresponding critical buckling forces for the bar with parabolic variation in the second moment of inertia along the bar axis.

3.4.2. Critical Buckling Forces for the Sinusoidal Variation in the Second Moment of Inertia of the Cross-Section along the Bar

Using the numerical model with 100 finite elements, for the bar with sinusoidal variation in the second moment of inertia, the first three eigenshapes of stability loss were obtained that corresponded to the first three values for the critical buckling force shown in Figure 19. It was observed that the smallest value of the critical force for stability loss *Fcr* was equal to 321.6489 kN.

3.4.3. Critical Buckling Forces for the Triangular Variation in the Second Moment of Inertia of the Cross-Section along the Bar

Using the numerical model with 100 finite elements for the bar with triangular variation in the second moment of inertia, the first three eigenshapes of stability loss were obtained that corresponded to the first three values for the critical buckling force shown in Figure 20. It was observed that the smallest value of the critical force for stability loss *Fcr*. was equal to 275.4322 kN.

**Figure 19.** The first three shapes of stability loss and corresponding critical buckling forces for the bar with sinusoidal variation in the second moment of inertia along the bar axis.

**Figure 20.** The first three eigenshapes of stability loss and corresponding values for the critical buckling forces for the bar with triangular variation in the second moment of inertia along the bar axis.

3.4.4. Critical Buckling Forces for the Trapezoidal Variation in the Second Moment of Inertia of the Cross-Section along the Bar

Using the numerical model with 100 finite elements for the bar with trapezoidal variation in the second moment of inertia, the first three shapes of stability loss corresponding to the first three values of the critical buckling force were obtained and are shown in Figure 21. These were obtained using the numerical model with finite elements of the bar

with trapezoidal variation in the second moment of inertia along the bar axis. The smallest value of the critical force for stability loss *Fcr*. was equal to 333.7152 kN.

**Figure 21.** The first three shapes of stability loss and corresponding values for the critical buckling force for the bar with trapezoidal variation in the second moment of inertia along the bar axis.

#### **4. Discussion**

In Figure 22, the values of the normalized critical buckling forces *Pcr*/*Pcr*<sup>0</sup> are comparatively analyzed for all the bars involved in this theoretical research in order to establish the best shapes of both bars with stepwise variable cross-sections and with continuous variable cross-sections, which can lead to a significant increase in the critical buckling force.

**Figure 22.** Comparison of the critical buckling force obtained by the FEM for all the cases of bars involved in this research.

By analyzing the results shown in Figure 22, it can be observed that the normalized critical buckling force of 3.43 obtained for the bar with stepwise variation in the crosssection corresponding to case STEPWISE420 was approximately equal to the normalized

critical buckling force of 3.427 obtained for the bar with sinusoidal variation in the second moment of inertia along the bar axis from *I* to the maximum value of 4*I*. This remark is very important as long as the bar with the stepwise variable cross-section is more easily obtained from a technological point of view.

Considering Figure 22, it may be also remarked that the normalized critical buckling force of 3.556 recorded for the bar with trapezoidal variation in the second moment of inertia along the bar axis was close to the normalized critical buckling force of 3.513 obtained for the bar with parabolic variation in the second moment of inertia along the bar axis.

For the bars with variable cross-sections involved in this research, the highest value of 3.69 for the normalized critical buckling force was recorded for the bar with stepwise variation in the cross-section corresponding to case STEPWISE430 (Figure 22).

#### **5. Conclusions**

The research presented in this paper is very important for the design of shapes of slender bars that can lose their stability under compression loads. Increasing the critical buckling force for columns with annular cross-sections by stepwise or continuous variation in the cross-sections along the bar axes was proposed. Moreover, the shapes of the slender bars were designed so that these bars lost their stability in the elastic field.

A particular case was considered for which the second moment of inertia varied between the values *I* and 4*I* along the bar axis. For bars with variation in the cross-sections in three steps, the normalized critical buckling force continuously increased by a seconddegree polynomial function reported in the paper related to the ratio *k*<sup>2</sup> between the lengths of the bar portions. It was shown that the geometry proposed for the bar with stepwise variation in the cross-section was more rational with respect to the bar with a constant cross-section, taking into account the ratio between the critical buckling force and the volume of the bar.

Among the bars involved in this research with continuous variations in the annular cross-sections, it was concluded that the parabolic or trapezoidal variations in the crosssections were the best solutions. These design solutions led to increases in the critical buckling force by 3.513 times and 3.556 times, respectively, compared with the bar with a constant cross-section along the axis. In fact, the normalized critical buckling forces obtained for both parabolic and trapezoidal variations were close to the one corresponding to case STEPWISE420 with stepwise variation in the cross-section, which was more technologically affordable.

The numerical analysis approach in this research was one of an appropriate design solution to optimize the variation in the cross-section along the axis of a bar by comparing the results in order to maximize the critical buckling force. The numerical model was validated by an analytical model for the bar with stepwise variation in the cross-section, and then the numerical approach was applied for bars with different methods of variation in the second moment of inertia along the bar axis. In fact, the continuous variation in the cross-section was approximated with stepwise variation in the cross-section considering that the bar consisted of a number of portions equal to the number of finite elements considered in the numerical analysis. The main advantage was that the MATLAB program used in this research was not a commercial one, and as a result, the algorithm could be further adapted for other types of variation in the cross-section along the bar axis.

At a time when the rapid printing of construction elements has taken off around the world in the field of construction, the results presented in this paper are of great interest in many applications, taking into account that these results could lead to increasing the critical buckling force for compressed slender bars. However, there is still a need for 3D-printing equipment for construction to be perfected, developed, and made more affordable in terms of cost so that the manufacture of bars with variable cross-sections is easy for construction elements whose shapes have been optimized.

Considering the design solutions for continuous variation in the cross-section presented in this article, it is possible to manufacture sleeves to be welded on compressed bars

of various structures (stadium roofs, steel bridge structures, trusses) to increase the critical buckling force. The sleeves may be welded without removing the bars from the structures. Because the research in this article was limited to slender bars with pin connections at both ends, the study may be continued for other boundary conditions in further research.

**Author Contributions:** Conceptualization, M.F.B.; formal analysis, M.F.B. and C.C.; investigation, M.F.B. and C.C.; methodology, M.F.B.; supervision, M.F.B.; validation, M.F.B. and C.C.; visualization, M.F.B. and C.C.; writing—original draft, M.F.B. and C.C.; writing—review and editing, C.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors acknowledge the support of Project No. 11002100 from the Transylvania University of Brasov and a construction company. The authors also acknowledge the support of the Transilvania University of Brasov for providing the software involved and for financial support regarding the publication fees of this article.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge the support of Project No. 11002100 from the Transylvania University of Brasov and a construction company. The authors also acknowledge the support of the Transilvania University of Brasov for providing the software involved and for financial support regarding the publication fees of this article.

**Conflicts of Interest:** The authors declare no conflict of interest.
