It is shown how the differential equations of linear elasticity can be replaced by a completely equivalent set of governing algebraic equations . Those general algebraic equations can be reduced to a set of equilibrium equations for the displacements that i s an algebraic version of Navier’s differential equations. The governing algebraic equations can also be used to define the compatibility equations in terms of a compatibility matrix that is the counterpart of a differential compatibility operator matrix found in the theory of linear elasticity . Remarkably, this algebraic compatibility matrix can be obtained from homogeneous solutions of the equilibrium matrix, thus enabling the stresses (or corresponding generalized forces) to be obtained directly, without introducing displacements. Unlike earl ier force -based methods this second- generation force - based approach can be used effectively at all levels of engineering analysis, including basic engineering courses . As an example, the forces in a statically indeterminate truss are obtained using simple extensions of the statically determinate truss methods found in a first course in engineering statics . This example also shows that when second generation force -based methods are implemented in software such as MATLAB ®, solutions can be obtained without any programming, making the methods ideal tools for learning how to solve equilibrium problems in an engineering curriculum. 1 Professor Emeritus, Department of Aerospace Engineering, Iowa State University, Ames IA 50010. E -mail: lschmerr@iastate.edu Introduction Engineers have long sought to analyze the equilibrium of structures to determine the forces present. With the developm ent of the equations of linear elasticity for small deformations, the requisite equations were available in differential form for the stresses and displacements . However, the inability to solve these differential equations except in simple cases led to the development of algebraic methods that relied instead on work -energy principles . Such formulations received a major boost in the 60’s with the development of displacement -based finite element methods. In fact, the displacement- based finite element method has largely overshadowed earlier force- based algebraic methods based on work- energy principles because of its generality and its suitability for being implemented with matrices on computers . A comprehensive and fascinating history of the various approaches used to solve equilibrium problems is given in Kurrer (2018). In this paper we will show how a second- generation force -based method can once again be at the forefront in the analysis of engineering structures and how the method can be used at all levels in an engineering curriculum, from the simple equilibrium problems considered in an engineering statics course to force -based finite element solutions. This force -based method will be shown to be based on a general formulation where a general set of governing algebraic equations of linear elasticity can be developed that directly corre spond to all the differential equations and differential relations found in the theory of elasticity. These governing algebraic equations, like the traditional displacement -based finite element method, will be developed from work -energy principles but wher e stresses and displacements are treated on an equal footing instead of a purely displacement -based approach, as found in most discussions of finite elements. From these general equations one can then develop reduced formulations that involve only the stresses (or equivalently the forces) or only displacements. These reduced formulations are also algebraic equivalents of their differential equation counterparts. The second -generation force -based method that emerges is notable for t hree reasons: 1) Forces ra ther than displacements are often the quantities of most interest in an engineering analysis. 2) This force -based method, unlike earlier force methods, is as suitable for computer implementation as are displacement -based methods. 3) The force- based formulation can serve as the basis for teaching engineers how to solve equilibrium problems directly for the forces present at all levels in a curriculum. Second -generation f orce-based methods have been used as an alternative finite element method by Patniak and co -authors (see Patniak et. al. 1990, 2004, and Kaljevic´ et. al. 1996a, 1996b, 1996c, where they are called the integrated force method ). They have also been used for strength of materials problems (e.g., Patniak and Hopkins 2004, and Schmerr 2021). We will not cover those aspects of the force-based approach here. Instead, we will first show how both force-based and displacement- based methods are really two different subsets of the general algebraic equations of linear elasticity . Then we will use a truss problem example, one of the classical structural problems that students encounter in an engineering statics course, to illustrate the force -based method but where we will consider a statically indeterminate truss, a problem all current statics texts do not cover because of the (mistaken) belief that such statically indeterminate problems are too complex to be included in statics (see Schmerr (2023), which is a website with materials in support of a forthcoming engineering statics text that will apply these force -based methods). Work-Energy Relations The differential equations for a linear elastic medium can be used prove the principle of virtual work, which in matrix-vector form can be written as (Schmerr, 2021): { } { } { } { } { }( ){ }T TT VVSe dV u f dV u T dSδσ δ δ = + ∫∫∫n (1) where { }T xx yy zz yz xz xyσ σσσσσσ= (2a) { }T xx yy zz yz xz xy e eee γγγ = (2b) { }T xyz u uuu= (2c) are the stress components, engineering strain components, and displacement components, respectively, written as column vectors and where { }uδare virtual changes of the displacements and { }eδare corresponding virtual changes of the strains. The vector ( ){ }Tnrepresents the stress vector acting on the surface S of a volume V whose outward unit normal is n. The body force vector { }fis { }T xyz f fff= (3) and for a linear elastic medium the stresses and strains are related as: { }[ ]{ }Ceσ= (4) where the matrix [C] is the 6 ×6 matrix of elastic constants. One can view this principle of virtual work as equivalent to satisfying the differential equations of equilibrium. In a similar fashion, we can use the differential equations to prove the princi ple of complementary virtual work given by (Schmerr 2021): { } { } { } { }( ){ }{ }TTT VV Se dV f u dV T u dSδσ δ δ = + ∫∫∫n (5) where { }δσare virtual changes of the stresses and { }( ){ },fTδδnare virtual changes of the body forces and surface stress vector, respectively. The principle of complementary virtual work is equivalent to satisfying the differential strain -displacement relations. Now, like the approach used in a standard finite element analysis, let us express the displacements at a point x in the volume V as ( ){ } ( ) { } u HU=xx (6) where ( )Hxis a matrix of “shape” functions and { }Uis a vector of discrete nodal displacements . In contrast to the standard finite element approach, we will also express the stresses using separate shape functions ( )[]Nx, i.e., ( ){ } ( ) { }NFσ=xx (7) in terms of discrete parameters, { F}, which we will call generalized forces since, for example, in strength of materials problems they typically represent either forces or moments. Unlike the displacements { }U, the generalized forces { F} are not defined at nodes but represent parameters that describe the stresses in the volume V . In the differential equations of linear elasticity, the stresses satisfy the equilibrium equations [ ]{ } { } { } 0E fσ+= (8) where the differential equilibrium matrix operator is [ ]/ 000/ / 0/ 0/ 0/ 00/ / / 0Ex zy yzx zyx∂∂ ∂∂ ∂∂  = ∂∂ ∂∂ ∂∂ ∂∂ ∂∂ ∂∂  (9) The stress shape functions ( )[]Nxwill also be chosen so that they represent solutions of the homogeneous equilibrium equations, i.e., [ ][ ] 0EN= . One can use, for example, stress functions to define these shape functions so this is not a major restriction (Patniak and Hopkins, 2004, Schmerr 2021). The principle of virtual work then becomes { }[ ] [ ] { } { } [ ]{ } { } [ ]( ){ }T TT T TT V VSU B N dV F U H f dV U H T dSδ δδ = + ∫ ∫∫n (10) where [ ] [ ] [ ]T E BH= (see Eq. (30)). Since the variations of the virtual displacements are arbitrary, we have [ ] [ ] { }[ ]{ }[ ]( ){ }T TT V VSB N dV F H f dV H T dS = + ∫ ∫∫n (11) which we will write as [ ]{ } { } { } { }s EF Q P P=+= (12) where [ ] [ ] [ ] { }[ ]{ } { }[ ]( ){ }T V T V T S SE B N dV Q H f dV P H T dS= = =∫ ∫ ∫n (13) The matrix [ E] is the equilibrium matrix for a set of equilibrium equations, the {Q } and {P S} represent known applied body forces and surface tractions at the nodes, and {P} is just the total known force vector. In the principle of complementary virtual work, Eq. (5), the virtual changes of the stresses also satisfy the homogeneous equilibrium equations so { } 0 fδ= and the principle can be written as { } { }( ){ }{ }T VSe dV T u dSδσ δ=∫∫n (14) Using the stress- strain relations in the form { }[ ]{ } eDσ= (15) where [ ] [ ]1DC−= is the elastic compliance matrix, and placing the shape functions into Eq. (14), we find: { }[ ] [ ][ ] { }( ){ }[ ]{ }TT T VSF N D N dV F T H dS Uδδ = ∫∫n (16) But since the variations of the stresses { }δσsatisfy homogeneous equilibrium equations where { } 0 fδ=, from Eq. (11) we have: [ ]{ }[ ] [ ] { }[ ]( ){ }TT VSE F B N dV F H T dSδ δδ= = ∫∫n (17) and placing the transpose of this result into Eq. (16) gives: { }[ ]{ } { } [ ]{ }T TTF GF F E Uδδ = (18) where we will call [ ] [ ] [ ][ ]T VG N D N dV=∫ (19) the flexibility matrix. Since the variations of the generalized forces are arbitrary, we have: [ ]{ }[ ]{ }TGF E U= (20) If we define a set of deformations (generalized strains), { ∆}, as { }[ ]{ }GF∆= (21) then Eq. (20) becomes { }[ ]{ }TEU∆= (22) Algebraic Equations of Linear Elasticity We now have a complete set of equations from work-energy principles by combining Eq. (12), Eq. (21), and Eq. (22) . They are the algebraic equations of linear elasticity: [ ]{ } { } { }[ ]{ } { }[ ]{ }TEF P EU GF= ∆= ∆= (23) These are completely analogous to the governing differential equations of linear elasticity, which are the equilibrium equations, the strain-displacement equations, and the stress- strain relations: [ ]{ } { } { } { }[ ]{ } { }[ ]{ }0E T Ef eu Ceσ σ+= = =  (24) In structural problems t he forces { F} can be generalized forces (forces or moments) and the displacements can be generalized displacements (displacements or rotations). The equations of Eq. (23) can form the foundation for developing both force- based and displacement -based finite element methods. However, they have even a broader significance since the equations of Eq. (23) are the basis for solving any statically indeterminate problem involving discrete forces and displacements where the deformations are small, and the deformable elements have linear elastic behavior. Equation (24) is a complete set of 15 equations for 15 unknowns (six stresses, six strains, and three displacements). Likewise, Eq. (23) gi ves a complete set since if we have N equilibrium equations for M unknown forces, the deformation- displacement equations are M equations for the unknown M deformations and N unknown displacements, and the deformation- force equations provide M additional equations (with no additional unknowns). Thus, we have a total of 2M+N equations for 2M+N unknowns. Displacement -based Formulation The governing algebraic equations for linear elastic structures, Eq. (23), can be used to formulate a set of equations only in terms of displacements. If one writes the deformation -displacement relationship in terms of the forces as { }[ ]{ }[ ] [ ]{ }11 TF G G EU−−= ∆= (25) then the equilibrium equations become [ ][ ] [ ] { } { }1 TEG E U P−= (26) One can compare Eq. ( 26) with the displacement equation that results from the governing differential equations by placing the stress- strain and strain -displacement relations into the equilibrium equations, yielding: [ ][ ][ ] { } { } { } 0T EEC uf +=  (27) This is essentially Navier’s equations (for a general anisotropic elastic material) where there are three equations for the three displacements. In Eq. (27) if there were N equilibrium equations then this equation is N equations for N displacements. We can write Eq. (26 ) in terms of a “pseudo” stiffness matrix pK, where: [ ][ ] [ ]1 T pK EG E−= (28) giving { } { }pKU P= (29) We call pKa pseudo stiffness matrix since in the ordinary displacement -based finite element analysis the stresses and strains appearing in the principle of virtual work, Eq. (1), are written instead in terms of only the displacement shape functions and their derivatives, i.e., ( ){ } ( ) { } ( ){ }[ ]( ){ }[ ]( ) { } ( ) { }TT EEu HU e u H UB U= = = =   xx x x xx (30) where ( )Bxis the matrix of differentiated displacement shape functions. Placing these expressions into the principle of virtual work, Eq. (1), gives { }[ ] [ ][ ] { } { } [ ]{ } { } [ ]( ){ }T TT T TT V VSU B C B dV U U H f dV U H T dSδ δδ = + ∫ ∫∫n (31) which must be satisfied for all virtual variations of the nodal displacements, leading to a set of linear algebraic equations : [ ]{ } { }KU P= (32) where [ ] [ ] [ ][ ]T VK B C B dV=∫ (33) is the stiffness matrix. This is obviously quite different in form from that of the pseudo stiffness matrix, Eq. (28), which is formed from the equilibrium matrix, Eq. (13), and its transpose. However, in strength of materials problems such as the axial extension and torsion of bars and the bending of beams the two stiffness matrices are identical (Schmerr 2021). Either of these stiffness matrices can be used to represent an element in a finite element method where many elements can be assembled into a global stiffness matrix and used to solve very complex problems. It is interesting to note that while the pseudo stiffness matrix is directly related in form to Navier’s equations the ordinary stiffness matrix is not. Second -Generation Force-based Formulation In statically determinate problems the number of unknown forces is equal to the number of equilibrium equations present so we can solve for the unknown forces directly from equilibrium. But, as is often the case, if we have more unknown forces than equations of equilibrium, then the problem is statically indeterminate, and we need to augment the equilibrium equations with additional equations. Let M equal the number of forces and let N < M be the number of equilibrium equations. The deformation-displacement relations then represent M equations for the deformations in terms of N displacements so there must be M −N relations between the deformations. These M −N relations are the compatibility equations, which we will write as [ ]{ } { } 0 S∆= (34) where [ S] is the (M −N)×M compatibility matrix. If we consider the equilibrium equations in conjunction with the M −N compatibility equations, Eq. (34), written in terms of the forces we have [ ]{ } { } [ ][ ]{ } { } 0EF P SGF= = (35) which are M equations for the M unknown forces. We can combine these two sets of equations into a single set of “system” equations [ ] [ ][ ]{ }{ } { }0EP F SG   −= −     (36) which we will write as the solvable system: { }{ } sys sys EFP = (37) The key, of course, is to be able to determine the compatibility matrix [ S]. The second -generation force -based m ethod has a uniquely efficient way for obtaining this matrix. From the deformation- displacement relations we must have [ ]{ }[ ][ ]{ } { } 0TS SE U∆= = (38) which can be satisfied for arbitrary displacements if [ ][ ] [ ] 0TSE= (39) Taking the transpose, we also find [ ][ ] [ ] 0TTES= (40) Equation (40) says that the columns of the transpose of the compatibility matrix are force solutions { }hFof the homogeneous equilibrium equations, i.e. [ ]{ } { } 0h EF= (41) and there will be M −N such solutions present in the compatibility matrix [ S]. Finding such homogeneous solutions is straightforward so we can form up and solve Eq. (36) to obtain the forces directly. Equation (36) is a force-based method for solving statically indeterminate problem s. We only need the equilibrium matrix, [E], and the flexibility matrix, [G], to determine the forces since we can determine the compatibility matrix, [S], from homogeneous solutions of the equilibrium equations. In the differential equation approach the compatibility equations analogous to Eq. (38) involve second derivatives of the strains and we have (Schmerr 2021): [ ]{ } { } 0Se= (42) where [ ]22 22 2 22 22 2 22 22 2 2 2 2 22 2 22 2 2 2 2 22 20 / // 0 0 / 0 / 0/ 0 / / 0 0 0/ 0 0 2/ / / / 2/ 0 0 / / / 0 2/ 0 / / /Sz y yz z x xz y x xy xy xz yz z yz x xy xz xz xy y yz − ∂∂ − ∂∂ ∂∂ ∂ − ∂∂ − ∂∂ ∂∂ ∂ − ∂∂ − ∂∂ ∂∂ ∂=∂∂ ∂ − ∂∂ ∂ − ∂∂ ∂ ∂∂ ∂∂ ∂ ∂∂ − ∂∂ ∂ − ∂∂ ∂ ∂∂ ∂ − ∂∂ ∂ ∂∂ − ∂∂ ∂  (43) If we combine the compatibility equations with the equilibrium equations, we find [ ]{ } { } { } [ ][ ]{ } { }0 0E Sf Dσ σ+= =  (44) where [ ] [ ]1DC−= is the elastic compliance matrix. The equations in Eq. (44) are in terms of the stresses only and are counterparts to the algebraic equations of Eq. (35) . This stress- based formulation, however, is nine equations for the six stresses so unlike Eq. ( 35) we have too many equations in Eq. (44) for a well -posed stress formulation. One way to solve this problem is to express the stresses in terms of stress potentials that automatically satisfy the equilibrium equations leading to an equal number of equations and unknowns. The Airy stress function is a well-known example of this type of function that is used to solve 2- D elasticity problems (Kaljevic´ et. al. 1996b, Schmerr 2021) . Note that there are also differential relations between the compatibility matrix operator and equilibrium matrix operator in the form of Eq. (39) and Eq. (40) given by [ ][ ] [ ] [ ][ ] [ ]0 0T SE TT ES== (45) so that all the differential equations or relations appearing in a theory of linear elasticity have algebraic counterparts. The equilibrium matrix [ E] can be used to represent an element in a finite element method where many elements can be assembled into a global equilibrium matrix. The compatibility matrix [S] can then be obtained from that global matrix and combined with the equilibrium equations in the same form as seen in Eq. (35) to solve very complex problems directly for the forces . Example : A Sta tically Indeterminate Truss As described in the Introduction t he second -generation force -based method has been successfully used in a finite element solution procedure in place of the traditional displacement- based approach . However, what has not been widely recognized is that it can also be used as a general method for solving any statically indeterminate problem involving linear elastic behavior. This makes the method an ideal tool for teaching how to solve equilibrium problems. For example, the forces present in statically indeterminate trusses can be obtained with procedures that are simple extensions of those normally considered for statically determinate problems in an engineering static s course (Schmerr 2023). We will give an example of such a truss problem that is solved with MATLAB ®. The truss in Fig. 1 is fixed at both A and C and is loaded by the forces F = 8 kN and P = 20 kN. In a pinned truss the members all experience only axial extension or compression, so the deformation- force relatio nship is simply gF∆= , where the deformation ∆ is the elongation and g is the flexibility given by g = L /AE, and where A is the cross -sectional area, E is Young’s modulus, and L is the member length. In this case the flexibility matrix [ G] will be a diagonal matrix containing the flexibilities along its major diagonal. Let AE be the same for all members. Determine the forces in the truss members and the reactions at points A and C. If we use the method of pins for this problem, we need to draw the free body diagrams of all the pins, as shown in Fig. 2. The eight equilibrium equations at pins A, B, C, and D are 0 3 /5 0 0 4 /5 0 0 3 /5 3 /5 8 0 0 4 /5 4 /5 0 0 3 /5 0 0 4 /5 0 0 16 0 0 12 0x x AD AB A y y AB A x AB Bc B y BD AB BC B x x C D BC C y y BC C x CD AD D y BD DF AF F F AF F FF F FF F F CF F F CF F FF FF= ++ = = += = − + += =−− − = = −− = = += = − += = −=∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ (46) There are nine unknowns, consisting of the five forces in the truss members and the four reaction forces at A and C. Thus, the problem is statically indeterminate. The compatibility matrix in this case is a vector. To solve this problem , we will use the symbolic algebra capabilities of MATLAB which allows us to write the e quilibrium equations of Eq. ( 46) directly in symbolic form : % symbolic member forces and reactions syms Fab Fad Fbc Fbd Fcd Ax Ay Cx Cy % force vector of all the unknowns T = [Fab Fad Fbc Fbd Fcd Ax Ay Cx Cy]; % equations of equilibrium at the pins Eq(1) = Ax + Fad +3*Fab/5 == 0; Eq(2) = Ay +4*Fab/5 == 0; Eq(3) = - 3*Fab/5 + 3*Fbc/5 + 8 == 0; Eq(4) = - Fbd -4*Fab/5 - 4*Fbc/5 == 0; Eq(5) = Cx - Fcd -3*Fbc/5 == 0; Eq(6) = Cy + 4*Fbc/5 == 0; Eq(7) = Fcd - Fad + 16 == 0; Eq(8) = Fbd - 12 == 0; Once we have these symbolic equations , we can use them to define the equilibrium equations in matrix form [E ]{T} ={P}. The MATLAB function equationsToMatrix gives us the equilibrium matrix [E] and the vector {P}: [E, P] = equationsToMatrix(Eq, T); Having [E ] we th en can use the MATLAB null function to obtain the transpose of the compatibility matrix [ S]T by solving the homogeneous equilibrium equations. Taking the transpose then gives us the compatibility matrix [ S] itself : S = null(E).'; Next, we need to def ine the flexibility matrix [ G]. We can omit a common flexibility term L r/EA for all the truss members, where we will let a reference length be L r = Lab = 5 m since such a common factor can always be eliminated from the homogeneous compatibility equations. Then the remaining terms in the flexibility matrix for the members are just ratios of the member lengths to Lr. We set the flexibilities of the reactions at A and C equal to zer o since there are no deformations associated with these forces at the fixed pins. % form up a vector of normalized flexibilities and zeros flex = [ 1 3/5 1 4/5 3/5 0 0 0 0]; % and place them on the diagonal of [G] with the MATLAB diag function G = diag(flex); We then combine the equilibrium and compatibility equations and form the elements of the system equations [ Esys]{F} = {Psys}. We solve for the unknown forces , where the components in {F} are the values of the symbolic tensions contained in {T} : Esys = [ E; S*G]; Psys = [P; 0]; F = Esys\ Psys; The vector F is still in symbolic form, so we convert it to a double precision numerical vector with the MATLAB function double: double(F) ans = -0.8333 % Fab = - 0.8333 kN 8.0000 % Fad = 8 kN -14.1667 % Fbc = - 14.1667 kN 12.0000 % Fbd = 12 kN -8.0000 % Fcd = - 8 kN -7.5000 % Ax = - 7.5 kN 0.6667 % Ay = 0.6667 kN -16.5000 % Cx = - 16.5 kN 11.3333 % Cy = 11.3333 kN You can verify that the reactions do satisfy equilibrium for the entire truss. Setting the flexibilities to be zero at the rigid supports is much like the approach used in displacement -based finite element methods where very large stiffness values (corresponding t o very small flexibilities) are used to approximately satisfy zero displacements at the supports. However, in the force -based method the flexibilities can be set to be identically zero . This second -generation force -based method is simple and direct in cont rast to obtaining the compatibility equations for a truss by displacements (rotations) (e.g., Williams and Dym 2005). Discussion and Conclusions It has been shown that the differential equations of linear elasticity, which have been used for solving problems since the late 19th century, can be replaced by a n equivalent set of algebraic equations of linear elasticity. These algebraic equations are then used to develop formulations based either on displacements or forces. The algebraic equations of linear elasticity form the basis for solving all types of problems, from the simple problems found in statics to complex problems using force -based or displacement -based finite element methods . The emphasis here has been on a force -based approach because of its many advantages that have already been discussed. It has also been shown how the symbolic algebra and numerical capabilities of software such as MATLAB can be used to obtain second -generation force -based solutions directly without programming by using MATLAB as a sophisticated “calculator” containing all the built -in functions needed to solve equilibrium problems . Thus, this second -generation force - based approach is also ideally suited for teaching students how to solve equilibrium problems in a modern and efficient manner .