1 American Academic Scientific Research Journal for Engineering, Technology, and Sciences ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 http://asrjetsjournal.org/ Linear Strain Triangle Mathematics:Stiffness, Stress and Consistent Load Vector Mohamed Abdalla Almheriegh * Associate Professor, Department of civil Engineering,Faculty of Engineering - Tripoli University,P O Box 82677, Tripoli, Libya Email: malmherigh@gmail.com Abstract Typically, linear strain triangle LST is widely used for elastic analysis not only due to its ease but also due to the good results it can lead to over the traditional constant strain triangle CST. In plane elastic analysis of two combined materials the LST is quite reliable and trust worthy in terms of stress results for the same element size, LST is very useful in modelling combined materials thus believed to be efficient and can easily take into consideration, self-weight of used materials, strains due to causes other than loading; such as moisture, temperature, creep and shrinkage are easily incorporated. A mathematical formulation of stiffness matrix, stress and the more rarely dealt with consistent load vector for loads distributed on element edge are proved and highlighted for young engineers who mostly dealing with readily used FE codes, the formula proved for LST element load vector is further extended and thus a consisted load vector for conical axe-symmetric shell element is introduced. Keywords: consistent load vector for conical axe-symmetric shell element or pressure vessel element; modeling of multiple materials; mathematics of finite element LST; stiffness matrix, stress and, consistent load vector. 1. Introduction The constant strain triangle CST provides good introduction that is easily understood and formulated. It is also a useful element for practical problems, but does require that rather small elements be used in regions where stress varies abruptly. This means that the number of such elements tends to be large when reliable stress values are needed. By using more advanced elements such as LST, in which the stress may vary within the element, much better results can be obtained with the same number of elements. ------------------------------------------------------------------------ * Corresponding author. http://asrjetsjournal.org/ American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 2 There is no limit to the degree of variation in stress that can be specified and, in this paper, only one step in this direction will be taken; provision will be made for the stress to vary linearly within the element a mathematical presentation of stiffens matrix, stress and the rarely dealt with consistent load vector is proved based on conservation of energy principle the formula is readily extended to axe-symmetric or pressure vessel element. 2. Displacement Functions Later it will be shown that the strain, and hence the stress varies linearly if the displacement function is of the form u= α1+α2x+ α3y+ α4x 2 + α5xy+ α6y 2 v= α7+α8x+ α9y+ α10x 2 + α11xy+ α12y 2 or, u v = Pα 1 and 1 x y x 2 xy y 2 0 0 0 0 0 0 P= 0 0 0 0 0 0 1 x y x 2 xy y 2 2 To determine the six-unknown coefficients in the function for u we must know six horizontal displacements at six specified points. The same is true for the unknowns in the v function. This can be accomplished by introducing nodes at the midpoints of the element sides as shown in figures (a, b and c) It will be noted that local axis has been established with origin at the centroid of the element. Dealing with these local axes does not alter the values in the stiffness matrix but does make it easier to perform the integration in the stiffness formula K= [ A -1 ] T ∫ B T DB dv A -1 The displacements at node 1, δ1 and δ7 are given as δ1 δ7 = P x = x1 . α 3 y= y1 1 x1 y1 x1 2 x1y1 y1 2 0 0 0 0 0 0 = 0 0 0 0 0 0 1 x1 y1 x1 2 x1y1 y1 2 α Similar equations can be written for δ components at all other nodes. When the δ components are put into numerical order these equations become Ᾱ 0 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 3 δ= 0 Ᾱ α 4 where 1 x1 y1 x1 2 x1y1 y1 2 1 x2 y2 x2 2 x2y2 y2 2 Ᾱ = 1 x3 y3 x3 2 x3y3 y3 2 5 1 x4 y4 x4 2 x4y4 y4 2 1 x5 y5 x5 2 x5y5 y5 2 1 x6 y6 x6 2 x6y6 y6 2 Then equation 4 becomes Ᾱ 0 δ= A α and A = 6 0 Ᾱ Solving for α gives α= A -1 δ 7 A is a12x12 matrix and could be inverted directly, but much time is saved by inverting the submatrix Ᾱ and then incorporate in A -1 The displacement function in 1 is then written u = P A -1 δ 8 v Now compatibility of the displacements is then checked at points located on the common boundary between adjacent elements. That is, we will determine if pairs of points that touched one another across a boundary before displacement, remain in contact after displacement and hence whether displacement causes gaps to appear. Considering only u displacement, the function in general is u= α1+α2x+ α3y+ α4x 2 + α5xy+ α6y 2 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 4 5 4 6 centroid 3 1 2 a- Local axis and Local Numbering System δ11 δ5 δ12 δ10 δ9 δ6 δ4 δ7 δ8 δ3 δ2 δ1 b- Nodal Displacement Components f11 f5 f10 f12 f4 f9 f6 f7 f8 f3 f2 f1 c- Forces Acting on Element However, we are interested in u only at points on a particular boundary. Let the boundary be parallel to y axis which means that all x components are constants and the u function takes the parabolic form ub =A +By +Cy 2 , on this boundary there are 3 nodes, hence 3 unknown values of ub are the constants A, B and C can be determined and shape of the displaced edge would be a uniquely determined parabola. If we now consider the u displacements of points just across the boundary, they too vary parabolically and are uniquely determined by the motion of three nodes. Since displacements are common at the three nodal points, the parabolic displacements on both sides of the boundary are identical. Hence there will be no tendency for gaps open or for overlap to occur. Similar treatment of vb would show that there is no slipping along the interface formed by the boundary. For an oblique boundary the manipulations would be more complex but the conclusions identical. Consequently, displacements are compatible at all element interfaces. American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 5 3. Stiffness of the Linear Strain Triangle In the strain formula Ɛ = Δ P A -1 δ in which the deferential operator matrix 𝜕 𝜕𝑥 0 Δ = 0 𝛿 𝛿𝑦 𝛿 𝛿𝑦 𝜕 𝜕𝑥 The polynomial matrix given earlier in 2 from which B = Δ P becomes B = 𝜕 𝜕𝑥 0 0 𝛿 𝛿𝑦 𝛿 𝛿𝑦 𝜕 𝜕𝑥 1 x y x 2 xy y 2 0 0 0 0 0 0 0 0 0 0 0 0 1 x y x 2 xy y 2 0 1 0 2x y 0 0 0 0 0 0 0 = 0 0 0 0 0 0 0 0 1 0 x 2y 9 0 0 1 0 x 2y 0 1 0 2x y 0 In the strain equation Ɛ = B A -1 δ The only variables occur in matrix B and hence are seen to have unit exponent, consequently strain varies linearly within the element. With matrix D given by 1 𝒗 0 D = 𝐸 1−𝑣 𝒗 1 0 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 6 0 0 1−𝑣 2 All matrices in the stiffness equation K above are known with the incremental volume dv for an element with a uniform thickness, t can be taken as dv =t dA, therefor K= 𝐴-1 T t ʃ B T D B dA A -1 10 Matrix B contains x and y terms as well as constants therefore when product of B T D B is considered, it is evident that we must integrate terms containing constants, x, y x 2 , xy, y 2 . These integrations lead to some intricate formulas except in the case where the origin is at the centroid, then the solutions are quite simple. Even these are difficult to compute when working with Cartesian coordinates; but by means of area coordinates following [1,2], through [7] formulas can be derived. With the use of these formulas, we found that ∫ B T D B dA = 𝑎𝑟𝑒𝑎 𝑥 𝐸 1−v2 x 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 𝒗 0 0 0 0 0 μ 0 0 0 0 μ 0 0 0 0 0 0 0 4ꭓꭓ 2ꭕу 0 0 0 0 0 2𝒗ꭓꭕ 4𝒗ꭓу 0 0 0 2ꭕу μꭓꭕ+уу 2ꭕу 0 0 0 2μꭓꭕ (μ+𝒗)ꭕу 2 𝒗уу 0 0 0 0 2μꭕу 4μуу 0 0 0 4μꭕу 2μуу 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 μ 0 0 0 0 μ 0 0 0 0 11 0 𝒗 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2μꭓꭕ 4μꭕу 0 0 0 4μꭓꭕ 2μꭕу 0 0 0 0 2𝒗ꭓꭕ (μ+𝒗)ꭕу 2μуу 0 0 0 2μꭕу ꭓꭕ+μуу 0 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 7 0 0 0 4𝒗ꭓу 2𝒗уу 0 0 0 0 0 2ꭕу 4уу Where μ = (1 – 𝒗 ) ̸ 2 ꭓꭕ = 1 ̸ 12 ( ꭕ1 2 + ꭕ3 2 + ꭕ5 2 ) ꭕу = 1/ 12 ( ꭕ1у1 + ꭕ3у3 + ꭕ5у5 ) уу = 1/ 12 ( у1 2 + у3 2 + у5 2 ) We are now able to evaluate the stiffness matrix for any element using the elastic constants of the material, the element stiffness and the coordinates of the nodes based on axes having origin at the centroid of the element. With A -1 evaluated, stiffness matrix is done using equations 5, 10 and11. 4. Stress Calculation in LST The stiffness of each element in turn is accumulated to give the stiffness matrix for the assembly of elements. Loads and constraints are applied and solution is run to give all displacements U, taking each element in turn the values in δ are extracted from U and used to calculate stress by σ = DB A -1 δ 12 It is evident that stress is linear within the element or function of x and y coordinated. Compared with constant strain triangle CST and for triangles meshed with same size, when compared with known solution the LST shows more accurate values given that the analyst is guarding against errors [3] and [4]. 5. Consistent Load Vector When concentrated loads occur, they are handled simply by placing a node at the point of load application and taking the load components as known force components. However, a distributed traction acting on an element boundary is more difficult to deal with. In the LST (linear strain triangle), a traction on an element edge could be treated as concentrated loads at the end nodes by making the tow systems statically equivalent. With the LST, there are three nodes on any edge and the principle of static equivalence is not sufficient to determine distribution uniquely. The solution to this problem will be found by resorting to energy principles. Consider an element boundary such as 1-2-3 in the shown figure, if a horizontal displacement occurs, the edge displacement will be parabolic in form and the new edge location may be taken as 1’-2’-3’ as shown by the broken lines. American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 8 The edge displacement is given by a u= ay 2 +by +c = [y 2 y 1] b c when this displacement takes place, the edge load, p, does work on the system. Over an incremental length of edge dy long the work done is a y 3 3’ 2’ 2 y L/ 2 L/ 2 1 1’ x P0 P = P0 + sy u = ay 2 + by + c American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 9 dWp = p u t dy = ( p0 + sy ) [y 2 y 1] b t dy c a = [p0y 2 + sy 3 p0y +sy 2 p0+sy] b t dy c and the work done on the whole edge is a Wp = ∫ dwp = t [p0L 3 /12 SL 3 /12 p0L] b c Consider now a concentrated load system, p1, p2 and p3 which would be used to replace the distributed system. These forces move through distances U1, U2 and U3 which can be written U1 u (y=-L/2) (-L/2) 2 -L/2 1 a U2 = u (y=0) = 0 0 1 b U3 u (y=L/2) (L/2) 2 -L/2 1 c The work done by the concentrated forces is given by U1 Wp = [P1 P2 P3] U2 U3 L 2 /4 -L/2 1 a = [P1 P2 P3] 0 0 1 b L 2 /4 L/2 1 c Making the two loading systems equivalent by equating the work done by the systems gives L 2 /4 -L/2 1 a [P1 P2 P3] 0 0 1 b L 2 /4 L/2 1 c a = t [p0L 3 /12 SL 3 /12 p0L] b c when the right-hand element is removed from both sides and both sides are post multiplied by the inverse of - L/2 L/2 American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 10 L 2 /4 -L/2 1 0 0 1 L 2 /4 L/2 1 Which is 2/L 2 -4/L 2 2/L 2 -1/L 0 1/L 0 1 0 We get 2/L 2 -4/L 2 2/L 2 [P1 P2 P3] = t [p0L 3 /12 SL 3 /12 p0L] -1/L 0 1/L 0 1 0 = t [p0L/6 - SL 2 /12 -p0L/3 + p0L p0L/6 + SL 2 /12] Hence P1 = (tL/6)( p0 - SL/2) P2 = (2tL/3)p0 P3 = (tL/6)(p0 + SL/2) Using p1, p2 and p3 to represent the intensity of the distributed load at points 1,2 and 3, noting that tL is the edge area, we get P1 = (1/6)p1 × edge area P2 = (2/3)p2 × edge area P3 = (1/6)p3 × edge area These are the components of the consistent load vector. They are statically equivalent to the distributed load and also do the same amount of work during a displacement. The work equality was based on an edge displacement function which is consistent with the assumed displacement function for all points in the element, hence, the designation “consistent load vector”. 7. Concluding Remarks A mathematical formulation of stiffness matrix, stress and a consistent load vector is proved for the linear strain triangle element LST simple but efficient in modelling multiple materials with good results compared to the simpler CST element and more difficult and time-consuming advanced elements. the formula extracted for American Academic Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2021) Volume 83, No 1, pp 1-11 11 consistent load vector for the LST is readily extended for the conical axe-symmetric pressure vessel element for convenience, [8] through [13]. References [1] Desai, C. S., and J. F. Adel, Introduction to the Finite Element Method, (New York: Van Nostrand Reinhold, 1972) [2] Zeinkiewicz. O. C.M, The Finite Element Method in Engineering Science, Second edition (New York: McGraw Hill, 1971) [3] William H. Bowes and Leslie T. Russel, Stress Analysis by the Finite Element Method for practicing Engineers, Lixengton Books. D.C. Heath and Company, Toronto, 1975 [4] Mohamed Almheriegh. “Initial Strain Problems Formulation of Combined Materials” IOSR Journals of Mechanical and Civil Engineering, IOSR Journals, International Organization of Scientific Research,. Volume: 11 issue: 1, pp 37-47 (Version -4) e-ISSN: 2278-1684 p-ISSN: 2320-234X. [5]. S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods, Springer Science + Buiness Media LLC. 3rd eds., 2008. [6]. McNeal R. H., (1978) A Simple Quadrilateral Shell Element, Computers and Structures, Vol. 8, pp. 175-183. [7]. Clough R. W. and Tocher J. L, (1965) Finite Element Stiffness Matrices for Analysis of Plate Bending, Proc. Conference on Matrix Methods in Structural Mechanic, WPAFB, Ohio, pp. 515-545. [8]. Green B. E., Strome D. R., and Weikel R. C.,( 1961) Application of the stiffness method to the analysis of shell structures, Procedures on Aviation Conference, American Society of Mechanical Engineers, Los Angeles, arch. [9]. Bazeley, G. P., Cheung Y. K., Irons B. M. and Zienkiewicz, O. C., (1966) Triangular Elements in Plate Bending, Confirming and Non – Confirming Solutions, Proc. 1st Conference on Matrix Methods in Structural Mechanics, pp. 547-576, Wright Patterson AF Base, Ohio. [10] Cleveland State University, FE course https://academic.csuohio.edu/duffy_s/CVE_512_8.pdf [11]. Cowper, G.R.: Gaussian quadrature formulas for triangles. Int. J. Numer. Methods Eng. 7, 405–408 (1973) [12]. Flores, F.G.: A two-dimensional linear assumed strain triangular element for finite deformation analysis. J. Appl. Mech. 73, 970–976 (2006) [13]. Cook, R.D.: On the Allman triangle and related quadrilateral element. Comput. Struct. 22, 1065–1067 (1986) https://academic.csuohio.edu/duffy_s/CVE_512_8.pdf