In this work1, we follow an optimization strategy for beam systems using heuristic techniques such as Quantum Annealing (QA) and Simulated Annealing (SA). The optimization is carried out in two stages: first, minimizing the Hamiltonian with respect to beam displacements to determine equilibrium configurations under external loads; second, iteratively optimizing the cross-sectional areas of the beam members to maximize system stifness while maintaining volume constraints. The approach combines modern quantum-inspired techniques with classical computational methods to eficiently explore the solution space. Preliminary numerical results are given at the end.
Contents 2. Combinatorial Optimization, QUBO and Adiabatic Quantum Optimization 3 3. Framing the problem 4 3.1. FEM Approximation: Shape Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3.2. Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4. Truss systems: No Bending 8 4.1. Mechanical Equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4.2. Optimization of Cross-Sections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4.3. Two-stage optimization strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4.4. Quantum annealing implementation (D-Wave) . . . . . . . . . . . . . . . . . . . . . . . 11
©2025Copyrightforthispaperbyitsauthors. UsepermittedunderCreativeCommonsLicenseAttribution4.0International(CCBY4.0).
A. A Coercivity Criterion for a Rank-One Perturbed Quadratic on the Nonnegative Orthant 13
The optimization of structural systems, particularly beam systems, is a longstanding and critical problem in engineering design. The goal is to enhance performance metrics such as strength, stifness, or material eficiency while minimizing resource use. Beam systems—composed of straight members connected at joints—are fundamental in applications ranging from bridges to spacecraft and micro-lattice structures.
Traditional simulation tools such as the Finite Element Method (FEM) provide accurate predictions of displacements and internal forces, but they become computationally demanding as the number of elements, load cases, and design variables grows. Likewise, gradient-based optimization techniques (e.g., Sequential Linear Programming, Sequential Quadratic Programming, or interior-point methods) can be very eficient near a local optimum yet often struggle with highly non-convex, discrete, or multi-modal design spaces typical of topology and sizing optimization. Constraint handling, sensitivity discontinuities, and the curse of dimensionality further complicate their use.
Heuristic and stochastic search methods—including Genetic Algorithms, Particle Swarm Optimization, Simulated Annealing (SA), and more recently quantum-inspired and quantum metaheuristics—have therefore gained traction. These methods are relatively insensitive to non-convexity and can escape shallow local minima, making them attractive for high-dimensional combinatorial formulations. Among them, Quantum Annealing (QA) stands out as a hardware-accelerated heuristic tailored to problems expressible as Quadratic Unconstrained Binary Optimization (QUBO) or, equivalently, Ising models. QA systems (e.g., D-Wave processors) implement an analog realization of adiabatic quantum computation, promising advantages in exploring rugged energy landscapes by exploiting quantum tunneling and massive parallelism at the hardware level. Our work builds on the recent paper [1] where optimization of truss systems with axial strength only was studied. 1.1. Why QUBO? Many engineering design problems ultimately boil down to choosing one option among a finite set for each of many coupled decisions: select or discard a member, assign a cross-sectional size from a catalog, toggle load paths, etc. Encoding such discrete choices as binary variables is natural. The challenge is to express the objective and constraints using only linear and quadratic terms in those binaries—which is precisely what a QUBO model requires. Formally, a QUBO minimizes
min
x∈{0,1}
x x + c x + 0,
(
A QUBO representation ofers several practical benefits beyond hardware compatibility: • Unified encoding of objectives and constraints, allowing heterogeneous physical requirements to be aggregated into a single cost function. • Sparsity exploitation: FEM-derived coupling is often banded or sparse, leading to structured matrices that can be embedded more eficiently on restricted hardware graphs. • Hybrid decomposition: QUBOs are amenable to classical pre- and post-processing (e.g., variable ifxing, roof duality bounds) and to hybrid algorithms that split the problem between classical CPUs/GPUs and quantum co-processors. 1.2. Importance of Quantum Computing for Engineering Optimization Quantum computing, broadly construed, is poised to augment classical computing in domains where combinatorial explosion or rugged objective landscapes limit classical heuristics. For structural optimization, the potential advantages include: • Scalability for discrete design spaces: Many sizing/topology problems are NP-hard. Quantum hardware that directly samples low-energy states could provide speedups or, at minimum, highquality candidate sets faster. • Massive parallel sampling: QA can return thousands of low-energy samples in a single run, enabling statistical analysis of design alternatives and robust optimization under uncertainty. • Hybrid quantum–classical workflows: Even without quantum advantage, integrating QA as a subroutine (e.g., for dificult subproblems or warm starts) can reduce wall-clock time and improve solution diversity. • Cross-disciplinary transfer: Techniques developed for logistics, finance, and machine learning
QUBOs readily translate to structural design, accelerating methodological innovation.
Nonetheless, current devices are noisy, have limited qubit counts and connectivity, and require careful formulation to realize any benefit. Therefore, problem reduction, parameter tuning, and classical post-processing remain integral.
In this work, we follow the procedure outlined in [1], where a two-step optimization routine was
outlined: (
The rest of the paper is organized as follows. Section 2 discusses the basics of QUBO. Section 3 formulates the structural optimization problem and derives its Hamiltonian representation, including constraint penalties. Section 4 gives the first results for truss systems where bending is not contemplated. Section moves on to presenting numerical tests on the Euler-Bernoulli Theory, where truss bending is taken into account.Section 5 concludes with a summary of future research avenues, including extensions to dynamic loading and reliability-based design. 2. Combinatorial Optimization, QUBO and Adiabatic Quantum
Optimization A wide class of combinatorial optimization problems can be approached with quantum algorithms. We briefly recall the standard formulation of such problems. Let : → R be a real-valued cost (or objective) function defined over a set of decision variables . The task is to find an element ⋆ ∈ that minimizes : ⋆ = arg min () ∈ subject to a collection of constraints ℓ() = 0 (ℓ = 1, . . . , ),
ℎ() ≤ 0 ( = 1, . . . , ).
A convenient strategy is to turn the constrained problem (
() = () + ∑︁ ℓ ℓ()2 + ∑︁ max[︀ ℎ(), 0]︀ .
ℓ Here ℓ > 0 and > 0 are penalty coeficients chosen large enough that any violation of the constraints is energetically disfavored.
When the decision space is binary, = B with B = {0, 1}, and the penalized energy is quadratic in the variables, the problem is known as a Quadratic Unconstrained Binary Optimization (QUBO).1
QUBO instances are directly related—indeed, equivalent—to Ising models, where binary spins ∈ {−1, 1} replace {0, 1} variables. A simple change of variables maps one formulation to the other; the explicit transformation is presented in the next section.
Consequently, many combinatorial problems can be reframed as finding the ground state of a quantum Hamiltonian, enabling the use of physics-based techniques. One particularly promising framework is Adiabatic Quantum Computation (AQC), which relies on the adiabatic theorem: if a system is evolved slowly enough, it remains in its instantaneous ground state.
In the literature, AQC is often used interchangeably with Quantum Annealing (QA). Strictly speaking, QA encompasses protocols that need not be fully adiabatic, but in this work we use the two terms synonymously.
We will not discuss Adiabatic Quantum Optimization any further in this contribution. For comprehensive reviews on Ising Models and Quantum Adiabatic Optimization, see, for example, [2, 3], and [4, 5] for the quantum adiabatic theorem. See also the original paper [6].
We now move on to framing the problem and seeing how it can be converted to QUBO.
Consider a beam element belonging to a 2D frame system in the Euler-Bernoulli theory of elasticity. The goal is to start from the continuous formulation (internal energy and external work) and derive the complete FEM discrete formulation, including axial and bending stifness.
The total potential energy of such a beam can be defined as:
Π = − • is the internal strain energy, • is the work of external forces. 1Higher-order terms can be reduced to quadratic form with ancillary (auxiliary) variables. For instance, a cubic term 123 can be replaced by 14 by defining 4 = 23 and enforcing the relation through a penalty such as 34 +23 −2 24 −2 34, which vanishes only when 4 = 23.
Here, the work done by external forces is taken positive. Thus, the work done on the system is negative. That explains the minus sign in the -term. For an element of length , the axial strain and the bending curvature are given by • axial strain: = • bending curvature: = ,
, where ′ is the axial displacement in the beam system, whereas is the transverse displacement. Here, we assumed that the beam axis is parallel to the -axis, while describes the vertical displacement in the beam system, i.e. along the -direction. Later, we will have to consider also the global frame − . We will denote deformations in the and directions and , respectively. Switching from a local to the global frame is rather simple, and will be shown in the following.
Thus, the internal energy is given by the sum of the strain energy and the bending energy: = • ′(): axial displacement along the local axis of the beam, • (): transverse displacement, • : Young’s modulus, • : cross-sectional area, • : second moment of area.
Notice that the second moments of area depend on the beam sections. They are indeed defined as Here, we take the beam cross sections to be circular. Other geometries can be contemplated. In the case of circular sections, the second moment of area along z reduces to the simple formula ∫︁
= 2 , =
2 . ∫︁
∫︁
= 2 = 4 4 = 2 4
. = ∑︁ ,
where the ’s are the global vertical displacements. 3.1. FEM Approximation: Shape Functions An analog formula holds for .
The work done by external forces in the form of external loads applied at the joints can be simply written as One of the goals of this paper is to solve the beam system, i.e. determine strains and stresses. We do so by treating each beam as a discrete entity, employing techniques from FEM (Finite Element Methods).
We think of the as a discrete element with two nodes (node 1 and node 2). Then, the local degrees of freedom for a 2D beam are:
dlocal = [′1, 1, 1, ′2, 2, 2]
• ′1, ′2: axial displacements at the nodes,
(
In the usual Euler -Bernoulli theory, the ’s are the derivatives of the vertical displacement: () = () . (10) In the following, as is usual within the FEM, we will treat them as separate variables. How to reconcile = / will be clear in a moment.
In the following, we will be working in the linear regime. Thus, we assume linear interpolation for the axial displacement: ′() = 1()′1 + 2()′2, 1 = 1 − , 2 = , = / Thus, we have for the axial strain:
For the vertical displacement, it is customary to use Hermite shape functions (ensuring 1 continuity) as interpolators:
() = 1 ()1 + 2 () 1 + 3 ()2 + 4 () 2 with: 2 () = ( − 2 3 () = 3 2 − 2 4 () = (− 2 + 2 3, 2 + 3), 3, 2 + 3), 2() = ∑4︁ 2 , 2 2
=1 where Klocal is given in the local system by the simple formula Let us now move on to the bending energy.
= u Klocalu, Klocal = [︂ 1 −1 where is the dimensionless variable: = .
Therefore, the curvature takes on the following form: with d = [1, 1, 2, 2] . Being the shape functions polynomial, the derivatives are easily computed.
The discrete axial energy is then easily evaluated: =
Such a form is typical of elastic systems: this is just an upshot of the linear approximation. It is often useful to write employing a matrix formulation. Defining u = [1, 2], we can set which, once again, can be conveniently rewritten in matrix notation. Using the Hermite shape functions (12) we get: with: =
It turns out to be convenient to write the beam axial and bending energy in terms of the complete local stifness matrix:
The discrete bending energy is made of a curvature term: (14) (15) system.
(16) ⎡
⎢ 0 Klocal = ⎣⎢⎢⎢⎢⎢− 00 0 −
0 12 3 6 2 0 12 3 6 2 0 6 2 4 0 6 − 2 2 − 0 0 0 0 − 0 12 3 6 − 2
0 12 3 6 − 2 0 ⎤ 6 2 ⎥ 2 ⎥ 0 ⎥⎥⎥ − 62 ⎥⎦ 4 dlocal = [′1, 1, 1, ′2, 2, 2] where the degrees of freedom are, once again,
At this point, it is useful to quote how we can go from a local beam system to the global − If the beam forms an angle with the global X-axis, we define the matrix
T = ⎢⎢⎢⎢ 00 ⎡ 0 0 0 0⎤ ⎢− 0 0 0 0 ⎥ 0 0 00⎥⎥⎥⎥ , 0 1 0 0 ⎢⎣ 0 0 0 − 0 ⎥⎦ 0 0 0 0 0 1 = cos , = sin The global stifness matrix is the obtained via:
Kglobal = T K localT This way we can express the beam energy in the global system, which is more suitable for what we will be doing in the following.
Assembling all elements, we find the rather compact form: Π(d) = 1 d Kd − d F
2 • K = stifness matrix assembled from all transformed elements, • F = vector of nodal loads (e.g., vertical weight at a node).
The equilibrium condition is, of course, obtained by minimizing through minimization: However, solving such an equation might be a challenging task, especially for convoluted beam systems. This is where the QUBO optimization kicks in.
The optimization of a beam system is governed by its mechanical equilibrium under applied loads and the constraints imposed by material usage. The problem can be expressed mathematically as given below.
The optimization problem is mapped to a QUBO model for compatibility with quantum annealing solvers. To reach mechanical equilibrium, the Hamiltonian to be minimized is = Π(), and can be split into two terms, one for the internal system energy and another one for the external work: (17) (18) (19) (20) (21) where 3 enforces the volume constraint as a penalty term: ⎛
⎞2 3 = ⎝
∑︁ − TOT + slack⎠ .
(,)∈
This is the second step of the optimziation process (Step 2). Here slack refers to the slack variables to enforce the inequality constraints. Should they not be introduced, we would be enforcing an equality constraint, which in the present formulation is just as good.
Step 1 and 2 are repeated in a loop until we do not appreciate a significant change of the cross sections anymore. In that case, convergence has been reached and the configuration can be thought of as definitive for our optimization task. In the end, we should see a deformed beam structure (strains and stresses) and beams with diferent sections, so to optimize material usage.
This is the first step of the optimization process (Step 1). Explict formulas for 1 are given below, according to the theory we are considering (trusses, Euler-Bernoulli or Timoshenko). 2 will be of the form Minimizing leads to equilibrium configuration for the beam system, i.e. leads to the system displacements d.
Once the d is found, we can proceed to find the optimal cross-sections distribution for the beam. To do so, we aim at maximizing the global stif amtrix. Thus, the hamiltonian this time reads: = 1 + 2.
2 = −d F.
= − 1 + 3, 4. Truss systems: No Bending
The equilibrium configuration of a beam system minimizes the total potential energy, which consists of the elastic strain energy stored in the members and the work done by external loads. To get started, we begin with a system where bending is negligible, i.e. () ≈ 0 in the local system. Thus, the only variables are the ′’s in the local system which are mapped to the global ’s and ’s. The Hamiltonian for this system is given by: • 1: Elastic strain energy of the beam members.
• 2: Work done by external forces.
The elastic strain energy 1 is expressed as: • : Young’s modulus of the material, • : Cross-sectional area of member (, ), • : Length of member (, ), • (, ): Displacement components of node .
The external work 2 is: 2 = ∑︁ P · u ,
∈ where P represents the external force vector at the node . 4.2. Optimization of Cross-Sections (22) (23) (24) (25) (26) (27) (28) Once the equilibrium position is found, i.e. we determined the displacement vectors ⃗, we can attempt to maximize stifness while satisfying volume constraints. Here, we follow closely reference [ 1] with minor diferences. The optimization is formulated as follows:
Maximize: Subject to:
∑︁ (,)∈
2 , ∑︁ ≤ TOT, (,)∈ where is the strain in member (, ) and TOT is the total allowable volume of material.
The cross section is to be varied only by a small amount (say no more than 20 %), or otherwise we would be to far afield wrt the equilibrium configuration found in the previous step. The process repeats till convergence is reached. In that case, we have established the equilibrium position of the beam system with optimal cross-sectional distribution.
To set the problem properly, we discretize the areas in the following fashion: each truss member is assigned an -bit binary vector {,}=−10 , which is decoded via random per-bit coeficients , into an additive area update. Specifically, we compute
−1 = ∑︁ , ,,
=0 Δ = (︀ 2 ˜ − 1 )︀ max, ˜ =
∑︀=−10 ,
(0) + Δ , = ∈ [0, 1], with max not exceeding 20% its original value. The new areas { } are obtained through annealing, by minimizing global system keeping the volume not exceeding a given value:
= − 1 + 3 with 1 as before, while 3 given by equation (20). This configuration—absolute step updates, randomized bit weights, and hard volume filtering—ensures both rich exploration of the design space and feasibility with respect to the volume constraint.
We discuss now the case of the truss system in Figure 1. The truss system has 12 nodes and 29 trusses, all with equal sections.
We then apply an external weight at the rightmost node of the bottom row. 4.3. Two-stage optimization strategy The optimization proceeded in two tightly coupled stages. In the first stage, for any prescribed set of external loads, we minimized the total Hamiltonian of the beam system with respect to the nodal displacement field. This energy-based formulation (strain energy + potential of external forces) yields the equilibrium configuration directly, and can be achieved with standard FEM methods or simulated annealing. We used both, reaching similar results in similar times. A standard Python package for FEM is PyNite and its FEModel3D module.
The resulting displacement vector served as the state variable that feeds the second stage. In that second stage, we iteratively updated the cross-sectional areas of the individual beam members to maximize the global structural stifness (equivalently, to minimize compliance) while enforcing a fixed total volume (mass) constraint. The sizing problem was cast as a constrained optimization.
Figure 2 shows the undeformed and deformed configurations of the truss. Member thickness is proportional to the optimized cross-sectional area, so thicker lines indicate elements that the sizing algorithm kept or enlarged because they carry higher axial forces (or are needed to satisfy stifness constraints). Conversely, thin members correspond to low-force paths and could be candidates for further reduction or even removal if allowed by the design constraints.
The gray drawing shows the original, undeformed geometry; all members are plotted with the same topology but already scaled in thickness by their final areas so you can read the importance of each bar directly on the baseline configuration. The red drawing is the deformed shape, magnified ( ×7.8 ) to make displacements visible. Blue arrows point from the initial to the magnified positions and illustrate the displacement vectors. The vertical load is applied at the bottom-right node, and the largest rotations and translations occur close to that point, as expected.
Reading the plot: • Load path: The thickest bars cluster along the lower chord and the diagonals connecting the loaded node to the supports, revealing the primary force flow. • Eficiency after optimization: Members that are thin across both configurations carry comparatively little force; the optimizer reduced their area to save material while keeping global stifness and strength within limits. • Deformation pattern: The amplified red outline highlights global sway and local joint rotations; because the scale factor is stated, the true displacements can be recovered if needed. In summary, line thickness encodes the optimized section magnitude, while the overlaid, scaled deformation shows how the structure responds to the applied load. Together they provide an immediate visual link between structural demand (forces) and structural response (displacements).
Figure 2 is consistent with results found in [1] and is what we see at step 30 of the iteration. However, going further down the double-step procedure we do not seem to see any rapid convergence, with a changing topology as well. We leave further exploration about convergence for the future. 4.4. Quantum annealing implementation (D-Wave) To eficiently explore the high-dimensional design space, the second-stage search was initialized and periodically re-seeded using samples from a D-Wave quantum annealer. The discrete sizing/topology problem was encoded in Quadratic Unconstrained Binary Optimization (QUBO) form, z∈{0,1} (z) + (z) + (z),
min where penalizes compliance (i.e., negative stifness), enforces the volume budget, and aggregates manufacturability/regularization terms. The binary vector z encodes discrete cross-section choices for each member; decoded continuous areas were recovered via a mapping A = (z).
The QUBO was submitted to the D-Wave Advantage_system4.1 processor with the following key settings: • Number of reads: 100 • Anneal time per read: default (no custom anneal time specified) • Anneal schedule: LINEAR (default, no custom breakpoints) • Chain strength: auto (relative to the largest QUBO coeficient) • Embedding method: automatic minor-embedding via minorminer • Gauge / spin-reversal transforms: OFF • Post-processing: NONE (chain breaks resolved via majority vote) • QPU access timeout: 5000ms
Raw annealer samples were filtered for feasibility ( (z) ≤ * = 4000) and ranked by objective value. The top candidate ( = 1) was decoded to continuous variables and used as a warm start for a gradient-based classical optimizer (e.g., sequential quadratic programming). This hybrid loop—QPU sampling → classical refinement—was performed at each iteration.
In this section we repeat verbatim what has been done in the previous section, but including bending as
well: transversal displacement is no longer considered negligible. The full hamiltonian will include,
for each truss, a term like the second integral in the equation (
The solution plotted in Figure 3 was obtained through Simulated Annealing as implemented by D-Wave. Both displacements’ magnitudes and curvature have been greatly exagerated to make them more visible.
One final remark concerns the optimization within the Euler–Bernoulli framework. Unlike the purely axial case, the area optimization here involves a quadratic bending contribution in the areas, since the stifness matrix contains terms proportional to ∼ 2. In the second optimization step this contribution enters with a minus sign. To ensure well-posedness of the problem, the optimization hyperparameters must therefore be chosen such that the volume-constraint term dominates for large , as discussed in Appendix A.
However, since in practice the optimization acts on small variations of the areas ( → + Δ with Δ small), one can equivalently solve the linearized problem, retaining only the terms proportional to Δ.
In the near future, it would be desirable to provide a comprehensive benchmark against established classical optimization procedures for the problem addressed in the present work. Further generalizations are also under consideration, such as extensions to more refined beam models — for instance, the Timoshenko theory outlined in the Appendix B — as well as the inclusion of nonlinear efects. Moreover, the case of dynamical loads, i.e., time-dependent forces acting on the truss structure, is currently being investigated. Such scenarios are not only relevant for applications with variable loads (e.g., a person walking on a truss), but also pave the way for studying wave propagation in continuous media within a QUBO framework. We plan to report on these developments in the near future.
Declaration on Generative AI During the preparation of this work, the authors used ChatGPT in order to: Grammar and spelling check, Paraphrase and reword. After using this tool, the authors reviewed and edited the content as needed and takes full responsibility for the publication’s content.
A. A Coercivity Criterion for a Rank-One Perturbed Quadratic on the
Nonnegative Orthant
Here, we prove that our problem in the general setting as a global minimum if some coeficients are chosen properly. Consider () = − ∑︁ − =1 ∑︁ 2 + ︁( ∑︁ − )︁ 2, =1 =1 ∈ R+ , where , , > 0, > 0, and > 0. Let = diag(1, . . . , ) and = (1, . . . , )⊤. Define * := max 1≤≤
. 2 (29) (30) Criterion • If > * , then is coercive on R+ ; i.e., () → +∞ as ‖‖ → ∞ with ≥ 0 . Therefore a global minimizer exists. • If < * , then is unbounded below on R+ . • If = * , then along any direction ≥ 0 attaining the maximum in (30) the quadratic coeficient vanishes. In this borderline case, so is bounded below if the coeficient of is nonnegative for every such : () = (︀ − ⊤ − 2 * ⊤)︀ + const,
− ⊤ − 2 * ⊤ ≥ 0.
Otherwise () → −∞ as → ∞.
Fix ≥ 0, ̸= 0, and consider the ray = ( ≥ 0). Then () = − ⊤ − 2 ⊤ + (︀ ⊤ − )︀ 2 = − ⊤ + 2 ︁(
︁) − ⊤ + (⊤)2 − 2( ⊤) + 2.
Let Thus:
quadratic is maximized at a vertex of the simplex, giving Since () = () for any > 0 , impose ⊤ = 1 and maximize ∑︀ 2 over ≥ 0 . This convex () := ⊤ (⊤)2
2 = * .
decides boundedness. • If < * , there exists ≥ 0 with negative 2-coeficient, giving () → −∞. • If > * , the 2-coeficient is positive for all ≥ 0, so () → +∞ and coercivity follows. • If = * , the 2-coeficient is zero for maximizers of (), and the sign of the linear coeficient * program with bounds ≥ 0. • The minimizer need not be interior; some coordinates may be zero. However, = 0 is not optimal when all , , , > 0, since (0) = − − 2 < 0.
• Once > is fixed, the global minimizer can be found by solving system for the quadratic B. Timoshenko Beam Theory: Derivation and FEM Discretization The Euler-Bernoulli theory that we have used throughout the paper can be generalized to the Timoshenko theory. The model takes into account shear deformation and rotational bending efects. We no longer have the definig relation / = (), but and are treated really as independet variables.
We give a streamlined setup in the following.
B.1. Notation and Setup • Beam axis along , transverse deflection along , width along . • Point in the cross–section: (, , ), neutral axis at = 0. • Unknown fields (functions of ): () : transverse displacement (along ), () : rotation of the cross–section about , 0() : axial displacement of the neutral axis (optional). • Material: (Young’s modulus), (shear modulus), shear correction factor .
• Section properties: area , second moment of area about the neutral axis.
B.2. Kinematics (Timoshenko Assumption) Cross–sections remain rigid in their own plane but may rotate independently of ′() (shear deformation allowed). The displacement field (small rotations) is The strains are then given by: (, ) = 0() − (), (, ) = 0, (, ) = (). = = = ′0() − ′(), + = ′() − ().
B.3. Constitutive Relations and Section Resultants Assuming linear elasticity, = , = .
Resultants (force/moment per unit length) are given by: () = () = () = ∫︁ ∫︁ ∫︁ = ′0(), = − ′(), = (︀ ′() − () )︀ .
For pure bending, set 0 = 0 and = 0.
Let () be the distributed transverse load (positive downward): Substituting and yields ︀( ( ′ − ) )︀ ′ + = 0, −( ′)′ − ( ′ − ) = 0.
The internal virtual work is given as usual: External virtual work (for distributed load only): ∫︁ 0 ∫︁
0 int = ︀[ ′ ′ + ( ′ − )( ′ − ) ︀] . ext =
() + (end force/moment terms).
Set int − ext = 0 for all admissible variations to obtain the weak form. (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41)
(42) (43) B.6. Two-Node Timoshenko Beam Element (1D FEM) Consider an element of length with nodes 1 and 2. Use linear (0) shape functions: 1() = 1 − , The interpolation of and really are independent: Derivatives (constant within the element for ′) are easily computed: The nodal DOFs are still: Bending Part =
⎥⎥ . 2]︀ . =
[︀ − ︀] [︀ − ︀] . = ⎢
− ⎢
. = 12 2
is
1 3 1 + ⎣⎢−12 ⎥⎥ .
Define whose matrix form is: Then Carrying out the integrals, we get: = 1 ⏟ A compact equivalent using = Mass Matrix (Optional) For density : and similarly for rotational inertia using with .
B.7. Assembly and Boundary Conditions Assemble the global stifness and force vectors in the usual way, apply essential boundary conditions on and/or , and solve
Kd = f .
Setting = 0 ⇒ ′ = recovers Euler–Bernoulli kinematics. In FEM this corresponds to → ∞ or directly constraining ′ = with suitable interpolation.