Nonuniform Torsion Without Shear Deformation Effect Using FEM

SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 8 NONUNIFORM TORSION WITHOUT SHEAR DEFORMATION EFFECT USING FEM Dang-Bao TRAN1,2, Jaroslav NAVRÁTIL1 1 Department of Structures, Faculty of Civil Engineering, VSB–Technical University of Ostrava, Ludvíka Podéště 1875/17, Ostrava, Czech Republic. 2 Department of Civil Engineering, Faculty of Architecture, Thu Dau Mot University, Tran Van On 06, Binh Duong Province, Vietnam. dang.b[email protected], jaroslav[email protected]. DOI: 10.35181/tces-2021-0006 Abstract. This paper presents the use of a finite element method (FEM) to analyze nonuniform torsion with an arbitrary cross-section with homogeneous elastic material without shear deformation effect. Beams are constrained by the most common types of supports, such as fixed, pinned, and roller, and are subjected to any applied torsional load, concentrated, and distributed. The presented FEM transforms the 3D analysis of nonuniform torsion beams into separated 2D cross-sectional and 1D modeling. The geometric constants of the cross-section are firstly derived from 2D FEM, which uses a 9-node isoparametric element. Then, a 1D FEM, which uses the Hermitian shape function, is developed for computing the twist angle, the derivative of twist angle, and stress resultants. Finally, the stress field is obtained from the local analysis on the 2D-cross section. A MATLAB program is executed to validate the numerical method through examples. The validation examples have proven the reliability of the author's numerical method for analyzing problems defined above. Keywords Nonuniform torsion, thin-walled structure, shear deformation effect, warping restraint, finite element method. 1. Introduction In engineering practice, we often encounter beam structures loaded in torsion. When the warping of a member's cross-section is not restrained, the stress field of a prismatic beam with homogeneous isotropic material can be derived from the Saint-Venant theory strictly [1-4]. In practice, because i) boundary conditions are imposed, ii) geometrical characteristics of the section, the warping of beam's cross-section is restrained, which leads to additional normal and shear stresses, which the Saint Venant theory does not take into account. Vlasov [5] was the first to formulate the problem of nonuniform torsion. Benscoter [6] improved Vlasov's theory, which neglects the shear deformation effect, leading to errors with closed cross-sections [7]. Many authors established FEM to consider nonuniform torsion taking into account the shear deformation effect [714]. However, the above research [7-14] used the approximations of Thin Tube Theory [5] to determine the bar's geometric constants, restricting the accuracy and applicability of the formulations. El Fatmi [15,16] proposed a beam theory with a nonuniform warping, including the effects of torsion and shear forces for arbitrary cross-sections made of homogeneous isotropic elastic material. Mokos and Sapountzakis [17] presented a nonuniform torsion theory considering the shear deformation effect for general cross-sections implemented by Boundary Element Method. The purpose of this paper is to establish a numerical method using FEM to solve the nonuniform torsion problem without the shear deformation effect of the prismatic beam with arbitrary cross-section unchanged throughout the length with homogeneous isotropic elastic material, based on displacement and strain fields derived from Sapountzakis, EJ et al. [17]. This paper is a preliminary research step for future research, which considers the shear deformation effect in nonuniform torsion. 2. A brief introduction of the theoretical formulas for nonuniform torsion Let us consider a prismatic beam with arbitrary crosssection, constant along the length L with the modulus of elasticity E, and shear modulus G. The longitudinal axis is the x-axis, and the cross-sections lie in the y–z plane. S is the shear center of the cross-section. The parallel system , S zzz=+ S yyy=+ intersects at O, the arbitrary point. CXYZ is the parallel system with Sxyz through the center SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 9 of gravity C. The multiply connected domain Ω is bounded by n curves, 1,Γ 2,Γ..., 1, n− Γ , n Γ as Fig. 1. Tangent vector t with associate coordinate s and normal vector n set up the right-handed system. The beam is enforced to the arbitrary loads, distributed torque mt(x), concentrated torque Mti(x), concentrated warping moment Mwi(x). The cross-section is assumed with no distortion. Fig. 1: Cross-section of the prismatic beam The displacement field is expressed as [17] (,,)uxyz ='(,) (,), PS xS S yz yz θφ φ + (1) (,) (), x vxz z x θ =− (2) (, ) (), x wxy y x θ = (3) where u, v, w are the axial and transverse displacements of the beam with respect to Sxyz x θ is the angle of twist (,) P Syz φ is the primary warping function with respect to shear center S derived from 20in on . P S PP SS yzyzn nnznyn yz φ φφ ∇= Ω ∂∂ +=− Γ ∂∂ (4) Moreover, the restrained warping function can be computed as follows by using the transformation of coordinates [18] 1 (,) (,) (,) . PP P SOSS O yz yz yZ zY yzdA A φφ φ Ω =+−−  (5) where (,) P Oyz φ is the warping function with respect to 𝑂𝑦𝑧 system coordinates, A is the area of the cross-section S S φ is the secondary warping function with respect to the shear center S, given as ''' 2() in 0on , SP x SS SS SS yz n Ex G nn yz θ φφ φφ ∇=− Ω ∂∂ +=Γ ∂∂ (6) The stress fields obtained from the theory of elasticity as [17] "() (,), P xx x S Ex yz σθφ = (7) xy τ ='() P xy P S x Gx z y τ φ θ  ∂−+   ∂  144424443 , S xy S S Gy τ φ ∂ ∂ 123 (8) xz τ ='() P xz P S x Gx z z τ φ θ  ∂++   ∂  144424443 , S xz S S Gz τ φ ∂ ∂ 123 (9) In addition, the shear stress in (8), (9) can be seen as composed of two parts: primary shear stress and secondary shear stress. The warping moment, the total twisting moment, the primary twisting moment, the secondary twisting moment are denoted by , w M , t M , P t MS t M obtained as [17] w M=P xx S d σφ Ω −Ω= '' , Sx EC θ − (10) ', PP PP P SS txy xz tx M zyd yz GI φφ ττ θ Ω   ∂∂ =−++Ω    ∂∂    =  (11) S t M= PP SS SS xy xz yz φφ ττ Ω  ∂∂ −− =   ∂∂  ''' (), S EC x θ − (12) , P S tt t MM M=+ (13) where t I is the torsional constant, determined as t I =22 , PP SS yzy z d zy φφ Ω  ∂∂ ++ − Ω   ∂∂   (14) S C are the warping constant determined by () 2. P SS Cd φ Ω =Ω  (15) The relation between the stress resultant, total twist moment Mt, and the distributed torque mt(x) is expressed as [17] (). tt Mmx x ∂=− ∂ (16) After some algebra, the equilibrium equation for the nonuniform torsion problem of a homogeneous isotropic bar without shear deformation effect is expressed as 42 42 . xx Stt dd EC GI m dx dx θθ −= (17) SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 10 3. FEM procedures 3.1. Angle twist, x θ Applying Galerkin’s method, one of the methods of weighted residual, the governing equation Eq. (17) is transformed to weak form as 42 42 0 0, Lxx Stt dd IEC GI mdx dx dx θθ η  =−−=     (18) where η is the test function. The beam is discretized into a number of finite elements. After some manipulation, the weak formulation of Eq. (18) can be expressed as 22 22 1ee nxx St LL i dd dd IEC dxGI dx dx dx dx θθ η η =  =+−     0 0, e L x tttw L ddd mdx GI M M dx dx dx θηη ηη  +− + − =     (19) where Le is an beam element domain n is the number of elements for the beam. The Hermitian polynomial interpolation function, which achieves the correct solution of the problem with refining mesh, proves more flexible than the Hyperbolic interpolation function, which gives the most accurate results with the most analytical solution [19]. We choose cubic functions for the spatial interpolation of the twist angle, x θ , in terms of nodal variables. To this end, we consider an element that has two nodes, one at each end. The twist angle can be expressed as 12 112 3 24 () () () () , xx xx x dd Hx Hx Hx Hx dx dx θθ θθ θ =+ + + (20) where 23 23 12 23 3 23 23 34 23 23 32 2 () 1 , () , 32 () , () , x xxx Hx Hx x l ll l xx xx Hx Hx ll ll =− + = − + =− =−+ (21) 12 12 ,,, x x xx dd dx dx θθ θθ are the nodal degrees of freedom. Inserting the Eq. (20), Eq. (21) to Eq. (19) results in the stiffness matrix of the beam element 12 , ee =+KK K (22) where 111 0, e L eT S EC dx= KBB (23) 222 0, e L eT t GI dx= KBB (24) where 2 22 2 3 12 4 12222 , dH dH dH dH dx dx dx dx  =   B (25) 3 12 4 2. dH dH dH dH dx dx dx dx  =  B (26) The corresponding element nodal degrees of freedom is 12 12 ,,, . T xx xx dd dx dx θθ θθ  =  e d (27) The third term in Eq.(19) results in the element force vector. For a generally distributed torque, we need to compute 1 2 03 4 () , e L t H H mx dx H H    =     e F (28) where e F is the element force vector. The last term in Eq. (19) is the boundary conditions of the total twisting moment and the warping moment at the two boundary points, x = 0 and x = L, of the beam. If these boundary conditions are known, the known total twisting moment and warping moment are included in the system force vector at the two boundary nodes. Otherwise, they remain unknown. However, the twist angle and the derivative of the twist angle are known as geometric for this case. Assembling the element stiffness matrix and vectors leads to the system matrix equation given below .=Kd F (29) 3.2. Warping function, , P S SS φφ Using Galerkin’s method, with test function 1()H η ∈Ω and applying the Gauss-Green theorem, the governing equations Eq. (4), Eq. (6) are transformed to weak form as PP SS d yy zz φφ ηη Ω  ∂∂ ∂∂ +Ω−   ∂∂ ∂∂  (𝑛𝑧−𝑛 𝑦)𝜂𝑑𝑠 =0. (30) ''' () 0. SS P SS x S Ex dd yy zz G φφ θ ηη φη ΩΩ  ∂∂ ∂∂ +Ω− Ω=   ∂∂ ∂∂   (31) The warping function values, P S φ , , S S φ in Eq. (30) and Eq. (31) are approximated by the FEM, which is presented in [3, 4]. SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 11 4. Validation examples In this section, the accuracy of the proposed FEM has been examined. A computer code is developed in the MATLAB R2015a software based on the formulations described in the previous sections. The obtained results are compared with the available works of the literature. 4.1. Example 1 As a first example, a thin-walled beam with an I-shaped cross-section shown in Fig. 2 is analyzed and compared with the results obtained by Tralli, A.M. [8] and Kim, N. I. et al. [11]. The I-section beam is clamped–free and the concentrated torsional moment Mx = 25 kNm is applied at its free end. The length of the beam is 5 m. The Young’s modulus and the Poisson’s ratio are E = 2 x 106 MPa and 𝜈= 0.3 , respectively. The first test was analyzed using 16 axial elements and 19 elements in the cross-section. The geometric constants of the cross-section determined from the present study are as follows 76 2.1559 10 S Cm − =× , 64 2.8643 10 t Im − =× . Fig. 2: A cantilever beam, the applied load, and the geometry of the I cross-section, units [mm]. Fig. 3: The variation of twist angle along the x-axis. In Table 1, the angle of twist at the free end predicted by the present study is given and compared with the results of other researchers. The variations of the angle of twist along the length of the thin-walled beam are shown in Fig. 3. It can be seen from Table 1 and Fig. 3 that the twist angle results obtained from the present study are in excellent agreement with other works. Tab.1: Comparison of the twist angle 𝜃 at the free end calculated by different methods Methods 𝜃 [𝑟𝑎𝑑] at x = L Present study 0.5171 Tralli [8] 0.5183 Kim et al [11] 0.5173 Vlasov [11] 0.5167 Fig. 4: Distribution of the axial stress 𝜎 in 2D cross-section at fix-end, unit MPa, Fig. 5: Maximum of normal stress 𝜎 along the axis x of the beam. Fig. 4 shows the contour plot of the axial stress 𝜎 at the clamped end of the I-section beam. Fig. 5 depicts the variation along the axis x of the normal stress 𝜎. Fig. 6 and Fig. 7 present the variation along the axis x of the warping moment and the torsional moments, respectively. It is observed from Fig. 4 that the variations of the axial stress are symmetric to the area center of the cross-section. In contrast to the flanges, no significant axial stress appears at the web of the I-section beam. The maximum value of 𝜎 = 515.911 MPa occurs at the corners of the SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 12 flanges. In Fig. 5, the variations of the maximum axial stress 𝜎 at the fixed-end of the cross-section along the length of the beam are shown and compared to the ones obtained from [8, 11]. It can be seen that the axial stress predicted by the present study is in good agreement with the results of Tralli, A. M. [8]. Fig. 6: The variation of warping moment Mw along the axis x of the beam. Fig. 7: The variation of torsional moments along the axis x of the beam. 4.2. Example 2 A cantilevered thin-walled beam with a channel crosssection depicted in Fig.8 is analyzed in this section. The depth, width, and thickness of the cross-section are h = 0.833 m, b = 0.917 m, t = 1/6 m, respectively. The length of the beam is L = 4 m and, it is made of a material with the modulus of elasticity E = 2 x 105 MPa and the shear modulus G = 0.77 x 105 MPa. The considered thin-walled beam is subjected to uniformly distributed twisting moment equal to mx = 4.07 x 103 kNm/m. The second example was analyzed using in the present study by employing 80 axial elements and 64 elements in the crosssection. The geometric constants of the cross-section determined from the present study are as follows 6 0.0059 S Cm=, 4 0.00407889 t Im=. Fig. 8: Cantilever beam, the applied load, and the geometry of the channel cross-section, units [mm]. Shakourzadeh, H. et al. [7] used FEM based on Benscoter's model to solve this problem and compared it with the results obtained from Vlasov's theory, which neglects the shear deformation effect. Meanwhile, Sapountzakis, E. J. et al. [20] used Analog Equation Method (AEM) and AEM with isogeometric analysis with 3 B-spline types: i) Quadratic B-spline, ii) Cubic B-spline), iii) Quadratic B-spline to investigate this problem. Fig. 9: The variation of twist angle along the x-axis. Fig. 10: The variation of derivative of the twist angle along the x-axis. SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 13 Fig. 9 and Fig.10 show the variation of the angle of twist and the derivative of twist angle along the x-axis, respectively. The value of the angle of twist and the derivative of twist angle at the free end, and the warping moment at the fixed end are presented in Table 2 and compared with the ones obtained from Shakourzadeh, H. et al [7], Sapountzakis, E. J. et al. [20]. Fig. 11: The variation of warping moment along the x-axis. Fig. 12: The variation of the torsional moment along the x-axis. Tab.2: Comparison of the twist angle, the derivative of twist angle, and warping moment between the methods Methods 𝜃 [rad] at x = L 𝜃 󰆒 [rad/m] at x = L Warping moment Mw [N.m2] at x = 0 Present study -0.0429 -0.0114 -19.083 x 106 SaintVenant model [7] -0.103 - - Vlasov model [7] -0.045 -0.012 -18.33 x 106 FEM-Benscoter model [7] -0.045 -0.011 -18.17 x 106 AEM (50NP) [20] -0.050 -0.009 -16.75 x 106 AEM(Quadratic Bspline) [20] -0.039 -0.006 -13.48 x 106 AEM(Cubic B-spline) [20] -0.061 -0.013 -18.29 x 106 AEM(Quaratic Bspline) [20] -0.046 -0.008 -15.5 x 106 Fig. 11 and Fig.12 depict the variation of the warping moment and the torsional moments along the axis, respectively. It can be seen from Table 2 that the twist angle analyzed by the present study is in good agreement with the results of Vlasov model [7], FEM-Benscoter model [7], and AEM(Cubic B-spline) [20]. The distribution of the axial stress 𝜎 and shear stress 𝜏 inside the channel cross-section at the position x = 0.5 m is shown in Fig. 13. The maximum normal stress 𝜎 = 657.865 MPa occurs on the tip of the flange. As expected, the distribution of the shear stress 𝜏, with the maximum value 205.981 MPa, decreases from the outer flange edge to the inner edge [21]. (a) (b) (c) Fig. 13: The normal stress 𝜎 and b) the shear stress 𝜏 b) the shear stress 𝜏 at x = 0.5 m from the fixed-end, units [MPa]. 5. Conclusion In this paper, a finite element method is developed to solve nonuniform torsion without the shear deformation effect of the prismatic beam with an arbitrary cross-section with homogeneous isotropic material. Two examples were performed for validating the accuracy of the present study. SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 14 The comparison of the results from the present study with the corresponding results in the literature shows that the present study can predict the responses of the non-uniform torsion with arbitrary cross-section without shear deformation effect accurately. Linear analysis of nonuniform torsion considering shear deformation would be one of the future research topics. Acknowledgments The works were supported by the Student Grant Competition VSB-TUO. The registration number of the project is SP2021/77 “Nonuniform torsion in prismatic beams with arbitrary cross-sections using FEM”. References [1] TIMOSHENKO, S. P and J. N. GOODIER. Theory of elasticity. New York, USA: McGraw-Hill Book Company, 1951. [2] PILKEY, W. D. Analysis and design of elastic beams: Computational methods. New York, USA: John Wiley & Sons, 2002. ISBN: 0-471-38152-7. [3] TRAN, D. B., J. NAVRÁTIL and M. ČERMÁK. An efficiency method for assessment of shear stress in prismatic beams with arbitrary cross‐sections. Sustainability (Switzerland). 2021. Vol. 13, no. 2, pp. 1– 20. ISSN 20711050. DOI: 10.3390/su13020687. [4] TRAN, D.-B. Torsional Shear Stress in Prismatic Beams With Arbitrary Cross-Sections Using Finite Element Method. Stavební obzor - Civil Engineering Journal. 2021. Vol. 30, no. 2. ISSN 1805-2576. DOI: 10.14311/cej.2021.02.0030. [5] VLASOV, V. Z. Thin walled elastic beams. Israel Program for Scientific Translations. Jerusalem, Israel. 1961. [6] BENSCOTER, S. U. A Theory of Torsion Bending for Multicell Beams. Journal of Applied Mechanics. 1954. Vol. 21, no. 1, pp. 25–34. ISSN 0021-8936. DOI: 10.1115/1.4010814. [7] SHAKOURZADEH, H., Y. Q. GUO and J. L. BATOZ. A torsion bending element for thin-walled beams with open and closed cross sections. Computers and Structures. 1995. Vol. 55, no. 6, pp. 1045–1054. ISSN 00457949. DOI: 10.1016/0045-7949(94)00509-2. [8] TRALLI, A. A simple hybrid model for torsion and flexure of thin-walled beams. Computers and Structures. 1986. Vol. 22, no. 4, pp. 649–658. ISSN 00457949. DOI: 10.1016/0045-7949(86)90017-9. [9] BACK, S. Y. and K. M. WILL. A shear-flexible element with warping for thin-walled open beams. International Journal for Numerical Methods in Engineering. 1998. Vol. 43, no. 7, pp. 1173–1191. ISSN 00295981. DOI: 10.1002/(SICI)10970207(19981215)43:7<1173::AID-NME340>3.0.CO;2-4. [10] Erkmen, R.E. and M. Mohareb. Torsion analysis of thin-walled beams including shear deformation effects. Thin-walled structures. 2006. Vol. 44, no. 10, pp.10961108. ISSN 0263-8231. DOI: 10.1016/j.tws.2006.10.012. [11] KIM, N. I. and M. Y. KIM. Exact dynamic/static stiffness matrices of non-symmetric thin-walled beams considering coupled shear deformation effects. ThinWalled Structures. 2005. Vol. 43, no. 5, pp. 701–734. ISSN 02638231. DOI: 10.1016/j.tws.2005.01.004. [12] MURÍN, J. and KUTIŠ, V. An effective finite element for torsion of constant cross-sections including warping with secondary torsion moment deformation effect. Engineering Structures. 2008. Vol. 30, no. 10, pp. 2716–2723. ISSN 01410296. DOI: 10.1016/j.engstruct.2008.03.004. [13] MURÍN, J., V. KUTIŠ, V. KRÁLOVIČ and T. SEDLÁR. 3D beam finite element including nonuniform torsion. Procedia Engineering. 2012. Vol. 48, p. 436–444. ISSN 18777058 DOI: 10.1016/j.proeng.2012.09.537. [14] MURÍN, J., M. AMINBAGHAI, V. KUTIŠ, V. KRÁLOVIČ, T. SEDLÁR, V. GOGA and H. MANG. A new 3D Timoshenko finite beam element including nonuniform torsion of open and closed cross sections. Engineering Structures. 2014. Vol. 59, pp. 153–160. ISSN 01410296. DOI: 10.1016/j.engstruct.2013.10.036. [15] EL FATMI, R. Non-uniform warping including the effects of torsion and shear forces. Part I: A general beam theory. International Journal of Solids and Structures. 2007. Vol. 44, no. 18–19, pp. 5912–5929. ISSN 00207683. DOI: 10.1016/j.ijsolstr.2007.02.006. [16] EL FATMI, R. Non-uniform warping including the effects of torsion and shear forces. Part II: Analytical and numerical applications. International Journal of Solids and Structures. 2007. Vol. 44, no. 18–19, pp. 5930–5952. ISSN 00207683. DOI: 10.1016/j.ijsolstr.2007.02.005. [17] SAPOUNTZAKIS, E. J. Bars under Torsional Loading: A Generalized Beam Theory Approach. ISRN Civil Engineering. 2013. Vol. 2013, pp. 1–39. ISSN 20905114. DOI: 10.1155/2013/916581. [18] GRUTTMANN, F., R. SAUER and W. WAGNER. Shear stresses in prismatic beams with arbitrary crosssections. International Journal for Numerical Methods in Engineering. 1999. Vol. 45, no. 7, pp. 865–889. ISSN 00295981. DOI: 10.1002/(SICI)10970207(19990710)45:7<865::AID-NME609>3.0.CO;2-3. [19] MOHAREB, Magdi and F. NOWZARTASH. Exact Finite Element for Nonuniform Torsion of Open Sections. Journal of Structural Engineering. 2003. Vol. 129, no. 2, pp. 215–223. ISSN 0733-9445. DOI: 10.1061/(asce)07339445(2003)129:2(215). [20] SAPOUNTZAKIS, E. J. and I. N. TSIPTSIS. Bsplines in the Analog Equation Method for the generalized beam analysis including warping effects. Computers and Structures. 2017. Vol. 180, pp. 60–73. ISSN 00457949. DOI: 10.1016/j.compstruc.2016.03.007. [21] GENOESE, A. et al. A mixed beam model with nonuniform warpings derived from the Saint Venànt rod. Computers & Structures. 2013. Vol. 121, p. 87–98. ISSN 0045-7949. DOI: 10.1016/j.compstruc.2013.03.017. SECTION BUILDING STRUCTURES & STRUCTURAL MECHANICS VOLUME: 21 | NUMBER: 2 | 2021 | DECEMBER © 2021 TRANSACTIONS OF VSB - TECHNICAL UNIVERSITY OF OSTRAVA CIVIL ENGINEERING SERIES 15