Two-Dimensional Stress Analysis of Isotropic Deep Beams Using the Finite Difference Method

: This paper presents an approach to the two-dimensional analysis of elastic isotropic deep beams using the finite difference method (FDM). Deep beams are subjected to in-plane loading and present a shear span to height ratio less than 2.50; consequently, Euler − Bernoulli beam theory and Timoshenko beam theory do not apply. Deep beams analysis is generally conducted using numerical methods such as the finite element method and to a lesser extent the FDM; the strut-and-tie model and the stress field method are also widely utilized. Analytical approaches usually make use of the Airy stress function, where stresses are formulated in terms of the stress function; however, the exact solution of this function satisfying all of the boundary conditions can hardly be found, even for simple cases. In this paper, deep beams were analyzed using the FDM. The FDM is an approximate method for solving problems described with differential equations. The FDM does not involve solving differential equations; equations are formulated with values at selected nodes of the structure. Therefore, the deep beam was discretized with a two-dimensional grid, and additional nodes were introduced at the boundaries and at positions of discontinuity (openings, brutal change of material properties, non-uniform grid spacing), the number of additional nodes corresponding to the number of boundary conditions at the node of interest. The introduction of additional nodes allowed us to apply the governing equations at boundary nodes and to satisfy the boundary and continuity conditions. An Airy stress function approach and a displacement potential function approach were considered in this study whereby strong formulations of equations (equilibrium, kinematic, and constitutive) were set. Stress and stability analyses were carried out with this model; furthermore, deep beams of varying stiffness, layered beams, and beams having openings were analyzed. For slender beams, the results obtained with the Airy stress function approach showed good agreement with those of Euler − Bernoulli beam theory, and for deep beams, the shapes of stress distributions were in good agreement with a proper understanding of the behavior of structures. On the other hand, the displacement potential function approach delivered unsatisfactory results, probably due to the use of an inefficient equation solver; a more powerful tool will be needed in future research for this purpose.


Introduction
This paper describes the application of Fogang's model [1] based on the finite difference method (FDM), used for the Euler−Bernoulli beam, to the isotropic deep beam. Deep beams are essentially subjected to in-plane loading and present a shear span to height ratio is less than 2.50; consequently, Euler−Bernoulli beam theory and Timoshenko beam theory do not apply. The analytical approach to deep beam analysis is based on the Airy stress function method; this method developed in 1862 by G.B. Airy [2] consists of introducing a new function, the Airy stress function. The stresses were formulated in terms of this new function and a new differential equation, a biharmonic equation, was obtained. The problem of determining the stresses was then reduced to that of finding the stress function solution of the biharmonic equation which satisfies the boundary conditions. Neau [3] developed a scheme for applying doubly infinite power series to the Airy stress function for isotropic bodies; problems in which boundary stresses can be described by means of power series are solvable by this method. Jayne et al. [4] solved the characteristic fourth-order partial differential equation for two-dimensional elastic anisotropic and orthotropic materials, using a doubly infinite power series. Ahmed et al. [5] used an ideal mathematical model, based on a displacement potential function, to formulate the problem; displacements and stresses were formulated in terms of this potential function and a new differential equation, a biharmonic equation, was obtained.
However, the exact solution of the differential equation satisfying the boundary conditions can hardly be found, even for simple cases. Numerical methods permit therefore to overcome solving the differential equations. Ismail et al. [6] performed a series of nonlinear finite element to evaluate the different design approaches available in the literature for design of reinforced concrete deep beam with large opening; three finite element models were developed and analyzed using the computer software ATENA.Vilar et al. [7] proposed a numerical solution to deep beams using the layerwise displacement theory; a finite element solution for deep beams based on a layerwise displacement field considering the full stress/strain tensors was provided. Sri Harsha et al. [8] gave the analytical investigation of reinforced concrete deep beams reinforced with horizontal and vertical web reinforcement; a formula using nonlinear finite element method by ABAQUS was proposed to define the shear strength of deep beams.
For the design of structural concrete, Schlaich et al. [9] proposed a generalization of the truss analogy in order to apply it in the form of strut-and tie-models to every part of any structure; it was described how strut-and tie-models, which condense all stresses in compression and tension members and join them by nodes, can be developed by following the path of the forces throughout a structure. Liu et al. [10] proposed a model for deep beams with rectangular openings that stems from a two-parameter kinematic theory for solid beams; the model was established based on an analysis of the shear behavior and failure modes of test specimens using nonlinear finite element and strut-and-tie models. Silveira et al. [11] proposed a solution based on the stress field method for the analysis, design, and detailing of deep beams. The stress field method, an alternative method to strut-and tie-method for concrete structures subjected to discontinuities, consists of finding the stresses acting in discrete area elements whereby non-linear elastic-plastic stress fields are used.
De Mello et al. [12] presented the stringer-panel method, an alternative procedure to strut-and-tie method; the structure is divided on two distinct elements: stringers, which absorb normal forces, and panels, which absorb shear forces by membrane action. Then the overall structure behavior is investigated by means of non-linear analysis.  In the classical analysis using the FDM, points outside the deep beam are not considered; the boundary conditions are   applied at the beam's boundaries and not the governing equations. Consequently, the non-application of the governing   equations at the boundaries have led to inaccurate results, making the FDM less interesting in comparison to other numerical methods such as the finite element method. In this paper, a model based on FDM was presented. This model consisted of formulating the differential equations with finite differences and introducing additional nodes outside the beam and at positions of discontinuity (openings, brutal change of stiffness's, non-uniform grid spacing). The introduction of additional nodes allowed us to apply the governing equations at the boundaries and to satisfy the boundary and continuity conditions. An Airy [2] stress function approach and a displacement potential function approach of Ahmed et al. [5] were considered whereby strong formulations of equations (equilibrium, kinematic, and constitutive) were set. In the Airy stress function approach stresses were formulated in terms of the Airy stress function whereas in the displacement potential function approach displacements and stresses were formulated in terms of the potential function.

Equations of the theory of elasticity
In this section the equations of the theory of elasticity are recalled. Figure 1 represents a deep beam with the axis convention (X, Y).

Figure 1 Deep beam with axis convention X, Y
Displacements in x-and y-directions are denoted by u(x,y) and v(x,y), respectively. The axial strains ε xx and ε yy and the shearing strain γ xy are related to the displacements as follows: Combining Equations (1a-c) yields the following compatibility condition between axial strains and shearing strain where D is the axial rigidity of the deep beam. The solutions considered in the present study and presented thereafter involve a two displacement function, an Airy stress function, and a displacement potential function.

Two displacement function
The two displacement function approach is governed by Equations (5a-b) which are applied at any node of the structure whereby the displacements are the unknowns. The geometric boundary conditions are directly formulated whereas stress related boundary conditions are satisfied using Equations (3a-c).

Airy stress function
Combining Equations (1a-b), (2a-b), and (6a-b) yields the following relationship between the Airy stress function (ASF) and the rates of change of displacements The analysis is then reduced of determining one single function, the ASF. However, satisfying the geometric boundary conditions is not easy with this solution since the displacements are not specified in terms of the stress function. If constant body forces p x and p y are applied, Equations (6a-c) can be modified using Säckel [13] as follows (8c)

Displacement potential function
In case no body forces were applied on the structure, Ahmed et al. [5] introduced a displacement potential function (DPF) ψ(x,y), defined in terms of displacement components. However, the Cartesian coordinates x, y and the displacements u(x,y) and v(x,y) of Ahmed et al. [5] have to be inverted in order to be consistent with the axis convention of the present study. Therefore, the displacements are related to the potential function are as follows  Similarly to the solution using the ASF, the analysis is reduced of determining one single function, the DPF; interestingly, the satisfaction of boundary conditions is facilitated in this solution, since displacements and stresses are specified in terms of the potential function.
If constant body forces p x and p y are applied, Equations (9a, b) of Ahmed et al. [5] can be modified as follows The expression of the shearing stress is unchanged whereas that of axial stresses becomes

Finite difference approximations for a deep beam 2.2.1 Fundamentals of finite difference approximations
The two displacement function approach will not be considered further for the stress analysis since two equations (the governing equations (5a-b)) have to be set at any node, whereas only one equation is set at any node for the solutions using the Airy stress function (ASF) and the solution using the potential function. First, Equation (7) is the governing equation for the ASF approach. This equation has fourth order derivatives; consequently, the stress function φ(x,y) is approximated around the node of interest i as a fourth degree polynomial in each direction. The unknown at any node being the value φ i of the stress function, the corresponding finite difference approximation (FDA) is denoted by ASF−FDA. Second, Equation (9c) is the governing equation in case the DPF ψ(x,y) is considered. The unknown at any node being the value ψ i of the potential function, the corresponding FDA is denoted by DPF−FDA

Airy stress function finite difference approximation
Given the grid spacings ∆x = h and ∆y = λh in x-and y-direction, respectively. The stress function φ(x,y) is approximated around the node of interest i as a fourth order polynomial in each direction; however, for simplification purpose, the first and second partial derivatives in x-direction and the mixed partial derivative ∂ 2 φ/∂x∂y are expressed using a second order polynomial hypothesis for φ(x,y). The FDAs are then given by Stress Analysis of Isotropic Deep Beams Using the Finite Difference Method In the stencil notation the factor associated to the node of interest is in brackets. The partial derivatives in y-direction are formulated similarly. The FDAs of the fourth order partial derivatives in x-direction are given by The FDA of the term ∂ 4 φ/∂x 2 ∂y 2 defined using Equation (10a) is expressed by means of the following stencil (11a)
Using a fourth order polynomial hypothesis for ψ(x,y), the third derivatives ∂ 3 ψ/∂x 3 at different nodes are expressed in terms of values of ψ(x,y) as follows Using Equations (10a, c), the FDAs of the partial derivatives ∂ 3 ψ/∂x∂y 2 and ∂ 3 ψ/∂x 2 ∂y are given by   2  2  2  1  1  2  2  2  2  2  1  1   1  1   10  1  2  1  ,  000  4 10 1 Stress Analysis of Isotropic Deep Beams Using the Finite Difference Method Figure 2 shows the node distribution of a deep beam having equidistant nodes with spacings ∆x and ∆y in x-and ydirection, respectively. The node of interest k and the surrounding nodes are represented, whereby n, s, e, and w stand for the directions north, south, east, and west, respectively, according to the directions in the stencil. The node k may even be at the boundary of the beam, however being not at angles.

Figure 2 Point of interest k and its surrounding points
The FDAs for the Airy stress function approach and the displacement potential function approach are determined.

Airy stress function FDA at an interior node
Given the grid spacings ∆x = h and ∆y = λh. The governing equation (Equation (7)) at a given node can be expressed by means of a stencil using Equations (11) and (11a) as follows (13) [ ] The axial/shearing stresses at a given node are expressed using Equations (6a-c) and (10a-b) as follows; Introducing following modified displacements, the displacements are calculated using Equations (8a-b), (10a, c) and (15) as follows

Displacement potential function FDA at an interior node
The governing equation (Equation (9c)) can be expressed by means of a stencil using Equations (13). At node i the FDA of u(x,y) and v(x,y) are formulated using Equations (9a-b) and (10a-b) as follows The FDA of the axial/shearing stresses are formulated using Equations (9d-f) and (12a-b) as follows (14) (       (7)) of the beam at node k is then described with the following 13-point stencil For other angle points of the structure, the governing equations are determined using the same principle, and are presented in Appendix A. The axial/shearing stresses are expressed using Equations (14); therefore the stress related boundary conditions are formulated.
In case of geometric boundary conditions i.e. u = 0 at the node k, u-displacements are introduced at nodes k, e, and ee; for v = 0 at the node k, v-displacements are introduced at nodes k, n, and nn, as represented in Figure 4.

Governing equation, moments, and Kirchhoff shear forces at nodes e, s
The equations of interior nodes applied. (20)

Displacement potential function in the vicinity of beam angles
The node distribution is the same as that of the ASF-FDA, as shown in Figure 3.

DPF Finite difference approximation at node k
The FDA of the governing equation (Equation (9c)) is expressed with Equation (19), and that of u(x,y) and v(x,y) using Equations (17). The FDAs of the axial/shearing stresses are formulated using Equations (9d-f) and (12a-b) as follows For other angle nodes of the structure, the governing equations are determined similarly, and are presented in Appendix A. The axial/shearing stresses expressed using Equations (12a-b) are displayed in Appendix B.

Governing equation, moments, and Kirchhoff shear forces at nodes e, s
The equations of interior nodes applied.   Figure 5 shows a skew edge with regular and additional nodes. The outer normal n to the skew edge makes an angle α with the +x-axis. Two additional nodes are associated to each edge node; therefore, governing equations can be applied at the edge nodes and boundary conditions be satisfied. The governing equations (Equation (7) or (9c)) can be expressed using Equation (13). Let the axial stresses normal to the skew edge and the shearing stresses tangential to the skew edge be denoted by σ n and τ t , respectively. The equilibrium equations on an infinitesimal edge element yield the following widely known relationship between σ n /τ t and the axial/shearing stresses in Cartesian system

ASF−FDA and DPF−FDA at skew edges
The axial/shearing stresses in Cartesian system are expressed using Equations (14) for the ASF approach and Equations (18) for the DPF approach. The boundary conditions depend on the displacements and the axial/shearing stresses σ n and τ t ; at a free skew edge i.e. they are given by σ n = 0 and τ t = 0.

Finite difference approximations of loadings
Let us determine the FDA of the distributed load in the case of a varying distributed load per unit area q(x). The FDA, denoted by q i , can be taken as the average load around the node of interest and is then expressed as follows: The load q i is used to satisfy the boundary conditions, namely σ yy = -q i at the node of interest.

Analysis at positions of discontinuity
Positions of discontinuity are positions of application of concentrated forces, supports, openings, and springs.

Concentrated force P at node i
The concentrated load P acting in y-direction can be converted into a load per unit area qi at the node of interest (25) d being the thickness of the beam. The boundary condition σ yyi = -q i is applied at the node of interest.

Airy stress function FDA
In case of a support, the boundary conditions (u = 0, v = 0) are satisfied using Equations (20) and (21). In case of a concentrated spring of stiffness K W acting in y-direction i.e, the boundary condition is given by and in case of an elastic Winkler foundation of stiffness k W in y-direction i.e, the boundary condition is given by The axial stress σ yyi is calculated using Equation (14b), and the displacement v i using Equations (21a-c).

Displacement potential function FDA
In case of a support, the boundary conditions (u = 0, v = 0) are satisfied using Equations (17). In case of a concentrated spring or an elastic Winkler foundation Equations (26a-b) applied further, whereby the axial stress σ yyi is calculated using Equation (18b) and the displacement v i using Equations (17b).

Local grid refinement
A local grid refinement can be considered at positions where a concentrated load is applied, or at a concentrated support or spring. Refinement nodes (ο) are then introduced around the node of interest k, as represented in Figure 6 for ASF−FDA and DPF−FDA. Furthermore, additional nodes (×) are introduced and the FDAs and derivatives are obtained using Lagrange interpolating polynomials (Equation (27b)).

Deep beams with openings
A deep beam having an opening is represented in Figure 7   The function F represents the Airy stress function or the displacement potential function in the following. Particular attention must be taken by the formulation of the governing equations, especially the term ∂ 4 F/∂x 2 ∂y 2 in the vicinity of angle nodes. The mixed partial derivative ∂ 4 F/∂x 2 ∂y 2 is expressed using Equation (11a), whereby for angle nodes 1 and 7 the expressions (F 2a + F 8a )/2 and (F 6a + F 11a )/2 are considered in the stencil; furthermore, at node 2 i.e. it involves the nodes 8, 2a, and 3a, and at node 8 it involves the nodes 2, 8a, and 9a.
Stress Analysis of Isotropic Deep Beams Using the Finite Difference Method

Analysis of isotropic deep beams of variable stiffness
The axial rigidity D (Equation (5c)) may vary continuously throughout the beam, and is then denoted by D(x,y).
Substituting D(x,y) into Equations (3a-c) yields the axial/shearing stresses as follows (28) The substitution of Equation (28) into (5a-b) yields the following governing equations (29) Equations (29) are developed and the corresponding FDAs are formulated as described in previous sections. The analysis continues similarly to that of beams of constant stiffness.

Analysis of layered beams
A layered beam presents an abrupt change of material properties apart from a line, as represented in Figure 8.

Stress Analysis of Isotropic Deep Beams Using the Finite Difference Method
The model developed by Fogang [1] is applied: an opening is realized in the beam along the line i (nodes i1, i2, …, i9) and additional nodes (×) are introduced in the opening, as represented in Figure 9. The Airy stress function approach and the displacement potential function approach are considered.

Displacement potential function FDA
The node distributions are represented in Figure 9.   The governing equations at angle nodes f1 and ll are modified as follows to take into account the nodes fld and lld.  where n xi , n yi , and n xyi were the axial/shearing forces per unit length at node i, D was the flexural rigidity, and d was the plate thickness. The axial/shearing forces per unit length are related to the stresses of the present study as follows (32c) Therefore, the axial/shearing forces per unit length at any node are determined according to the present study. Thereafter the buckling analysis of the deep beam, in fact a plate buckling analysis, is carried out using Fogang [14].

Two point supported deep beam and subjected to a distributed load
A deep beam resting on a fixed and a rolling support and subjected to a uniformly distributed load per unit area p =1.0 kN/m², as shown in Figure 11, was analyzed. An 8 x 8 grid was considered, such that a = 10.0m = 8∆x and b = 8∆y. The axial stresses σxx at nodes 5, 15…, and 85 and the shear stresses τxy at nodes 3, 13…, and 83 are listed in Table1   and Table 2, respectively, depending on λ = b/a = ∆y/∆x. The calculations are conducted using the Airy stress function (ASF) and the displacement potential function (DPF). Results with ASF are also displayed in graphs. Details of the results are presented in the Supplementary file "Two point supported deep beam subjected to a distributed load."

Two point supported deep beam with an opening and subjected to a distributed load
A deep beam with an opening resting on a fixed and a rolling support and subjected to a uniform load per unit area p =1.0 kN/m², as shown in Figure 12, was analyzed. An 8 x 8 grid was considered, such that a = 8∆x and b = 8∆y. Details of the calculation are presented in the Supplementary file "Two point supported deep beam with opening subjected to a distributed load." Table 3 displays the axial stresses σ xx at mid-span of the beam, at nodes represented in Figure 12.

Conclusion
The finite difference method based model developed in this paper provided a solution to the stress and stability analyses of deep beams. This model consisted of formulating the differential equations with finite differences and introducing additional nodes outside the beam and at positions of discontinuity (openings, brutal change of stiffness's, non-uniform grid spacing). The introduction of additional nodes permitted to apply the governing equations at the boundaries and to satisfy all of the boundary and continuity conditions. An Airy stress function approach and a displacement potential function approach were considered together with strong formulations of equations (equilibrium, kinematic, and constitutive). By the Airy stress function approach, stresses were formulated in terms of the stress function but geometric boundary conditions were not directly formulated; as result, stresses throughout the structure and displacements in the vicinity of supports were delivered. In the displacement potential function approach, displacements and stresses were formulated in terms of the potential function; so all of the boundary conditions, stress related and geometric, were conveniently expressed. Deep beams of varying stiffness, layered beams, and beams having openings were analyzed with the model. The results obtained using the Airy stress function approach were in agreement with a proper understanding of structural behavior; unfortunately, the displacement potential function approach delivered unsatisfactory results, probably due to the use of an inefficient equation solver.
The following aspects were not addressed in this study but could be analyzed with the model in future research: