Extension of Castigliano’s method for isotropic beams

In the present contribution Castigliano’s theorem is extended to find more accurate results for the deflection curves of beam-type structures. The notion extension in the context of the second Castigliano’s theorem means that all stress components are included for the computation of the complementary strain energy, and not only the dominant axial stress and the shear stress. The derivation shows that the partial derivative of the complementary strain energy with respect to a scalar dummy parameter is equal to the displacement field multiplied by the normalized traction vector caused by the dummy load distribution. Knowing the Airy stress function of an isotropic beam as a function of the bending moment, the normal force, the shear force and the axial and vertical load distributions, higher-order formulae for the deflection curves and the cross section rotation are obtained. The analytical results for statically determinate and indeterminate beams for various load cases are validated by analytical and finite element results. Furthermore, the results of the extended Castigliano theory (ECT) are compared to Bernoulli–Euler and Timoshenko results, which are special cases of ECT, if only the energies caused by the bending moment and the shear force are considered. It is shown that lower-order terms for the vertical deflection exist that yield more accurate results than the Timoshenko theory. Additionally, it is shown that a distributed load is responsible for shrinking or elongation in the axial direction.


Introduction
Until the end of the nineteenth century, sophisticated graphical methods (e.g., by Christian Otto Mohr, Maurice Levy and Karl Culmann) were applied to analyze frame and truss structures, see Kurrer [1]. Then, around 1900, energy-based principles became more and more popular for modeling and dimensioning of structures. The driving force behind this rather new mathematical concept for analyzing constructions was the field of thermodynamics that had become an independent scientific discipline with a strong mathematical base some decades before. The introduction of thermodynamic concepts and quantities in structural design was a major challenge, which was not accepted in the scientific community for some decades, see the historical summary of Kurrer [2]. One of the energy theorems which was better recognized was Castigliano's method. The origin of Carlos Alberto Castigliano method dates back to 1879 and his diploma thesis [3]. The first of his theorem states that the partial derivative of the strain energy with respect to a certain displacement is equal to the force applied at this point. The second theorem of Castigliano states that the partial derivative of the complementary strain energy with respect to a certain force is equal to the displacement at the same location and in the direction of this force. Menabrea's minimum principle of the elastic energy is a special case of Castigliano's second theorem when the displacement vanishes. This allows for the calculation of redundant forces of statically indeterminate beams.
Within the linear regime (i.e., small strains and small deflections), Castigliano's second theorem provides a very efficient way for the calculation of mechanical quantities for engineering constructions. Instead of solving differential equations, one sets up the strain energy. Hence, it is still a very popular method in undergraduate classes because even without computational assistance, it is possible to find approximate solutions for frame structures like bridges by hand within an acceptable amount of time. Several practical examples are shown in the elementary textbooks from Ziegler [4], Parkus [5], Timoshenko and Young [6] and the contribution from Charlton [7]. The main advantage of Castigliano's method is that statically indeterminate beams (i.e., with multiple kinematical constraints) can be solved in a time-efficient manner. By considering the energy due to the shear stress, the shear effect on the structural deflection can be easily taken into account, see Timoshenko [6].
In particular for static deflections of frame constructions, it is highly advantageous to use energy methods if analytical solutions are required: Finding the bending moment distribution and, in addition, the shear and the normal force distribution, one can easily compute an analytical expression for the deflection by Castigliano's theorem. If the concept of dummy or phantom force is introduced, the deflection can be calculated at any desired location, see Dym [8]. The theorem is modified by Fu [9] who computes the deflection of a rectangular plate with two adjacent built-in edges and two adjacent edges free. In Ziegler and Irschik [10], it is shown that Maysel's formula can be understood as a generalization of Castigliano's method. A generalization of the reciprocity theorem is given in Irschik [11].
In this work, the knowledge of the Airy stress function is required for determining the strain energy. An elegant method is presented by Boley and Tolins [12], who apply the method of successive approximation for the bipotential equation: In an iterative manner, they find a more accurate solution for the stress function by each iteration. Irschik [13] extended this method by finding the stress function of the fourth iteration and thus finds a solution for a statically indeterminate (i.e., clamped-hinged) beam. Gahleitner and Schoeftner [14] modified the Boley-Tolins approximation method for anisotropic materials: the compatibility equation, which is a more complicated form of the bipotential equation in case of anisotropy, is iteratively solved for the Airy stress function. Analytical results perfectly agree with two-dimensional FE results in ABAQUS for a clamped-hinged beam. The iterative Boley procedure is applied to the charge equation of electrostatics in order to compute the electric potential as a function of the displacement field, see Krommer and Irschik [15]. In [16] Schoeftner and Benjeddou successfully applied Boley's iterative procedure in order to solve the compatibility equation and the charge equation of electrostatics simultaneously. For a mechanically and an electrically loaded piezoelectric bimorph, higher-order solutions for the deflection curves and for the electric potential are found.
This contribution is organized as follows: The Castigliano's theorem is extended in such a manner that a relation between the strain energy and the integral over stress-weighted deflection is derived. This is denoted as extended Castigliano's theorem or ECT. The goal of this contribution is to find more accurate formulae for the deflection curves of isotropic beams with rectangular cross section for which the computation of the complementary strain energy is required. Results from Boley and Tolins [12] for the Airy stress function of isotropic beams are regarded. Analytical formulae for the horizontal and the vertical deflection and the cross section rotation are obtained, which can be easily calculated if the bending moment, the shear force and the axial force distribution are known. It is shown that these results go beyond the Bernoulli-Euler (BE) and also the Timoshenko (TS) beam theory and properly take into account the thickness-to-length ratio as an important parameter. The ECT results are compared to analytical and also finite element results for various load cases and boundary conditions.

Extension of Castigliano's theorem (ECT)
The second theorem of Castigliano gives a relation between the complementary strain energy U * and the deflection v i in the direction of the external force P i . Mathematically speaking, the partial derivative of the complementary strain energy with respect to the force yields the deflection at the same location in the direction of this force, see, e.g., Ziegler [4] and Parkus [5], It is noted that this formulation is correct on beam level, but not in continuum mechanics when a single force might cause singularities. Hence, this formulation requires some adaptations. Assuming small deformations, so that the linearized theory of elasticity holds, the complementary strain energy U * , which is a function of stresses, can be replaced by the strain energy U . For beam-type structures, the usual way to calculate the strain energy is to find the bending and the shear force distribution as a function of the applied external load. Considering the energy due to bending only, Eq. (1) yields the same result as the Bernoulli-Euler theory. If the shear energy is also considered, one finds the same results as Timoshenko's theory which gives more accurate results for thicker beams, see Parkus [5]. In this section, the second theorem of Castigliano is extended in such a manner that it yields reliable results also for moderately and for thick beam structures, i.e., with a thickness-to-length λ < 1/2. Furthermore, beside the axial and the shear stress components σ x x and σ xy , the stress component σ yy (=stress in the thickness direction) is also considered in this contribution. Hence, the outcome in this contribution will be even more accurate than Timoshenko's result.
2.1 Derivation of the extended Castigliano theorem: ECT Figure 1 shows three configurations of a deformable body with thickness b: the reference or initial configuration, when the body is unloaded (left), the intermediate configuration after the first load increase, and the current or final configuration after the second load increase (right). The total load for load case I and load case II are identical, and only the order of the load application is reversed.

Load case I
For load case I the body is subjected to the external load P i t * i first. Here P i is a scalar parameter, which is quasi-statically increased from 0 to the end value P i , and t * i is the normalized traction. The dimensions of P i and t * i are N and m −2 . Integrating the product with respect to a certain area dS = bds, one observes that P i t * i is an external traction with dimension Nm −2 . For the second load step, P i is increased from P i to P i + ΔP i . The total energy stored in the elastic body over the first and the second load steps (denoted by U (P i ) and ΔU (P i )) must be equal to the total work of the external force W ext + ΔW ext : Fig. 1 Load case I (top): first, the undeformed elastic body (left) is subjected to the traction load P i t * i , and then the load is increased to (P i + ΔP i )t * i . Load case I I (bottom): now the loads are applied in reverse order. First the incremental load is applied ΔP i t * i , and then, the load is increased from In Eq. (2), u is the displacement field in the intermediate configuration (after the first load step), and Δu + u is the total displacement field in the current configuration (after the second load step).

Load case I I
For load case II, the loading order is reversed: first, the incremental load ΔP i is applied which causes the displacement field Δu and the strain energy ΔU + . Then, the load is increased from ΔP i to ΔP i + P i . The additional displacement for the second load step is u (which equals the deflection of the first load step for load case I ), and the additional strain energy is U + . Again the total energy stored is equal to the total external work performed: Combining Eqs. (2) and (3), it is clear that the total work performed and the final strain energy must match. Considering that holds for an incremental load increase ΔP i → ∂ P i , one finds Equation (5) states that the partial derivative of the strain energy with respect to P i equals the weighted displacement over the traction surface ∂ B i . For practical applications, Eq. (5) is not suitable because the deflection cannot be evaluated outside of ∂ B i . To circumvent this severe restriction, one introduces a so-called dummy force. Hence, replacing P i t * i by P r t * r + P d t * d over the surfaces ∂ B r and ∂ B d , where the subscript r and d denote real and dummy load tractions, and assuming an infinitesimal load step ∂ P d only for the dummy load (but not for the real load) and considering P d /P r → 0 (i.e., the dummy load is negligible), one finds the more applicable formula This means that the weighted displacement over the area of the dummy traction ∂ B d can be calculated by the partial derivative of the strain energy with respect to P d . It is noted that the strain energy U (P r , P d ), that is a function of the real and of the dummy load, must be computed first before differentiation is performed: One must not switch the mathematical operations (i.e., limit operation and partial derivative) in Eq. (6).  Figure 2 shows a beam which is loaded by the dummy forces P dx , P dy (in the x and y-directions) and the dummy moment P dz at x = x P in the z-direction. The reaction forces and the moment at the vertical edges Hence, the stress resultants N d , Q d = M d,x and M d can be interpreted as well as the loads P dx , P dy and P dz at ∂ B P as external forces (Fig. 2a). In our case the traction vectors over the lower and upper surfaces t d low = t d upp = 0 (at t d low = t d upp = 0) are zero and hence, the corresponding surfaces can be disregarded.
For the loads at the left and at the right boundary and at x = x P , the following conditions hold between the stress resultants and the stress components: From a didactical point of view, the new expression ∂ B d u · t d ds is considered instead of the weighted displacement ∂ B d u · t * d dS in Eq. (6). It holds t d = P di t * d . For a beam with rectangular cross section and constant width b (hence the differential is ds = dS/b), Eq. (6) can be rewritten as Approximating the displacement components by the Taylor series at y = 0, the first integral at the right-hand side in Eq. (11) yields with the thickness-to- Equation (13) is further simplified by neglecting λ-dependent terms: At and for the dummy load triple (P dx , P dy , If classical boundary conditions are assumed, where either the kinematical or the dynamical boundary conditions vanish at x = 0 and L, Eqs. (14) and (15) vanish. Evaluating the original integral of Eq. (16) and considering that only one of the dummy forces is active (while the others are zero) and that t d P = P di t * Equation (17) shows that the dummy load method yields the deflection or rotation at location x P .
Iteration Airy stress (normal force)

Complementary strain energy as function of the bending moment and normal force
The strain energy for a two-dimensional elastic system under plane-stress assumptions reads This relation can be expressed in terms of the stresses if an isotropic material law is used. The stresses are derived from the Airy stress function ψ(x, y): Substituting Eq. (19) into Eq. (18) and considering isotropic material, one finds that the strain energy only depends on the Airy stress function. For a rectangular isotropic beam, Boley and Tolins derive an approximation of the stress function by an iterative routine, see [12]: The Airy stress function is a function of the bending moment and the normal force and their even-numbered spatial derivatives. The main results are summarized in Table 1, where ψ M (x, y) and ψ N (x, y) denote the stress functions which depend: -on the bending moment M(x) (caused by a distributed load q(x) at the upper surface at y = c, see case I in [12]) and -on the normal force distribution N (x) (caused by the surface traction n(x)/2 in axial direction at the upper and lower surface y = ±c, see case I I I in [12]).
The subscript 1 in ψ M 1 (x, y) means that the stress function is directly related to the bending moment, whereas the subscripts 2, 3, , y), . . . denote the higher-order stress functions which are proportional to the second, the fourth or the sixth derivative of the bending moment. Table 2 also shows the well-known linear relation between the bending moment and the axial stress on beam level ψ M 1 yy = σ M 1 x x = M(x)y/I (2 nd column, 2 nd row). It can be shown that the higher-order terms are so-called self-equilibrated stresses which do not produce any moment, shear force or axial force, i.e., In an analogous manner, this also holds for the Airy stress ψ N (x, y) generated by axial loads and an external normal force. The total stress function is with Inserting Eqs. (19), (20) and (21) and Table 1 into Eq. (18), one finds the strain energy as a function of the thickness-to-length ratio λ = 2c/L in the following form: with In Eq. (22) the term O(λ 5 ) is the truncation error, i.e., the error of the strain energy is proportional to λ 5 . If one also considers the fourth iteration of Boley's results (not shown in Table 1

Partial derivative of the strain energy with respect to the dummy load
The bending moment and the normal force depend on P dx , P dy or P dz and are split into portions M d (x, P dy , P dz ) and N d (x, P dx ) which are generated by the dummy loads only, and portions M r (x), N r (x) which are caused by the real load only: The equilibrium conditions on beam level relate the vertical load q(x) and the horizontal load n(x) to the bending moment and the normal force, see Ziegler [4] and Parkus [5]: Evaluating the left-hand side of Eq. (6) and considering Eq. (17), one finds for the horizontal deflection curve (P dy = P dz = 0) and for the vertical deflection curve For the rotation one also finds (P dx = P dy = 0) It is noted that the partial derivatives ∂ M d /∂ P di and ∂ N d /∂ P di in Eqs. (26)-(28) are replaced by the quotient M d /P di and N d /P di (i = x, y, z): One observes that these relations are independent of P di . Equations (26)-(28) also show two surprising results for the deflections which are briefly discussed in the following: -horizontal deflection u 0 (x P ): The first term in Eq. (26) shows the fundamental term from the normal force distribution N r (x) that is also known from the elementary theory. The second line shows that the vertical load q(x) causes a horizontal deflection of the midline, and its effect is proportional to the product of thickness ratio and Poisson ratio λν. Higher-order derivatives of the normal force (i.e., n ,x = −N r,x x = 0) have a lower-order influence (∝ λ 2 ). -vertical deflection v 0 (x P ): The first and the second terms in Eq. (27) show the influences of the bending moment (Bernoulli-Euler, ∝ λ 0 ) and of the shear force (Timoshenko correction, ∝ λ 2 ). It is interesting to note that the load q(x) also causes an additional vertical deflection due to the Poisson effect ν = 0 which has the same order of magnitude like Timoshenko's shear correction. It can be shown that this is a consequence from V σ yy ε yy dV that is typically neglected if the principle of virtual work is applied if the variation of the thickness deformation vanishes (i.e., δε yy = 0 → V σ yy δε yy dV = 0). The last line shows the λ 4 -dependency from the higher load derivations q ,x (x) and q ,x x (x), but from a practical point of view they are negligible (see Sect. 3).

Verification of ECT with finite element or analytical results
For the verification of the presented theory, a cantilever (Sect. 3.1) and a simply supported beam (Sect. 3 Table 2 as a function of the non-dimensional length ξ = x/L for various load cases (either M 0 or F 0 acting at ξ = 1 or the load distributions q(ξ ) = q 0 and q(ξ ) = q 1 (1 − ξ)). These results are evaluated from Eqs. (26)-(28), which also show the contribution from each of the energy terms (U B , U S , U q , U coupl , U λ 4 ).

Deflections for a constant load q(x) = q 0
Figures 3a-c show the normalized vertical deflectionv = v 0 /v BE , the axial deflection u 0 at the free end (ξ = 1) and the normalized cross section rotationφ = ϕ/ϕ BE at ξ = 0.5. The analytical results for the tip deflections (Bernoulli-Euler (BE), Timoshenko (TS), ECT (up to order λ 2 ), see also Table 2 for load case q 0 ) and according to the 2D-results under plane stress assumptions (Eqs. (52)-(54)) read The results for the normalized cross section rotation at ξ = 0.5 read and the axial deflection is As expected, one observes from Fig. 3a that all theories that include the thickness-to-length ratio deviate from the Bernoulli (BE) theory: In particular for thick beams (i.e., 1/λ < 10), this influence becomes relevant. The zoom figure clarifies that the Timoshenko (TS) theory (blue) does not converge to the analytical solution (grey), while the results from the ECT do (red, magenta) because the energy U q also contributes. For a constant load distribution, the fourth-order effect does not occur ∂U λ 4 /∂ P dy = 0. One observes from Eq. (29) that the analytical 2D-solution contains a λ 4 -term that has practically no relevance; hence, ECT and analytical results agree within λ 2 -accuracy. The axial tip deflection is shown in Fig. 3b. The elementary BE and TS-theories do not predict any deflections, but the partial derivative ∂U q /∂ P dy due to the coupling energy does (u 0 ≈ −7.143 · 10 −13 m), see Eq. (31). It is noted that the analytical 2D-result equals the ECT result.
The cross section rotation is identical for BE and TS theories, see Eq. (31). Considering also the portion due to the normal stress σ yy , one finds that the ECT results (red, magenta) agree quite well with the analytical result. Figure 4 shows the results when a linear decreasing load is assumed. The analytical results are given in Table 3. The results for the deflection and the rotation with the ECT are very accurate and match as good as the previous results (see Fig. 4), but here the term ∂U λ 4 /∂ P dy , which is a λ 4 -term, also contributes to the vertical deflection, but this small difference is negligible and cannot be resolved even by the zoom in Fig. 4a: Hence, for an acceptable agreement with analytical results, it is enough to consider terms up to order λ 2 (i.e., U B , U S and U q ). For the axial deflection (Fig. 4b) and the rotation (Fig. 4c), the ECT results match best with the analytical 2D-results: u 0 (ξ = 1) = −q 1 L 2 λν/4E A ≈ −3.572 · 10 −13 m, see Table 2 and Eq. (59).

Numerical example: simply supported beam
Next a simply supported beam is considered. Two load cases are considered: a linear increasing load distribution q(x) = q 1 x/L and a parabolic load q(x) = q 2 x 2 /L 2 . The results of the ECT (see Table 3) are compared to finite element (FE) results. A Q8-element (eight-node quadrilateral element with quadratic ansatzfunctions for the displacement) is used for the finite element mesh. The total number of elements is 2000 for the moderately thick beam (i.e., L = 0.5 m, c = 0.05 m → λ = 1/5, 100 × 20 elements in the x and y direction (Fig. 5)).  Figure 5a shows the displacement field from the FE calculation for a moderately thick beam (λ = 1/5) on which the load q(x) = q 1 x/L acts. The deflection curves v 0 (ξ ) and u 0 (ξ ) are plotted in Figs. 5b and c. In the middle of the beam, the vertical deflection is v F E (ξ = 0.5, η = 0) ≈ 2.54 × 10 −11 m, which is in very good agreement with ECT results (red, magenta, Fig. 5b): v 0 (ξ = 0.5) ≈ 2.53 × 10 −11 m. In comparison, the Timoshenko theory is also close to the FE solution v T S (ξ = 0.5) ≈ 2.56 × 10 −11 m, unlike the BE theory which underestimates the deflection as expected v B E (ξ = 0.5) ≈ 2.33 × 10 −11 m. The horizontal deflection curve of the beam axis (Fig. 5c) vanishes for the elementary BE and Timoshenko theories (black and blue). The ECT and the FE results show a parabolic horizontal deflection of the centerline

Deflections for a linear load q(x) = q 1 x/L
see also Table 3 (4 th row, 3 rd column), which is due to U coupl (energy coupling of distributed load and normal force).

Deflections for a parabolic load q(x)
Finally, the influence of the thickness ratio for the simply supported beam for q(x) = q 2 x 2 /L 2 is studied. 200 finite elements are used in the x-direction (i.e., Δx = 0.0025). The element shape is (almost) a square element with Δy ≈ Δx, depending on the thickness ratio, Δy might be adjusted (e.g., λ = 1/3: 200 × 66 = 13200 elements). Figure 6a shows the normalized vertical deflectionv = v 0 /v BE at ξ = 0.5. For λ = 1/5 the FE result (grey-triangle) is 8.56 % larger than the BE result. Both ECT results are very close to the FE-solution (+8.69 % (red, λ 2 -accurate results) and +8.68 % (magenta, λ 4 -accurate results including U λ 4 ). The Timoshenko result (+9.82 %, blue) overestimates these results. Nonzero horizontal deflections u 0 (ξ = 0.5) are only predicted by the ECT results that contain the coupling term U coupl (Fig. 6b), i.e., the horizontal deflection is constant, u 0 ECT (ξ ) ≈ −1.49 · 10 −14 , and does not depend on the thickness ratio. For very thick beams, one observes that the FE results underestimate these results. Without going into detail, the consideration of the higher-order self-equilibrated stresses (see Table 1 in Boley-Tolins [12]) that become more relevant for thick beams, and a finer mesh size may increase the agreement between numerical and analytical results. From a practical point of view, such thick structures are often not considered as beams. Nevertheless, the range validity of the depends on the quantity to be investigated: For the vertical deflection, the accuracy of the ECT and Timoshenko results is very high even for thick beams. But the Timoshenko theory always yields a vanishing horizontal deflection of the x-axis, while the ECT results yield a reliable output as long as λ < 1/4 holds.
3.3 Numerical example: clamped-hinged beam with q(x) = q 0 Beside Sects. 3.1 and 3.2, where statically determinate beams are investigated, also a clamped-hinged beam is considered here. For the solution of the statically indeterminate beam, the problem is divided into a cantilever with constant load distribution q(x) = q 0 and into a cantilever with a yet unknown single force F 0 . The condition v(L) = 0 determines the redundant force F 0 (see Table 2 for ECT and sections A.2 and A.3 for the plane stress results). Figure 7a shows the normalized vertical deflection in the middle as a function of the thickness-to-length ratio and the deflection curve for a thick beam (Fig. 7b) of a clamped-hinged beam. All results except the BE results are in good agreement. One observes that the influence of the compressive stress σ yy (which is not covered by the Timoshenko result) additionally allows the extended Castigliano results to converge to the analytical plane stress solution (for λ = 1/5: error is 0.2 % (magenta, red) compared to the 0.8 % error by Timoshenko (blue), see zoom in Figs. 7a and b). As expected, the BE theory yields a much stiffer configuration beam (error ≈ 28 %).  Table 2 for ECT and section A for the plane stress results). Figure 8 shows the normalized vertical deflection at ξ = 0.5 as a function of the thickness ratio and the deflection curve for λ = 1 5 . In comparison with the statically determinate results (e.g., Figs. 3a or 4a), where the relative deviation of the vertical deflection is only 1 % for λ = 1/10, the BE theory underestimates the results (12 % error) here. The influence of shear is much higher than for the other configurations, whereas the deformation by the compressive stress σ yy is rather low (i.e., Timoshenko results (blue) are close ECT results (red, magenta) and to the analytical result (grey)).

Conclusion
In this present contribution, an extended Castigliano theorem (ECT) is developed in order to find more accurate results for the deflection curves of beam-type structures, i.e., for the horizontal and the vertical deflection and the cross section rotation. Introducing a so-called dummy load, which acts in addition to a real load that deforms an elastic body, one finds that the partial derivative of the (complementary) strain energy with respect to the dummy load is equal to the weighted deflection over the surface where the dummy load acts. This is a generalization of Castigliano's second theorem. For isotropic beam-type structures, the strain energy can be expressed by the Airy stress function for which accurate relations to normal force, shear force, bending moment and the load distributions exist. The analytical outcome shows that the results for the deflection curves with ECT go beyond the Timoshenko theory because the normal stress in thickness direction is properly taken into account. Statically determinate (a cantilever and a simply supported beam) and statically indeterminate beams (clamped-hinged and clamped-clamped beam are considered for the verification of the ECT results that are compared to finite element and analytical 2D-results (plane stress). Special cases of the ECT formulation are the results from Bernoulli-Euler and the Timoshenko, if only the bending and the shear energy are regarded. Particularly for statically determinate beams, the obtained ECT results are even more reliable than the Timoshenko theory.
σ xy (x, y) = 0, (36) σ yy (x, y) = 0. (37) The horizontal deflection, the vertical deflection and the rotation angle are A.2: Solutions resulting from the force F 0 at x = L The Airy stress function and the stresses read The horizontal deflection, the vertical deflection and the rotation angle are A.3: Solutions resulting from the constant distributed load q(x) = q 0 The Airy stress function and the stresses read The horizontal deflection, the vertical deflection and the rotation angle are A.4: Solutions resulting from the linear decreasing load q( The Airy stress function and the stresses read The horizontal deflection, the vertical deflection and the rotation angle are