# Developed Mathematical Model for Indeterminate Elements with Variable Inertia and Curved Elements with Constant Cross-Section

Authors: Mohamed A. El Zareef; Mohamed E. El Madawy; Mohamed Ghannam
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Issues such as analysis of indeterminate structural elements that have variable inertia as well as a curved shape still have no closed form solution and are considered one of the major problems faced by design engineers. One method to cope with these issues is by using suitable the finite element (FE) software for analyzing these types of elements. Although it saves time, utilization of FE programs still needs professional users and not all engineers are familiar with it. This paper has two main objectives; first, to develop simple mathematical models for analyzing indeterminate structural elements with variable inertia and that have a curved shape with constant cross section, this model is much easier to be used by engineers compared to the FE model. For simplicity and saving time, a MATLAB program is developed based on investigated mathematical models. The force method combined with numerical integration technique is used to develop these models. The developed mathematical models are verified using the suitable FE software; good agreement was observed between the mathematical and the FE model. The second objective is to introduce a mathematical formula to determine the accurate number of divisions that would be used in the mathematical models. The study proves that the accuracy of analysis depends on the number of divisions used in the numerical integration. The optimum number of divisions is obtained by comparing the output results for both FE and developed mathematical models. The developed mathematical models show a good agreement with FE results with faster processing time and easier usage.


---

## Body

## 1. Introduction

Beams with variable sections have a lot of practical applications in concrete and in steel structures. Tapered beams made of welded plate are an example of variable sections beams that are used in steel structures; these types of beams are used mainly in rafter and columns of frames or as continuous beams, and the main advantage of these type of section is they are economical where the beams have variable heights along their length so that the beam capacity equals the acting moment [1, 2]. In concrete structures, beams with variable sections are used in many structural applications and in many precast concrete buildings [3].Designing and analysis of variable section beams are considered one of the problems that engineers face. Due to lack of knowledge, using the mean value of stiffness or ignoring variation of stiffness along element length are used in the design process of such elements, which is not accurate and can lead to progressive collapse of the structure. In some special structures that have elements with variable stiffness, it is essential to find an accurate analysis method for these elements. One of the main issues in designing a beam with variable inertia is that the section which is subjected to maximum stress need not necessarily be at the midspan of the beam due the variation in beam depth, which leads to changing the inertia along the beam length [4].Many researchers were motivated to study the behaviour of beams with variable sections under different loading scenarios. Boiangiu et al. [5]; Ece et al. [6] and Felsoufi and Azrar [7] studied the dynamic behaviour and the free body vibration of beams with variable cross section. Yang et al. [8] and Dumitrache [9] studied the shear stress distribution and analysis for variable sections beam. Yang et al. [8] found that design codes do not accurately calculate the stress distribution in both the cracked and elastic stage using the effective shear force method.Koo Lee et al. [10], Rojas and Ramirez [11], and Sapountzakis and Tsiatas [12] investigated the elastic linear behaviour of beams with variable sections using the classical beam theory, and others studied the nonlinear behaviour of such beams [13, 14]. However, most of the proposed design equations are still complicated to be studied by the engineers.Another structural element that has a beautiful and attractive architecture view is arches. Arches were previously used by the Persian, Egyptian, Babylonian, Greek, and Islamic civilizations for different types of structures, but the ancient Romans were the first to use them as bridges. Arches have been used in some bridges in China and have been used widely in Islamic mosques, castles, and palaces.Nowadays, arches are used in different structural applications as railway and roadway bridges, pedestrian bridges, or as part of a building, hangers with a large open space without intermediate columns. The basic structural theory for arch is that it converts all stress into compressive and eliminates tensile stresses. This is useful especially for materials which carry less tensile stress.Similar to beams with variable sections, engineers face difficulties in analyzing arches and determining the critical cross section and also calculating straining action along the arch length. Due to the importance of arches in our structures, many researches study the behaviour of arches under different loading conditions with different shape profiles. Sonavane [15] applied basic principles of the flexibility method in the analysis and design of symmetrical circular arches under different loading conditions. The results obtained from the analysis were compared with the finite element (FE) model results. MATLAB program was used to determine the internal forces in the circular arch design.Ghannam and Najm [16] investigated the buckling behaviour of Islamic arches under different loading conditions. King and Brown [17] conducted the manufacturing and design of curved steel members using the British code (BS 5950–1 : 2000). King and Brown [17] included worked examples to clarify the design of arched steel used in different applications.Designers tend to use FE software as one of the most effective methods that can be used in the analysis of arches or straight beams with variable cross sections. However, there is still no exact mathematical model for solving such problems.Based on the above literature, it can be concluded that the main problem that faces engineers in designing structural elements that have variable inertia as well as that have a curved shape is the complication of the design methods that are available for these types of structures. Besides, not all engineers are familiar with using FE programs. Therefore, the main aim of this paper is to cover that gap in this area and to provide simple mathematical analytical models that would simplify the analysis and design of structural elements that have variable inertia and a curved shape. In this paper, the flexibility and numerical integration methods are used to develop simplified mathematical analytical models. These investigated models are used to evaluate the straining actions for straight beams with variable cross section and curved beams that have constant moment of inertia. To ease and save time of the solution process, simple MATLAB programs are established based on the investigated mathematical models. The output results from these models are verified using the suitable finite element software [18]. It was found that there is a good agreement between the developed program and the result of the FE program.By comparing the output results from both the FE software and developed mathematical models, it is observed that the number of divisions used in the numerical integration of the mathematical models play a main role to obtain accurate results.Therefore, the second objective of this paper is to present mathematical formulas for curved and straight beams with variable inertia to determine the optimum number of divisions used in the analytical model of each. A numerical analysis of the output results is applied to investigate these formulas. By adjusting the number of divisions used in the developed models, the output results show a good agreement with FE results.


## 2. Flexibility Method (Force Method)

The flexibility or force method was originally developed by James Clerk Maxwell in 1874. This method is very common and a powerful method for the analysis of statically indeterminate structures without limitations as in other structural analysis methods.The method is based on transforming the structure into a statically determinate system and calculating the value of redundant forces required to restore the geometric boundary conditions of the original structure. The force method (also called the flexibility method or method of consistent deformation) is used to calculate reactions and internal forces in statically indeterminate structures. The method can be used for either plane or space structures subjected to loads or forced deformations. It is well known that the analysis of beams with variable stiffness along length in a closed form solution cannot be obtained. Therefore, it is important to modify the existing method (flexibility method) for elements with constant stiffness to be suitable for beams with variable stiffness. Additionally, the exact analysis method for the curved beams with constant section is studied in this paper to fill the gap in the literature.The flexibility and the stiffness for an element at a certain point and direction are reverse. The inverse of stiffness is flexibility. The relationship between flexibility and stiffness can be expressed as follows:(1)rij=δij−1,where rij is the stiffness matrix and δij is the flexibility matrix of an element. Flexibility matrix for an element is generally determined by the following equation:(2)δij=∫Mix⋅MjxEIx⋅dl+∫Miy⋅MjyEIy⋅dl+∫Qix⋅QjxGAr⋅dl+∫Qiy⋅QjyGAr⋅dl+∫Ni⋅NjEA⋅dl+∫Miz⋅MjzGIt⋅dl.Equation (2) can be simplified according to the type of structure. For beams, equation (2) can be simplified into the following equation, where the bending moment Mx is the governing staining action for the beam cross section:(3)δij=∫Mix⋅MjxEIx⋅dl,where (Mix,Mjx) is the corresponding bending moment diagrams for force at ends (i, j), respectively, (E) is the modulus of elasticity, and (Ix) is the cross section moment of inertia. The solution of the numerical integration can be obtained using trapezoidal rules presented in the following equation:(4)δi=1/2+∑i=2n−1δ+δi=n/2n.The equations of the force method can be written in the matrix form as follows:(5)δijδijδjiδjjfifj=θiθj,where fi and fj are the corresponding member force at the ends i and j, respectively, θi and θj are the member rotation at the ends i and j, respectively, and δij is the deflection in the direction of the i redundant due to a unit load in the direction of the j redundant.


## 3. Proposed Analytical Model and Verifications

### 3.1. Application 1: Beam with Variable Inertia

#### 3.1.1. Analytical Model for Beams with Variable Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the beam due to clockwise unit rotation separately at both ends of the beam as shown in Figure1.Figure 1
The proposed model for solving beam with variable section.Givenθi=1.0 (rotation at each end of the beam separately),(6)hx=0=1.25h,hx=L/2=h,hx=L=0.75h,hx=h5−2x/L4,b=constant,where h is the depth of the beam at the midspan, hx is the depth of the beam at distance x from the left side of the beam, and b is the beam width.(7)Ix=bhx312,Ix=L/2=I=bh312,Ix=I5−2x/l364.The equations of the force method:(8)δ11δ12δ21δ22X1X2=θ1θ2,δij=∫01Mx1⋅Mx2EIx⋅dx,Mx1=1−xl,Mx2=xl,xl=ξ,Mξ1=1−ξ,Mξ2=ξ.By transformation fromx to ξ, δij equations take the form presented in the following equation:(9)δ11=∫011−ξ2⋅64lEI5−2ξ3dξ=64lEI∫01f1⋅dξ,f1=1−ξ25−2ξ3,δ12=−∫011−ξ⋅ξ⋅64lEI5−2ξ3dξ=−64lEI∫01f2⋅dξ,f2=1−ξ⋅ξ5−2ξ3,δ22=∫01ξ2⋅64lEI5−2ξ3dξ=−64lEI∫01f3⋅dξ,f3=ξ25−2ξ3.For the numerical integration of functionZi, assume the interval ξ is divided into n parts. The values of the function Zi are given in Table 1 using n=8 parts.Table 1
Values of the functionZi using n=8 parts.
ξ
f
1
f
2
f
3
0
8 ∗ 10−3
0
0
1/8
7.144 ∗ 10−3
1.021 ∗ 10−3
0.146 ∗ 10−3
2/8
6.179 ∗ 10−3
2.058 ∗ 10−3
0.686 ∗ 10−3
3/8
5.088 ∗ 10−3
3.053 ∗ 10−3
1.832 ∗ 10−3
4/8
3.906 ∗ 10−3
3.906 ∗ 10−3
3.906 ∗ 10−3
5/8
2.667 ∗ 10−3
4.444 ∗ 10−3
7.407 ∗ 10−3
6/8
1.458 ∗ 10−3
4.373 ∗ 10−3
13.12 ∗ 10−3
7/8
0.452 ∗ 10−3
3.186 ∗ 10−3
22.303 ∗ 10−3
1
0
0
37.037 ∗ 10−3The values of the integrals can then be calculated using equation (4), and the solution of the integrals is given below:(10)∫01f1⋅dξ=3.862∗10−3,∫01f2⋅dξ=2.755∗10−3,∫01f3⋅dξ=8.49∗10−3.Thus, the force method equation in numbers and their solutions:(11)lEI0.2471−0.176−0.1760.5433δij.X1X2Xi=θ1θ2θi,X1X2Xi=EIl5.26081.70421.70422.3927rij.θ1θ2θi.The stiffness matrixrij represents the stiffness matrix for a rod or beam with variable cross-section stiffness EIx. Figure 2 shows the bending moment diagrams corresponding to rod ends unit rotation. As shown in this figure, the values of bending moments differ significantly compared to the traditional method for constant stiffness.Figure 2
Bending moment diagrams due to unit rotation at ends.


#### 3.1.2. Finite Element Model and Verification for the Analytical Model

This section presents a case study which demonstrates the usage of the proposed model to solve a statically indeterminate beam with variable cross section. The beam has a total length of 4 m and fixed end moments at both ends; the cross section at the first end is 0.25 × 0.5 m rectangular cross section and 0.25 × 0.3 m at the other end of the beam; and the modulus of elasticity was assumed to be 2.1E + 10 Pa.The model was first simulated numerically using the finite element (FE) general software ABAQUS [18]. The beam was simulated using solid element C3D8R, which is an 8-node linear brick with 3 translational degrees of freedom at each node [18]. The model meshing is shown in Figure 3.Figure 3
FE model developed by ABAQUS program.The meshing was used with maximum aspect ratio of 2; boundary condition at the beams end was applied to two reference points at beams ends. Each reference point is coupled with the beams end surface. The loads applied were a rotation of 1 rad at each end. The rotation was separately applied for each end in two different steps. The results and the staining actions produced from the rotation at each end is separated from the result of rotation at the other end (they are not accumulated or added together). The results of the ABAQUS model is indicated in Table2. It should be mentioned that the larger section is referred as End 1, the smaller section is referred as End 2, and the member end moment calculated by the FE model is referred as MFE.Table 2
Comparison between the FE model results and the proposed model results.
Deformation type
End number
M
FE (N·m)
M
model (n=8) (N·m)
M
FE/Mmodel
Rotation at End 1 = 1 rad
End 1
36649000
36860593.34
0.994
End 2
12222900
11962246.09
1.022
Rotation at End 2 = 1 rad
End 1
12222900
11962246.09
1.022
End 2
17063300
16765176.8
1.018A sensitive analysis has been done in order to investigate the effect of the number of division on the accuracy of the proposed model for different beam height ratio at both its ends. Figure4 shows an example for the sensitive analysis in the case of height ratio at different ends being 0.6. In Figure 4, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of division used in the analysis, and the number of division used in this sensitive analysis is 2, 4, 6, 8, 1, and 12. It should be mentioned that the proposed model was programmed using the MATLAB program [19], which makes it easier to use the proposed model for different numbers of divisions. It can be observed the use of number of nine divisions gives a good result that is in reasonable agreement with the FE model results.Figure 4
Sensitivity analysis for number of division effect on the result of the end moments. End moment at (a) End 1 due to rotation at End 1 and (b) End 2 due to rotation at End 2.
(a)
(b)Table2 shows a comparison between the result obtained by the FE model and that obtained by the proposed model. In Table 1, MFE refers to the end moment obtained by the FE model at both Ends 1 and 2 due to rotations = 1 rad at both ends separately, where Mmodel refers to the end moment obtained by the proposed model using 8 divisions at Ends 1 and 2. As can be seen from the table, there is a very good agreement between the FE and the proposed model, which shows that the proposed model can accurately predict the end moment of indeterminate beams with variable cross section.Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between beam depths at both ends (h2/h1). The obtained number of divisions will give the best result that is in good agreement with the FE model; the proposed equation is shown in the following equation. The R-squared value (R2) for this equation is 0.99, which shows the accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 5.(12)n=6.2098h2h1−0.875,where n is the required number of integral divisions, h1 is the height of the beam at the first end, and h2 is the height of the bema at the second end where h1 ≥ h2.Figure 5
Proposed equation for prediction the required number of integral divisions.


### 3.2. Application 2: Curved Element with Constant Inertia

#### 3.2.1. Analytical Model for Curved Element with Constant Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the circular arch due to clockwise unit rotation at both ends of the arch separately and a horizontal displacement at end (B) as shown in Figure6.Figure 6
The statical system of the circular arch and the proposed main system.Illustrative diagram for the solution using the proposed model is presented in Figure7. Step-by-step solution using the proposed model is indicated in equations (13)–(18).Figure 7
Bending moment diagram due to redundant.Givenzi=1 (Z1 and Z2 are the rotation at both ends, respectively, and Z3 is the horizontal displacement at End B), EI, r=6m, f=2m, yc=r−f=4m, cosφ0=yc/r=2/3, φ0=48.19°=0.8411rad, l=2r⋅φ0=10.09m, s=r⋅φs, φ=2φ0, Mx1=1−s/l, Mx2=−s/l, s/l=ξ=φs/φ, Mx3=1·ys=z⋅cosφ0−φs−cosφ0.(13)δ11=∫011−ξ2EI⋅l⋅dξ=lEI∫011−2ξ+ξ2⋅dξ=lEI⋅13,(14)δ12=−∫011−ξ⋅ξEI⋅l⋅dξ=−lEI∫01ξ−ξ2⋅dξ=−lEI⋅16,(15)δ13=∫011−ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=l⋅rEI∫001−ξcosφ0−φs−cosφ0⋅dξ=l⋅rEI∫01f1⋅dξ,(16)δ22=∫01ξ2EI⋅l⋅dξ=lEI⋅13,(17)δ23=−∫01ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=−l⋅rEI∫01ξcosφ0−φs−cosφ0⋅dξ=−l⋅rEI∫01f2⋅dξ,(18)δ33=∫01r2⋅cosφ0−φs−cosφ02EI⋅l⋅dξ=l⋅r2EI∫01cosφ0−φs−cosφ02⋅dξ=l⋅r2EI∫01f3⋅dξ.For the numerical integration of functions, the interval is divided into 8 pieces. Function valuesfi are given in Table 3.Table 3
Function valuesfi at different sections on the arch.
ξ
f
1
f
2
f
3
0
0
0
0
1/8
0.1233
0.0176
0.0198
2/8
0.1846
0.0615
0.0606
3/8
0.1946
0.1167
0.0969
4/8
0.1667
0.1667
0.1111
5/8
0.1167
0.1946
0.0969
6/8
0.0615
0.1846
0.0606
7/8
0.0176
0.1233
0.0198
1
0
0
0The numerical integration was solved using the method of averages:(19)ω=yH/2+∑i=17yi+yk/28.The values of the integrals:(20)∫01f1⋅dξ=0.10813,∫01f2⋅dξ=0.10813,∫01f3⋅dξ=0.05823.The equations of the force method in the numbers:(21)lEI13−160.10813r−1613−0.10813r0.10813r−0.10813r0.05823r2·X1X2X3=100010001·z1z2z3.Their solution provides(22)100010001·X1X2X3Xj=EIl8.21313−2.21313−19.48178r−2.213138.2131319.48178r−19.48178r19.48178r90.085r2rij.z1z2z3zj.To construct a single diagramMzj, it is necessary to fix the unit diagrams Mxj and put them together:(23)Mzj=∑Mxi⋅Xzj.DiagramsMzj are shown in Figure 8. Ordinates of the individual diagrams Mxi, corrected diagrams Mxi⋅Xzj, and individual diagrams Mzj are shown in Table 4.Figure 8
Schematic diagram for the proposed model solving arch with constant section and the corresponding bending moment diagram.Table 4
Function valuesMxi at different sections on the arch.
ξ
M
x
1
M
x
2
M
x
3
M
x
1
⋅
X
1
z
1
M
x
2
⋅
X
2
z
1
M
x
3
⋅
X
3
z
1
M
z
1
M
x
1
⋅
X
1
z
2
M
x
2
⋅
X
2
z
2
M
x
3
⋅
X
3
z
2
M
z
2
M
x
1
⋅
X
1
z
3
M
x
2
⋅
X
2
z
3
M
x
3
⋅
X
3
z
3
M
z
3
0
1
0
0
8.2131
0
0
8.2131
−2.2131
0
0
−2.2131
−
19.4818
/
r
0
0
−3.2470
1/8
7/8
−1/8
0.4227
f
7.1865
0.2766
−
8.2347
f
/
r
4.7182
−1.9365
−1.0266
8.2347
f
/
r
−0.2182
−
17.0466
/
r
−
2.4352
/
r
38.078f/r2
−1.1315
2/8
6/8
−2/8
0.7387
f
6.1598
0.5533
−
14.3904
f
/
r
1.9163
−1.6598
−2.0533
14.3904
f
/
r
1.0837
−
14.6113
/
r
−
4.8704
/
r
66.5422f/r2
0.4493
3/8
5/8
−3/8
0.9340
f
5.1382
0.8299
−
18.1958
f
/
r
−1.0213
−1.3832
−3.0799
18.1958
f
/
r
1.6021
−
12.1761
/
r
−
7.3057
/
r
84.1385f/r2
1.4274
4/8
4/8
−4/8
f
4.1066
1.1066
−
19.4818
f
/
r
−1.2808
−1.1066
−4.1066
19.4818
f
/
r
1.2808
−
9.7409
/
r
−
9.7409
/
r
90.085f/r2
1.7578
5/8
3/8
−5/8
0.9340
f
3.0799
1.3832
−
18.1958
f
/
r
−1.6021
−0.8299
−5.1332
18.1958
f
/
r
1.0213
−
7.3057
/
r
−
12.1761
/
r
84.1385f/r2
1.4274
6/8
2/8
−6/8
0.7387
f
2.0533
1.6598
−
14.3904
f
/
r
−1.0837
−0.5533
−6.1598
14.3904
f
/
r
−1.9163
−
4.8704
/
r
−
14.6113
/
r
66.5422f/r2
0.4493
7/8
1/8
−7/8
0.4227
f
1.0266
1.9365
−
8.2347
f
/
r
0.2182
−0.2766
−7.1865
8.2347
f
/
r
−4.7182
−
2.4352
/
r
−
17.0466
/
r
38.078f/r2
1.1315
1
0
−1
0
0
2.2131
0
2.2131
0
−8.2131
0
−8.2131
0
−
19.4818
/
r
0
−3.2470


#### 3.2.2. Finite Element Model and Verification for the Circular Arch Analytical Model

This section shows a case study, which verifies the accuracy of the proposed model to solve a statically indeterminate curved beam (arch) with constant cross section. The arch has a radius of 6 m and a height of 4 m with fixed end moments at both ends, the cross section of the arch was 0.25 × 0.4 m, and the modulus of elasticity was assumed to be 2.1E + 10 Pa. An outline of the arch is shown in Figure 9.Figure 9
FE model developed by ABAQUS program for circular arch.Similar to the previous case, the model is first numerically simulated using the finite element (FE) general software ABAQUS [18]. The arch is simulated using solid element C3D8R.A sensitive analysis has been done in order to investigate the effect of the number of divisions on the accuracy of the proposed model. The result of the sensitive analysis is presented in Figure10. In Figure 10, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of divisions used in the analysis; the number of divisions used in this sensitive analysis is 4 to 16. Similar to the first case, the proposed model was programmed using the MATLAB program [19]. It can be observed the use of number of divisions equal to twelve gives result that is in reasonable agreement with the FE model results. At the number of 12 divisions in Figure 10, the ratio of MFE/Mmodel is equal to 0.96, 0.98, and 1.14 at End 1, middle, and End 2 of the arch, respectively, due to rotation at End 1. Due to horizontal displacement at End 2, the ratio of MFE/Mmodel is equal to 1 and 0.95 at End 1 and middle of the arch, respectively. It should be noted that the solution is done faster using the developed program compared to the FE model. For instance, for structural element with variable section, the time needed to build and execute the FE model is 10 times that required to input parameters and get solution from the developed numerical model.Figure 10
Sensitive analysis for number of division effect on internal forces of a circular arch. (a) Internal forces due to rotation at End 1 and (b) internal forces due to horizontal displacement at End 2.
(a)
(b)Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between arch height and its radius (f/r). The obtained number of division has good agreement with the FE model, and the proposed equation is shown in equation (16). The R-squared value (R2) for this equation is 0.98, which shows the good accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 11.(24)n=55.01fr2−13.22fr+8.55.Figure 11
Proposed equation for predicting the required number of integral divisions.


## 3.1. Application 1: Beam with Variable Inertia

### 3.1.1. Analytical Model for Beams with Variable Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the beam due to clockwise unit rotation separately at both ends of the beam as shown in Figure1.Figure 1
The proposed model for solving beam with variable section.Givenθi=1.0 (rotation at each end of the beam separately),(6)hx=0=1.25h,hx=L/2=h,hx=L=0.75h,hx=h5−2x/L4,b=constant,where h is the depth of the beam at the midspan, hx is the depth of the beam at distance x from the left side of the beam, and b is the beam width.(7)Ix=bhx312,Ix=L/2=I=bh312,Ix=I5−2x/l364.The equations of the force method:(8)δ11δ12δ21δ22X1X2=θ1θ2,δij=∫01Mx1⋅Mx2EIx⋅dx,Mx1=1−xl,Mx2=xl,xl=ξ,Mξ1=1−ξ,Mξ2=ξ.By transformation fromx to ξ, δij equations take the form presented in the following equation:(9)δ11=∫011−ξ2⋅64lEI5−2ξ3dξ=64lEI∫01f1⋅dξ,f1=1−ξ25−2ξ3,δ12=−∫011−ξ⋅ξ⋅64lEI5−2ξ3dξ=−64lEI∫01f2⋅dξ,f2=1−ξ⋅ξ5−2ξ3,δ22=∫01ξ2⋅64lEI5−2ξ3dξ=−64lEI∫01f3⋅dξ,f3=ξ25−2ξ3.For the numerical integration of functionZi, assume the interval ξ is divided into n parts. The values of the function Zi are given in Table 1 using n=8 parts.Table 1
Values of the functionZi using n=8 parts.
ξ
f
1
f
2
f
3
0
8 ∗ 10−3
0
0
1/8
7.144 ∗ 10−3
1.021 ∗ 10−3
0.146 ∗ 10−3
2/8
6.179 ∗ 10−3
2.058 ∗ 10−3
0.686 ∗ 10−3
3/8
5.088 ∗ 10−3
3.053 ∗ 10−3
1.832 ∗ 10−3
4/8
3.906 ∗ 10−3
3.906 ∗ 10−3
3.906 ∗ 10−3
5/8
2.667 ∗ 10−3
4.444 ∗ 10−3
7.407 ∗ 10−3
6/8
1.458 ∗ 10−3
4.373 ∗ 10−3
13.12 ∗ 10−3
7/8
0.452 ∗ 10−3
3.186 ∗ 10−3
22.303 ∗ 10−3
1
0
0
37.037 ∗ 10−3The values of the integrals can then be calculated using equation (4), and the solution of the integrals is given below:(10)∫01f1⋅dξ=3.862∗10−3,∫01f2⋅dξ=2.755∗10−3,∫01f3⋅dξ=8.49∗10−3.Thus, the force method equation in numbers and their solutions:(11)lEI0.2471−0.176−0.1760.5433δij.X1X2Xi=θ1θ2θi,X1X2Xi=EIl5.26081.70421.70422.3927rij.θ1θ2θi.The stiffness matrixrij represents the stiffness matrix for a rod or beam with variable cross-section stiffness EIx. Figure 2 shows the bending moment diagrams corresponding to rod ends unit rotation. As shown in this figure, the values of bending moments differ significantly compared to the traditional method for constant stiffness.Figure 2
Bending moment diagrams due to unit rotation at ends.


### 3.1.2. Finite Element Model and Verification for the Analytical Model

This section presents a case study which demonstrates the usage of the proposed model to solve a statically indeterminate beam with variable cross section. The beam has a total length of 4 m and fixed end moments at both ends; the cross section at the first end is 0.25 × 0.5 m rectangular cross section and 0.25 × 0.3 m at the other end of the beam; and the modulus of elasticity was assumed to be 2.1E + 10 Pa.The model was first simulated numerically using the finite element (FE) general software ABAQUS [18]. The beam was simulated using solid element C3D8R, which is an 8-node linear brick with 3 translational degrees of freedom at each node [18]. The model meshing is shown in Figure 3.Figure 3
FE model developed by ABAQUS program.The meshing was used with maximum aspect ratio of 2; boundary condition at the beams end was applied to two reference points at beams ends. Each reference point is coupled with the beams end surface. The loads applied were a rotation of 1 rad at each end. The rotation was separately applied for each end in two different steps. The results and the staining actions produced from the rotation at each end is separated from the result of rotation at the other end (they are not accumulated or added together). The results of the ABAQUS model is indicated in Table2. It should be mentioned that the larger section is referred as End 1, the smaller section is referred as End 2, and the member end moment calculated by the FE model is referred as MFE.Table 2
Comparison between the FE model results and the proposed model results.
Deformation type
End number
M
FE (N·m)
M
model (n=8) (N·m)
M
FE/Mmodel
Rotation at End 1 = 1 rad
End 1
36649000
36860593.34
0.994
End 2
12222900
11962246.09
1.022
Rotation at End 2 = 1 rad
End 1
12222900
11962246.09
1.022
End 2
17063300
16765176.8
1.018A sensitive analysis has been done in order to investigate the effect of the number of division on the accuracy of the proposed model for different beam height ratio at both its ends. Figure4 shows an example for the sensitive analysis in the case of height ratio at different ends being 0.6. In Figure 4, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of division used in the analysis, and the number of division used in this sensitive analysis is 2, 4, 6, 8, 1, and 12. It should be mentioned that the proposed model was programmed using the MATLAB program [19], which makes it easier to use the proposed model for different numbers of divisions. It can be observed the use of number of nine divisions gives a good result that is in reasonable agreement with the FE model results.Figure 4
Sensitivity analysis for number of division effect on the result of the end moments. End moment at (a) End 1 due to rotation at End 1 and (b) End 2 due to rotation at End 2.
(a)
(b)Table2 shows a comparison between the result obtained by the FE model and that obtained by the proposed model. In Table 1, MFE refers to the end moment obtained by the FE model at both Ends 1 and 2 due to rotations = 1 rad at both ends separately, where Mmodel refers to the end moment obtained by the proposed model using 8 divisions at Ends 1 and 2. As can be seen from the table, there is a very good agreement between the FE and the proposed model, which shows that the proposed model can accurately predict the end moment of indeterminate beams with variable cross section.Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between beam depths at both ends (h2/h1). The obtained number of divisions will give the best result that is in good agreement with the FE model; the proposed equation is shown in the following equation. The R-squared value (R2) for this equation is 0.99, which shows the accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 5.(12)n=6.2098h2h1−0.875,where n is the required number of integral divisions, h1 is the height of the beam at the first end, and h2 is the height of the bema at the second end where h1 ≥ h2.Figure 5
Proposed equation for prediction the required number of integral divisions.


## 3.1.1. Analytical Model for Beams with Variable Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the beam due to clockwise unit rotation separately at both ends of the beam as shown in Figure1.Figure 1
The proposed model for solving beam with variable section.Givenθi=1.0 (rotation at each end of the beam separately),(6)hx=0=1.25h,hx=L/2=h,hx=L=0.75h,hx=h5−2x/L4,b=constant,where h is the depth of the beam at the midspan, hx is the depth of the beam at distance x from the left side of the beam, and b is the beam width.(7)Ix=bhx312,Ix=L/2=I=bh312,Ix=I5−2x/l364.The equations of the force method:(8)δ11δ12δ21δ22X1X2=θ1θ2,δij=∫01Mx1⋅Mx2EIx⋅dx,Mx1=1−xl,Mx2=xl,xl=ξ,Mξ1=1−ξ,Mξ2=ξ.By transformation fromx to ξ, δij equations take the form presented in the following equation:(9)δ11=∫011−ξ2⋅64lEI5−2ξ3dξ=64lEI∫01f1⋅dξ,f1=1−ξ25−2ξ3,δ12=−∫011−ξ⋅ξ⋅64lEI5−2ξ3dξ=−64lEI∫01f2⋅dξ,f2=1−ξ⋅ξ5−2ξ3,δ22=∫01ξ2⋅64lEI5−2ξ3dξ=−64lEI∫01f3⋅dξ,f3=ξ25−2ξ3.For the numerical integration of functionZi, assume the interval ξ is divided into n parts. The values of the function Zi are given in Table 1 using n=8 parts.Table 1
Values of the functionZi using n=8 parts.
ξ
f
1
f
2
f
3
0
8 ∗ 10−3
0
0
1/8
7.144 ∗ 10−3
1.021 ∗ 10−3
0.146 ∗ 10−3
2/8
6.179 ∗ 10−3
2.058 ∗ 10−3
0.686 ∗ 10−3
3/8
5.088 ∗ 10−3
3.053 ∗ 10−3
1.832 ∗ 10−3
4/8
3.906 ∗ 10−3
3.906 ∗ 10−3
3.906 ∗ 10−3
5/8
2.667 ∗ 10−3
4.444 ∗ 10−3
7.407 ∗ 10−3
6/8
1.458 ∗ 10−3
4.373 ∗ 10−3
13.12 ∗ 10−3
7/8
0.452 ∗ 10−3
3.186 ∗ 10−3
22.303 ∗ 10−3
1
0
0
37.037 ∗ 10−3The values of the integrals can then be calculated using equation (4), and the solution of the integrals is given below:(10)∫01f1⋅dξ=3.862∗10−3,∫01f2⋅dξ=2.755∗10−3,∫01f3⋅dξ=8.49∗10−3.Thus, the force method equation in numbers and their solutions:(11)lEI0.2471−0.176−0.1760.5433δij.X1X2Xi=θ1θ2θi,X1X2Xi=EIl5.26081.70421.70422.3927rij.θ1θ2θi.The stiffness matrixrij represents the stiffness matrix for a rod or beam with variable cross-section stiffness EIx. Figure 2 shows the bending moment diagrams corresponding to rod ends unit rotation. As shown in this figure, the values of bending moments differ significantly compared to the traditional method for constant stiffness.Figure 2
Bending moment diagrams due to unit rotation at ends.


## 3.1.2. Finite Element Model and Verification for the Analytical Model

This section presents a case study which demonstrates the usage of the proposed model to solve a statically indeterminate beam with variable cross section. The beam has a total length of 4 m and fixed end moments at both ends; the cross section at the first end is 0.25 × 0.5 m rectangular cross section and 0.25 × 0.3 m at the other end of the beam; and the modulus of elasticity was assumed to be 2.1E + 10 Pa.The model was first simulated numerically using the finite element (FE) general software ABAQUS [18]. The beam was simulated using solid element C3D8R, which is an 8-node linear brick with 3 translational degrees of freedom at each node [18]. The model meshing is shown in Figure 3.Figure 3
FE model developed by ABAQUS program.The meshing was used with maximum aspect ratio of 2; boundary condition at the beams end was applied to two reference points at beams ends. Each reference point is coupled with the beams end surface. The loads applied were a rotation of 1 rad at each end. The rotation was separately applied for each end in two different steps. The results and the staining actions produced from the rotation at each end is separated from the result of rotation at the other end (they are not accumulated or added together). The results of the ABAQUS model is indicated in Table2. It should be mentioned that the larger section is referred as End 1, the smaller section is referred as End 2, and the member end moment calculated by the FE model is referred as MFE.Table 2
Comparison between the FE model results and the proposed model results.
Deformation type
End number
M
FE (N·m)
M
model (n=8) (N·m)
M
FE/Mmodel
Rotation at End 1 = 1 rad
End 1
36649000
36860593.34
0.994
End 2
12222900
11962246.09
1.022
Rotation at End 2 = 1 rad
End 1
12222900
11962246.09
1.022
End 2
17063300
16765176.8
1.018A sensitive analysis has been done in order to investigate the effect of the number of division on the accuracy of the proposed model for different beam height ratio at both its ends. Figure4 shows an example for the sensitive analysis in the case of height ratio at different ends being 0.6. In Figure 4, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of division used in the analysis, and the number of division used in this sensitive analysis is 2, 4, 6, 8, 1, and 12. It should be mentioned that the proposed model was programmed using the MATLAB program [19], which makes it easier to use the proposed model for different numbers of divisions. It can be observed the use of number of nine divisions gives a good result that is in reasonable agreement with the FE model results.Figure 4
Sensitivity analysis for number of division effect on the result of the end moments. End moment at (a) End 1 due to rotation at End 1 and (b) End 2 due to rotation at End 2.
(a)
(b)Table2 shows a comparison between the result obtained by the FE model and that obtained by the proposed model. In Table 1, MFE refers to the end moment obtained by the FE model at both Ends 1 and 2 due to rotations = 1 rad at both ends separately, where Mmodel refers to the end moment obtained by the proposed model using 8 divisions at Ends 1 and 2. As can be seen from the table, there is a very good agreement between the FE and the proposed model, which shows that the proposed model can accurately predict the end moment of indeterminate beams with variable cross section.Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between beam depths at both ends (h2/h1). The obtained number of divisions will give the best result that is in good agreement with the FE model; the proposed equation is shown in the following equation. The R-squared value (R2) for this equation is 0.99, which shows the accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 5.(12)n=6.2098h2h1−0.875,where n is the required number of integral divisions, h1 is the height of the beam at the first end, and h2 is the height of the bema at the second end where h1 ≥ h2.Figure 5
Proposed equation for prediction the required number of integral divisions.


## 3.2. Application 2: Curved Element with Constant Inertia

### 3.2.1. Analytical Model for Curved Element with Constant Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the circular arch due to clockwise unit rotation at both ends of the arch separately and a horizontal displacement at end (B) as shown in Figure6.Figure 6
The statical system of the circular arch and the proposed main system.Illustrative diagram for the solution using the proposed model is presented in Figure7. Step-by-step solution using the proposed model is indicated in equations (13)–(18).Figure 7
Bending moment diagram due to redundant.Givenzi=1 (Z1 and Z2 are the rotation at both ends, respectively, and Z3 is the horizontal displacement at End B), EI, r=6m, f=2m, yc=r−f=4m, cosφ0=yc/r=2/3, φ0=48.19°=0.8411rad, l=2r⋅φ0=10.09m, s=r⋅φs, φ=2φ0, Mx1=1−s/l, Mx2=−s/l, s/l=ξ=φs/φ, Mx3=1·ys=z⋅cosφ0−φs−cosφ0.(13)δ11=∫011−ξ2EI⋅l⋅dξ=lEI∫011−2ξ+ξ2⋅dξ=lEI⋅13,(14)δ12=−∫011−ξ⋅ξEI⋅l⋅dξ=−lEI∫01ξ−ξ2⋅dξ=−lEI⋅16,(15)δ13=∫011−ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=l⋅rEI∫001−ξcosφ0−φs−cosφ0⋅dξ=l⋅rEI∫01f1⋅dξ,(16)δ22=∫01ξ2EI⋅l⋅dξ=lEI⋅13,(17)δ23=−∫01ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=−l⋅rEI∫01ξcosφ0−φs−cosφ0⋅dξ=−l⋅rEI∫01f2⋅dξ,(18)δ33=∫01r2⋅cosφ0−φs−cosφ02EI⋅l⋅dξ=l⋅r2EI∫01cosφ0−φs−cosφ02⋅dξ=l⋅r2EI∫01f3⋅dξ.For the numerical integration of functions, the interval is divided into 8 pieces. Function valuesfi are given in Table 3.Table 3
Function valuesfi at different sections on the arch.
ξ
f
1
f
2
f
3
0
0
0
0
1/8
0.1233
0.0176
0.0198
2/8
0.1846
0.0615
0.0606
3/8
0.1946
0.1167
0.0969
4/8
0.1667
0.1667
0.1111
5/8
0.1167
0.1946
0.0969
6/8
0.0615
0.1846
0.0606
7/8
0.0176
0.1233
0.0198
1
0
0
0The numerical integration was solved using the method of averages:(19)ω=yH/2+∑i=17yi+yk/28.The values of the integrals:(20)∫01f1⋅dξ=0.10813,∫01f2⋅dξ=0.10813,∫01f3⋅dξ=0.05823.The equations of the force method in the numbers:(21)lEI13−160.10813r−1613−0.10813r0.10813r−0.10813r0.05823r2·X1X2X3=100010001·z1z2z3.Their solution provides(22)100010001·X1X2X3Xj=EIl8.21313−2.21313−19.48178r−2.213138.2131319.48178r−19.48178r19.48178r90.085r2rij.z1z2z3zj.To construct a single diagramMzj, it is necessary to fix the unit diagrams Mxj and put them together:(23)Mzj=∑Mxi⋅Xzj.DiagramsMzj are shown in Figure 8. Ordinates of the individual diagrams Mxi, corrected diagrams Mxi⋅Xzj, and individual diagrams Mzj are shown in Table 4.Figure 8
Schematic diagram for the proposed model solving arch with constant section and the corresponding bending moment diagram.Table 4
Function valuesMxi at different sections on the arch.
ξ
M
x
1
M
x
2
M
x
3
M
x
1
⋅
X
1
z
1
M
x
2
⋅
X
2
z
1
M
x
3
⋅
X
3
z
1
M
z
1
M
x
1
⋅
X
1
z
2
M
x
2
⋅
X
2
z
2
M
x
3
⋅
X
3
z
2
M
z
2
M
x
1
⋅
X
1
z
3
M
x
2
⋅
X
2
z
3
M
x
3
⋅
X
3
z
3
M
z
3
0
1
0
0
8.2131
0
0
8.2131
−2.2131
0
0
−2.2131
−
19.4818
/
r
0
0
−3.2470
1/8
7/8
−1/8
0.4227
f
7.1865
0.2766
−
8.2347
f
/
r
4.7182
−1.9365
−1.0266
8.2347
f
/
r
−0.2182
−
17.0466
/
r
−
2.4352
/
r
38.078f/r2
−1.1315
2/8
6/8
−2/8
0.7387
f
6.1598
0.5533
−
14.3904
f
/
r
1.9163
−1.6598
−2.0533
14.3904
f
/
r
1.0837
−
14.6113
/
r
−
4.8704
/
r
66.5422f/r2
0.4493
3/8
5/8
−3/8
0.9340
f
5.1382
0.8299
−
18.1958
f
/
r
−1.0213
−1.3832
−3.0799
18.1958
f
/
r
1.6021
−
12.1761
/
r
−
7.3057
/
r
84.1385f/r2
1.4274
4/8
4/8
−4/8
f
4.1066
1.1066
−
19.4818
f
/
r
−1.2808
−1.1066
−4.1066
19.4818
f
/
r
1.2808
−
9.7409
/
r
−
9.7409
/
r
90.085f/r2
1.7578
5/8
3/8
−5/8
0.9340
f
3.0799
1.3832
−
18.1958
f
/
r
−1.6021
−0.8299
−5.1332
18.1958
f
/
r
1.0213
−
7.3057
/
r
−
12.1761
/
r
84.1385f/r2
1.4274
6/8
2/8
−6/8
0.7387
f
2.0533
1.6598
−
14.3904
f
/
r
−1.0837
−0.5533
−6.1598
14.3904
f
/
r
−1.9163
−
4.8704
/
r
−
14.6113
/
r
66.5422f/r2
0.4493
7/8
1/8
−7/8
0.4227
f
1.0266
1.9365
−
8.2347
f
/
r
0.2182
−0.2766
−7.1865
8.2347
f
/
r
−4.7182
−
2.4352
/
r
−
17.0466
/
r
38.078f/r2
1.1315
1
0
−1
0
0
2.2131
0
2.2131
0
−8.2131
0
−8.2131
0
−
19.4818
/
r
0
−3.2470


### 3.2.2. Finite Element Model and Verification for the Circular Arch Analytical Model

This section shows a case study, which verifies the accuracy of the proposed model to solve a statically indeterminate curved beam (arch) with constant cross section. The arch has a radius of 6 m and a height of 4 m with fixed end moments at both ends, the cross section of the arch was 0.25 × 0.4 m, and the modulus of elasticity was assumed to be 2.1E + 10 Pa. An outline of the arch is shown in Figure 9.Figure 9
FE model developed by ABAQUS program for circular arch.Similar to the previous case, the model is first numerically simulated using the finite element (FE) general software ABAQUS [18]. The arch is simulated using solid element C3D8R.A sensitive analysis has been done in order to investigate the effect of the number of divisions on the accuracy of the proposed model. The result of the sensitive analysis is presented in Figure10. In Figure 10, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of divisions used in the analysis; the number of divisions used in this sensitive analysis is 4 to 16. Similar to the first case, the proposed model was programmed using the MATLAB program [19]. It can be observed the use of number of divisions equal to twelve gives result that is in reasonable agreement with the FE model results. At the number of 12 divisions in Figure 10, the ratio of MFE/Mmodel is equal to 0.96, 0.98, and 1.14 at End 1, middle, and End 2 of the arch, respectively, due to rotation at End 1. Due to horizontal displacement at End 2, the ratio of MFE/Mmodel is equal to 1 and 0.95 at End 1 and middle of the arch, respectively. It should be noted that the solution is done faster using the developed program compared to the FE model. For instance, for structural element with variable section, the time needed to build and execute the FE model is 10 times that required to input parameters and get solution from the developed numerical model.Figure 10
Sensitive analysis for number of division effect on internal forces of a circular arch. (a) Internal forces due to rotation at End 1 and (b) internal forces due to horizontal displacement at End 2.
(a)
(b)Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between arch height and its radius (f/r). The obtained number of division has good agreement with the FE model, and the proposed equation is shown in equation (16). The R-squared value (R2) for this equation is 0.98, which shows the good accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 11.(24)n=55.01fr2−13.22fr+8.55.Figure 11
Proposed equation for predicting the required number of integral divisions.


## 3.2.1. Analytical Model for Curved Element with Constant Inertia

This section describes the proposed analytical model through an application in which the internal moments at supports are determined for the circular arch due to clockwise unit rotation at both ends of the arch separately and a horizontal displacement at end (B) as shown in Figure6.Figure 6
The statical system of the circular arch and the proposed main system.Illustrative diagram for the solution using the proposed model is presented in Figure7. Step-by-step solution using the proposed model is indicated in equations (13)–(18).Figure 7
Bending moment diagram due to redundant.Givenzi=1 (Z1 and Z2 are the rotation at both ends, respectively, and Z3 is the horizontal displacement at End B), EI, r=6m, f=2m, yc=r−f=4m, cosφ0=yc/r=2/3, φ0=48.19°=0.8411rad, l=2r⋅φ0=10.09m, s=r⋅φs, φ=2φ0, Mx1=1−s/l, Mx2=−s/l, s/l=ξ=φs/φ, Mx3=1·ys=z⋅cosφ0−φs−cosφ0.(13)δ11=∫011−ξ2EI⋅l⋅dξ=lEI∫011−2ξ+ξ2⋅dξ=lEI⋅13,(14)δ12=−∫011−ξ⋅ξEI⋅l⋅dξ=−lEI∫01ξ−ξ2⋅dξ=−lEI⋅16,(15)δ13=∫011−ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=l⋅rEI∫001−ξcosφ0−φs−cosφ0⋅dξ=l⋅rEI∫01f1⋅dξ,(16)δ22=∫01ξ2EI⋅l⋅dξ=lEI⋅13,(17)δ23=−∫01ξ⋅rcosφ0−φs−cosφ0EI⋅l⋅dξ=−l⋅rEI∫01ξcosφ0−φs−cosφ0⋅dξ=−l⋅rEI∫01f2⋅dξ,(18)δ33=∫01r2⋅cosφ0−φs−cosφ02EI⋅l⋅dξ=l⋅r2EI∫01cosφ0−φs−cosφ02⋅dξ=l⋅r2EI∫01f3⋅dξ.For the numerical integration of functions, the interval is divided into 8 pieces. Function valuesfi are given in Table 3.Table 3
Function valuesfi at different sections on the arch.
ξ
f
1
f
2
f
3
0
0
0
0
1/8
0.1233
0.0176
0.0198
2/8
0.1846
0.0615
0.0606
3/8
0.1946
0.1167
0.0969
4/8
0.1667
0.1667
0.1111
5/8
0.1167
0.1946
0.0969
6/8
0.0615
0.1846
0.0606
7/8
0.0176
0.1233
0.0198
1
0
0
0The numerical integration was solved using the method of averages:(19)ω=yH/2+∑i=17yi+yk/28.The values of the integrals:(20)∫01f1⋅dξ=0.10813,∫01f2⋅dξ=0.10813,∫01f3⋅dξ=0.05823.The equations of the force method in the numbers:(21)lEI13−160.10813r−1613−0.10813r0.10813r−0.10813r0.05823r2·X1X2X3=100010001·z1z2z3.Their solution provides(22)100010001·X1X2X3Xj=EIl8.21313−2.21313−19.48178r−2.213138.2131319.48178r−19.48178r19.48178r90.085r2rij.z1z2z3zj.To construct a single diagramMzj, it is necessary to fix the unit diagrams Mxj and put them together:(23)Mzj=∑Mxi⋅Xzj.DiagramsMzj are shown in Figure 8. Ordinates of the individual diagrams Mxi, corrected diagrams Mxi⋅Xzj, and individual diagrams Mzj are shown in Table 4.Figure 8
Schematic diagram for the proposed model solving arch with constant section and the corresponding bending moment diagram.Table 4
Function valuesMxi at different sections on the arch.
ξ
M
x
1
M
x
2
M
x
3
M
x
1
⋅
X
1
z
1
M
x
2
⋅
X
2
z
1
M
x
3
⋅
X
3
z
1
M
z
1
M
x
1
⋅
X
1
z
2
M
x
2
⋅
X
2
z
2
M
x
3
⋅
X
3
z
2
M
z
2
M
x
1
⋅
X
1
z
3
M
x
2
⋅
X
2
z
3
M
x
3
⋅
X
3
z
3
M
z
3
0
1
0
0
8.2131
0
0
8.2131
−2.2131
0
0
−2.2131
−
19.4818
/
r
0
0
−3.2470
1/8
7/8
−1/8
0.4227
f
7.1865
0.2766
−
8.2347
f
/
r
4.7182
−1.9365
−1.0266
8.2347
f
/
r
−0.2182
−
17.0466
/
r
−
2.4352
/
r
38.078f/r2
−1.1315
2/8
6/8
−2/8
0.7387
f
6.1598
0.5533
−
14.3904
f
/
r
1.9163
−1.6598
−2.0533
14.3904
f
/
r
1.0837
−
14.6113
/
r
−
4.8704
/
r
66.5422f/r2
0.4493
3/8
5/8
−3/8
0.9340
f
5.1382
0.8299
−
18.1958
f
/
r
−1.0213
−1.3832
−3.0799
18.1958
f
/
r
1.6021
−
12.1761
/
r
−
7.3057
/
r
84.1385f/r2
1.4274
4/8
4/8
−4/8
f
4.1066
1.1066
−
19.4818
f
/
r
−1.2808
−1.1066
−4.1066
19.4818
f
/
r
1.2808
−
9.7409
/
r
−
9.7409
/
r
90.085f/r2
1.7578
5/8
3/8
−5/8
0.9340
f
3.0799
1.3832
−
18.1958
f
/
r
−1.6021
−0.8299
−5.1332
18.1958
f
/
r
1.0213
−
7.3057
/
r
−
12.1761
/
r
84.1385f/r2
1.4274
6/8
2/8
−6/8
0.7387
f
2.0533
1.6598
−
14.3904
f
/
r
−1.0837
−0.5533
−6.1598
14.3904
f
/
r
−1.9163
−
4.8704
/
r
−
14.6113
/
r
66.5422f/r2
0.4493
7/8
1/8
−7/8
0.4227
f
1.0266
1.9365
−
8.2347
f
/
r
0.2182
−0.2766
−7.1865
8.2347
f
/
r
−4.7182
−
2.4352
/
r
−
17.0466
/
r
38.078f/r2
1.1315
1
0
−1
0
0
2.2131
0
2.2131
0
−8.2131
0
−8.2131
0
−
19.4818
/
r
0
−3.2470


## 3.2.2. Finite Element Model and Verification for the Circular Arch Analytical Model

This section shows a case study, which verifies the accuracy of the proposed model to solve a statically indeterminate curved beam (arch) with constant cross section. The arch has a radius of 6 m and a height of 4 m with fixed end moments at both ends, the cross section of the arch was 0.25 × 0.4 m, and the modulus of elasticity was assumed to be 2.1E + 10 Pa. An outline of the arch is shown in Figure 9.Figure 9
FE model developed by ABAQUS program for circular arch.Similar to the previous case, the model is first numerically simulated using the finite element (FE) general software ABAQUS [18]. The arch is simulated using solid element C3D8R.A sensitive analysis has been done in order to investigate the effect of the number of divisions on the accuracy of the proposed model. The result of the sensitive analysis is presented in Figure10. In Figure 10, the vertical axis shows the ratio between the end moment obtained by the FE model (MFE) and the end moment obtained by the proposed model (Mmodel). The horizontal axis shows the number of divisions used in the analysis; the number of divisions used in this sensitive analysis is 4 to 16. Similar to the first case, the proposed model was programmed using the MATLAB program [19]. It can be observed the use of number of divisions equal to twelve gives result that is in reasonable agreement with the FE model results. At the number of 12 divisions in Figure 10, the ratio of MFE/Mmodel is equal to 0.96, 0.98, and 1.14 at End 1, middle, and End 2 of the arch, respectively, due to rotation at End 1. Due to horizontal displacement at End 2, the ratio of MFE/Mmodel is equal to 1 and 0.95 at End 1 and middle of the arch, respectively. It should be noted that the solution is done faster using the developed program compared to the FE model. For instance, for structural element with variable section, the time needed to build and execute the FE model is 10 times that required to input parameters and get solution from the developed numerical model.Figure 10
Sensitive analysis for number of division effect on internal forces of a circular arch. (a) Internal forces due to rotation at End 1 and (b) internal forces due to horizontal displacement at End 2.
(a)
(b)Regression analysis was used to propose a simple equation to predict the required number of integral division (n) based on the ratio between arch height and its radius (f/r). The obtained number of division has good agreement with the FE model, and the proposed equation is shown in equation (16). The R-squared value (R2) for this equation is 0.98, which shows the good accuracy of this model to predict the best number of integral divisions. This is illustrated in Figure 11.(24)n=55.01fr2−13.22fr+8.55.Figure 11
Proposed equation for predicting the required number of integral divisions.


## 4. Conclusions

From the litterateur review, it was concluded that the complication of the available design methods is considered the main problem that engineers face in designing structural elements having variable inertia as well as a curved shape; besides, not all engineers are familiar with using FE programs. The objective of this paper is to cover the gap in this area and provide simple mathematical analytical models that would simplify the designing and analysis of structural elements that have variable inertia and a curved shape. The developed program was obtained using the MATLAB software, and the output results from this program were verified using the FE software [18].The following conclusions can be drawn from this paper:(1)
Significant proposed analytical model based on the flexibility method was developed for straight beams with variable inertia and curved beams (circular arch) with constant inertia.(2)
Good agreement has been found between the proposed model and the numerical model which gives more confidence in using the proposed model in structural analysis for beams with variable cross sections and arch with constant cross section.(3)
A new program was developed using the MATLAB software based on the proposed model. The new program can analyze beams with variable sections and arches with constant section consuming less time compared to the ABAQUS program, which makes the proposed model more reliable compared to the other method of analysis. For structural element with variable sections, the time needed to build and execute the FE model is 10 times that required to input parameters and get solution from the developed numerical model.(4)
The number of divisions used in the investigated analytical models play a main role to get best agreement with FE solution.(5)
New mathematical formulas were investigated for straight beams with variable inertia and curved beams with constant inertia to obtain the optimum number of divisions used in the proposed model.(6)
The developed mathematical models for analyzing straight elements with variable inertia and a curved shape with constant inertia could be a helpful method used in teaching structural analysis of such elements for undergraduate students.(7)
In the future research, the proposed model will be modified to be used with arched beam with variable cross sections.



---

*Source: 1602565-2019-04-07.xml*