Exact and Approximate Stiffness Matrix and Nodal Load Vector for a Beam Finite Element with Linearly Varying Stiffness along Its Length

Introduction. In modern engineering practice, specifically in the context of sustainable development and resource optimization, structures with variable stiffness — for example, beams with variable cross-sections — are widely used [1]. Such solutions significantly reduce the material consumption and dead weight of structures, which directly improves their cost-effectiveness and environmental friendliness by reducing material consumption [2]. In addition, the mathematical model of a beam of variable stiffness serves as a tool for the qualitative description of the stress-strain state of reinforced concrete elements, in which the change in stiffness occurs due to the nonlinear behavior of the material [3][4] and the formation of cracks [5][6].

In statically indeterminate beams, which are common in structural mechanics, the distribution of internal forces depends on the stiffness ratio of the elements. Neglecting the variable nature of stiffness can cause significant errors in determining internal forces and, consequently, the incorrect cross-section selection, which compromises both the reliability and economic viability of the project. Figure 1 shows a clear example: taking into account variable stiffness in a beam (Fig. 1 b) has led to a radical change in the bending moment diagram compared to the calculation model of constant stiffness (Fig. 1 a).

The analysis of beams with variable cross-sections is traditionally performed using numerical methods, in particular, the finite element method (FEM). Despite the existence of analytical solutions [7], they are often represented as complex series [8] or special functions [9], which complicates their practical implementation in engineering software. Moreover, such solutions often lack universality for arbitrary laws of stiffness change and boundary conditions [10][11].

The stress-strain state of beams with variable cross-sections has been studied by various authors for a long time. In the context of numerical calculations using FEM, some of the first works [12][13] can be mentioned. The most common approach to accounting for variable stiffness is a discrete approximation using finite elements of constant stiffness (Fig. 1 b), whose value is taken equal to the average in the section [14][15]. Despite its simplicity, this approach does not always ensure the required accuracy of calculations.

More accurate results are obtained by using a linear approximation of the stiffness within the finite element. This approach is described in this article. At the same time, classical methods for constructing stiffness matrices for elements with variable parameters [16], based on variational approaches [17] or approximations [18], may not provide the exact fulfillment of the differential equations of equilibrium and boundary conditions. This causes errors with significant changes in stiffness along the length of the element or in extreme cases, for example, when modeling zones with a sharp drop in stiffness (formation of cracks, plastic hinge).

Thus, the relevance of this study is determined by the need to develop an efficient and more accurate method for calculating beams with variable stiffness.

The objective of this study is to obtain and verify stiffness matrices for a beam finite element with a piecewise linear law of stiffness change.

To reach this objective, the following tasks were formulated:

1. To obtain an approximate stiffness matrix on the basis of a variational formulation of the problem and a linear approximation of the stiffness within the finite element.
2. To develop an exact stiffness matrix for a beam finite element with linearly varying stiffness using an analytical solution to the differential bending equation.
3. To conduct a comparative analysis of the accuracy of the calculation results using the proposed stiffness matrices and the classical approach with discrete stiffness approximation on a series of test examples; to determine the accuracy by comparing to the results of the analytical solution.

Materials and Methods

Variational Method for Obtaining the Stiffness Matrix. Variational methods [19] are the generally accepted and most versatile methods for obtaining the stiffness matrix of a finite element. The accuracy of the resulting matrix depends on how well the input approximation of the desired function reflects the solution to the differential equation.

For the beam under consideration, the classical Euler-Bernoulli hypotheses are applied. The change in the flexural stiffness of the cross-section D_z along the length of the finite element is assumed to be a linear function:

(1)

where D₀ — initial stiffness; L — length of the finite element; ψ₀ and ψ₁ — coefficients that take into account the change in initial stiffness D₀ at the beginning and end of the element, respectively. In this case, ψ₀ ≥ 0 and ψ₁ ≥ 0.

The finite element method is considered in the form of displacements. The unknowns are the deflection and rotation angle of the section v_i and φ_i (Fig. 2 a) at the beginning x = 0 and end x = L of the element.

To approximate the displacement function v(x), Hermite polynomials are used [1][4]:

(2)

where N_i — shape functions (Hermite polynomials) written below:

(3)

Approximation (2) corresponds to the following boundary conditions:

(4)

The prime in (4) denotes the derivative of function v(x) in respect to x.

From formulas (3), it follows that classical shape functions do not take into account the distribution of bending stiffness along the length of the element, which, in turn, negatively impacts the accuracy of the calculation results. To improve the accuracy of calculations, when using a variational formulation of the problem, the double approximation method can be applied [20].

The bending moment is associated with the curvature by the relation:

(5)

Here, D_z — bending stiffness of the beam cross-section according to (1).

Taking into account (2) and (3):

(6)

The potential energy of bending deformation is determined by the formula:

(7)

Taking into account (6) and (1), formula (7) takes the form:

(8)

where

(9)

If we introduce the vector of unknowns , then (8) can be rewritten in the following form:

(10)

Here, [K] — stiffness matrix of the element, determined taking into account (8) and (9):

(11)

This case ψ₀ = ψ₀ = 1 corresponds to the classical beam stiffness matrix.

The positive directions of the nodal loads and internal forces are shown in Figure 2 b, from which it follows:

(12)

Let us calculate the work of external forces.

(13)

In matrix form

(14)

where  — vector of concentrated nodal loads;  — vector of nodal loads from distributed loads. The components of these vectors have the form (15):

(15)

As we can see, the vector of nodal forces from distributed loads obtained in (15) does not take into account the change in stiffness along the length. This is a drawback of the variational approach in the context of the problem under consideration.

The total energy of the system, taking into account (10) and (14), is determined by the expression:

(16)

Applying the principle of minimum potential energy, we obtain the condition for the stationarity of functional (16):

(17)

The basic equation of the finite element method follows from (17):

(18)

Direct Method for Obtaining the Stiffness Matrix. The resulting stiffness matrix (11) has a significant drawback. If we assume v₀ = φ₀ = 0, as well as ψ₀ = 0 and ψ₁ = 1 (which corresponds to zero stiffness at the left clamped edge), the moment at the left edge will be determined by the formula:

(19)

Note that (19) considers a limiting case — a mathematical idealization, where the stiffness at the support is not simply small but asymptotically tends to zero. Thus, in accordance with (19), the moment at the left edge does not vanish, which contradicts formula (5). The reason lies in the inaccurate approximation of the deflection function by formula (2). To eliminate this discrepancy, an exact stiffness matrix was found through direct integrating the differential equation for beam bending, taking into account the change in stiffness along the length by formula (1).

A similar approach was used by other authors. In [21], the solution is obtained in the form of power series, therefore, it cannot be considered closed. In [22], a direct integration method is used to obtain the stiffness matrix of an element with a power law of stiffness change for stability problems. In [23], an arbitrary law of stiffness change along the length is considered. However, for the practical use of the presented results, integration is required. Works [24][25] are devoted to solving stability problems. In these papers, the law of stiffness change along the length is given in the form of a polynomial of the second and fourth degree, while the vector of nodal forces from distributed loads is not provided. A common drawback of works that take into account second and higher order polynomials in the function of stiffness change along the length is the complexity of calculating the stiffness matrix. For example, according to [26], to obtain the stiffness matrix coefficients, it is required to calculate expressions of the form sin[ε ln (λ)] (the original notations are preserved), which complicates the analysis of limiting cases and the implementation of solutions in commercial software packages.

To obtain an accurate stiffness matrix, the differential equation for plane transverse bending of a beam with variable stiffness, known from the course on strength of materials, is considered:

(20)

Given the form (1), the general solution to the homogeneous equation (20) will contain a logarithm. We eliminate the dimensional quantities under the logarithm sign by replacing the variable using the formula:

(21)

Considering (21):

(22)

As a result, general solution [27] of equation (20) takes the form:

(23)

where α = ψ₀ / (ψ₀ – ψ₁), v^* — particular solution to the inhomogeneous equation. Due to its cumbersome nature, the particular solution is not given.

Solution (23) is valid for the cases ψ₀ ≠ ψ₁, ρ ≠ α, since these cases generate uncertainties that are resolved through limiting in formulas (31), (32) and (33).

The unknown integration constants are determined from the boundary conditions:

(24)

In (24), the prime denotes the derivative with respect to ρ. Multiplier L at the rotation angles in the boundary conditions (24) is present due to the introduced change of variable (21), which is reflected in (22). After determining the integration constants from (24), the moment and shear force in the beam, taking into account the introduced change of variables, are calculated using the formulas:

(25)

The coefficients of the stiffness matrix, taking into account (25) and (12), can be found from the expressions:

. (26)

Matrix form (26) is equivalent to the one obtained earlier (18). The stiffness matrix in the case under consideration will take the form:

(27)

In (27), the following notations are introduced for the sake of compactness:

(28)

To implement the proposed stiffness matrix in software, it is critical to be able to transform distributed loads into a nodal force vector. The exact expressions for vector are given below.

(29)

The main advantage of the obtained matrix (27) and vector (29) is the simplicity of implementation, as well as the absence of series and special functions in the calculations.

When solving certain problems, it is more convenient to relate section stiffness coefficients ψ₀ and ψ₁ to the curvature. In this case, it is necessary to know how the curvature is related to the displacements and rotation angles of the nodes of the element. Formula (5) can be used for this. Knowing how the moment is expressed at the beginning and end of the element according to formulas (26), we obtain:

(30)

Cases ψ₀ = 0, ψ₁ = 0 and ψ₀ = ψ₁ = ψ generate uncertainties of various types in matrix (27) and vector (29). Let us develop them.

Case ψ₀ = 0.

(31)

Case ψ₁ = 0.

(32)

From (31) and (32), it is evident that the columns and rows of the stiffness matrix, as well as the components of the nodal force vector from distributed loads associated with the bending moments at the left and right ends of the finite element, respectively, are set to zero — which corresponds to equation (5). The vanishing of the bending stiffness can occur, for example, when tensile reinforcement breaks, or the compressed zone of a reinforced concrete beam fails in any section. In the case of a statically indeterminate beam, the structure can continue to function during such failures. In commercial software packages, when such situations arise, hinges are typically activated in the calculation model.

Case ψ₀ = ψ₁ = ψ. This case corresponds to the classical stiffness matrix.

(33)

In practical calculations, to reduce computational errors that may accumulate during the calculation of the stiffness matrix components and factors (28), it is recommended to use formulas (31), (32) and (33) not when strictly reaching the limiting case, but when approaching it. For example, set the threshold conditions for using (27) and (29) in the form ψ_i > 10⁻⁴ and |1 – ψ₀ / ψ₁| > 10⁻³. The specified values are approximate and depend on the specific implementation features.

Figure 3 shows graphs of three elements of the stiffness matrix, calculated using formulas (11) and (27) for ψ₁ = D₀ = L = 1.

The results presented in Figure 3 show that, for ψ₀ = 0 and the corresponding boundary conditions, the stiffness matrix coefficients associated with the moment at the left edge vanish in the exact matrix (27) — this corresponds to the physical meaning. Approximate stiffness matrix (11) does not have this property. With a significant difference in stiffness, at the beginning and end of the element (ψ₀ / ψ₁ < 0.2), values of the coefficients of the exact and approximate matrices can differ by up to 100% (excluding from consideration the limiting case ψ₀ → 0). In such situations, it is recommended to use exact stiffness matrix (27). Similar conclusions can be drawn for the right edge of the finite element.

Research Results. To demonstrate the results obtained using the proposed stiffness matrices, calculation examples are performed for beams of length 2L, shown in Figure 4.

Due to symmetry with respect to x = L, half of the beam is considered in the solution. The boundary conditions under the action of a concentrated force (diagram in Fig. 4 a) are as follows:

(34)

Similarly for the case of distributed load (diagram in Fig. 4 b):

(35)

The following stiffness distribution function along the beam length is considered:

(36)

s₀ and s₁ — beam stiffness coefficients at the beginning and middle of the span, that is, D_z(0) = D₀s₀ and D_z(L) = D₀s₁. In this case, D_z(L / 2) = D₀. It is important to note that (36) determines the law of stiffness change along the length of the beam as a whole, while (1) — along the length of the finite element.

Function (36) is selected on the base of the method of load application and the boundary conditions. In this case, for a beam of constant stiffness, the moments at the support and in the middle of the beam reach maximum values. If the beam material is reinforced concrete, cracks and a decrease in flexural stiffness should be expected at these points [27]. At a quarter span, the moments are close to zero and, consequently, the stiffness is close to the initial D_z(L / 2) = D₀.

Figure 5 shows the graphs of the stiffness distribution along the length of the beam, as well as the approximation methods under consideration for two alternatives: option 1 — s₀ = 0.001, s₁ = 1.0; option 2 — s₀ = 0.6, s₁ = 0.2.

From the results presented in Figure 5, it is evident that the linear approximation, even with five finite elements, provides satisfactory agreement (the deviation does not exceed 5%) with the initial function of stiffness change along the length for the examples under consideration. High accuracy is maintained when approximating the stiffness with a piecewise linear function for variable cross-sections of complex shape [28]. The value of the average stiffness for discrete approximation was calculated using the formula:

(37)

The results of calculating the coefficients for the stiffness matrices of finite elements are given in Table 1. For discrete approximation, ψ₀ = ψ₁ = ψ.

Since an elastic problem is being solved, the results of calculating the stress-strain state components depend linearly on the load magnitude and initial stiffness D₀. Therefore, in the calculation, these quantities are assumed to be equal to unity.

The solution to equation (20) with boundary conditions (34) and the law of stiffness change (36) is found analytically for values s₀ and s₁ under consideration. Due to its cumbersome nature, this solution is not presented here. It was not possible to obtain an analytical solution for arbitrary values s₀ and s₁.

The results from the obtained analytical solution were then compared to the results of the finite element calculation. In the basic calculation, the beam was divided into five finite elements. To study the effect of the number of finite elements on the accuracy of the results, a series of calculations were performed in which their number varied.

Based on the calculation results, the error value was determined in comparison with the analytical solution. To formulate conclusions regarding the required number of elements, a 5% error threshold was adopted. Conclusions regarding the feasibility of applying this approach in practical calculations were based on the assumption that the threshold accuracy was achieved with fewer than 60 elements — in this case, the finite element length was 50 mm. Using a finer granularity is impractical for calculating beams in building structures — it significantly increases computational costs and calculation time. The maximum number of elements in the calculation examples was set to 400 (with a finite element length of 7.5 mm). With a larger number of finite elements, the solver implemented by the author began to accumulate computational errors.

Example 1. Solution to the problem of bending of a clamped beam loaded at the center with a concentrated force. A sketch of the calculation scheme is shown in Figure 4 a. The stiffness distribution is adopted according to option 1: s₀ = 0.001 and s₁ = 1. The graph of the stiffness function for this option is shown in Figure 5 a. Figures 6 and 7 show the results of calculating the deflections and bending moments, as well as the calculation error.

Example 2. Example 1, considered previously, represents a case close to the limit, since the stiffness at the left edge of the beam tends to zero. For a more objective assessment of the accuracy of the approaches under consideration, a similar series of calculations was performed for the case of a smooth change in stiffness. Figures 8 and 9 present the results for option 2 of the stiffness change, namely for s₀ = 0.6 and s₁ = 0.2. The graph of the stiffness function for this option is shown in Figure 5 b.

Example 3. The problem of bending of a clamped beam under a distributed load. The sketch is shown in Figure 4 b. The calculation results in Figures 10 and 11 were obtained for option 1 of the stiffness distribution, i.e., for s₀ = 0.001 and s₁ = 1.

Discussion. Based on the analysis of the results shown in Figures 6 and 7, the following conclusions were formulated. In cases of significant changes in stiffness, as in Example 1, only the exact stiffness matrix (27) provided obtaining results comparable with an analytical calculation. In this case, the threshold accuracy of 5% was attained with five finite elements for displacements (Fig. 6 b) and with six finite elements (Fig. 7 b) for bending moments. Using the approximate stiffness matrix (11), similar accuracy was attained only with 90 and 180 finite elements for deflections and moments, respectively. Discrete approximation allowed us to obtain the threshold accuracy for deflections with 220 finite elements. For moments, the required accuracy for discrete approximation was not attained even with 400 elements.

The calculation results for a distributed load, shown in Figures 10 and 11 for Example 3, are similar to those shown in Figures 6 and 7 for Example 1. The key difference is that the accuracy of the results when using the approximate matrix (11) and the discrete approximation has decreased. Thus, it can be concluded that the calculation results are affected not only by the stiffness matrix — exact or approximate — but also by the vector of nodal forces from distributed loads.

Based on the calculation results for a smooth change in stiffness, shown in Figures 8 and 9 for Example 2, the following conclusions were drawn. The approximate (11) and exact (27) finite element stiffness matrices for a smooth change in beam stiffness along its length allowed for comparable results with fewer elements. For example, to attain an error of 5% when calculating deflections using the proposed matrices, no more than four finite elements were required, compared with eight for discrete approximation. Thus, even in the case of a smooth change in stiffness, a linear approximation within a finite element allowed for a two-fold reduction in the degree of discretization of the calculation model, while maintaining comparable calculation accuracy.

The results presented in Figure 8 b and Figure 9 b have a unique feature. The results obtained by the approximate stiffness matrix (11) with the same number of elements showed a smaller deviation from the analytical solution than the results obtained by the exact matrix (27). This was due to the convexity (Fig. 5) of the stiffness change function along the length. In this case, the area of the figure bounded above by the approximation was smaller than the initial. Consequently, the deflection value obtained by the linear approximation was greater (0 а). From the results of calculating the coefficients of the approximate (11) and exact matrix (27), shown in Figure 3, it was clear that the approximate matrix (11) had greater stiffness, which partially compensated for this effect and reduced the total error. Below is matrix [ΔK], whose coefficients were calculated from the formula where где  — coefficients of matrix (11);  — coefficients of matrix (27). The calculations were performed for parameters ψ₀ = 0.664, ψ₁ = 0.2, D₀ = 1.0 and L = 0.6, which corresponds to the fifth finite element for option 2 of the stiffness distribution.

(38)

The average value of the coefficients of matrix (38) is 1.098. Therefore, the stiffness of finite element 5 will be approximately 10% higher if the approximate stiffness matrix is used (11).

Conclusion. This paper proposes and verifies an approach to constructing stiffness matrices for beam finite elements with linear stiffness change along their length. Both approximate (based on a variational solution to the problem) and exact (by analytical integration of a differential equation) stiffness matrices and corresponding nodal load vectors are obtained. The key advantages of the resulting matrices are the following.

1. Analytical form of approximate (11) and exact (27) stiffness matrices is presented, as well as approximate (15) and exact (29) vectors of nodal forces from distributed loads, which makes them easily implemented in existing FEM programs — both in the author's implementation in Python and in commercial packages.
2. The presented verification examples (Fig. 6–11) allow us to state that even with a relatively rough discretization (a small number of elements), the linear approximation describes adequately the stress-strain state of beams with both smooth and abrupt changes in stiffness (including modeling of zones of significant decrease in stiffness).

Linear stiffness approximation provides an optimal balance between accuracy, complexity, and computational resources, making the proposed method a practical tool for analyzing real-world structures with variable stiffness. A comparison of calculation results using an approximate and exact stiffness matrix shows that the approximate matrix allows for attaining the required accuracy with high precision and lower computational costs for beams with smoothly varying stiffness along their length. For rapid changes in stiffness along the length, as well as in extreme cases, to obtain significantly higher accuracy while maintaining the same number of elements, it is recommended to use the exact (27) stiffness matrix. Discrete approximation is only applicable in cases of smooth changes in stiffness along the length.

The resulting stiffness matrices open up possibilities for more accurate and efficient analysis of real structures, including reinforced concrete elements, taking into account physical nonlinearity — this determines the direction of further research. It is also planned to obtain a geometric stiffness matrix for solving stability problems and a mass matrix for dynamic analysis.