Numerical computation of shear stress in a complex cross section subjected to a shear force and application to airfoils

– The aim of this work is to develop a new numerical computational program for determining the distribution of the shear stress in a complex cross section subjected to the application of a shear force and applications to airfoils, in light to determine the position and the value of the maximum stress, and consequently the determination of the shape factor seen their practical importance. In the general case, there may be two shear forces along the two principal axes of inertia, and consequently two shape factors can be found from those two directions. The calculation of the geometric characteristics of the section is necessary in this case. The discretization is based on the development of a new technique by dividing the upper part of the section in adjacent triangles common at a point to determine the static moment of this part and the length of the band width. The validation of the program is certiﬁed by the convergence of the numerical accurate results to those for the circular section. The exact solution exists for the circle. In this case, the relative error for high discretization approaches to zero. The application will be for airfoils, considering their practical interests in the ﬁelds of aeronautics and industrial mechanics.


Introduction
The dimensions of the structures play a very important role to know the magnitude of internal forces coming into being, and therefore the value of the considered stress [1][2][3][4][5][6][7][8].If one of the dimensions of the structure is large compared to the transverse dimensions, one can consider the structure as a beam.In the bending theory, bending moment is always accompanied with the shear force that gives direct normal stress σ to the section, while the shear force T gives a stress τ working in the plane of the section in the direction of T [1][2][3][7][8][9][10][11].
The calculation of τ due to the shear force T is neglected if the longitudinal length of the beam is very large compared with the transverse dimension [6,8,9,[12][13][14][15][16][17][18].But if the length is smaller comparing with the transverse dimension [19], which is the case for many of the structures relating to the construction of an aircraft, in this case we can not neglect the stress τ due to the effect of shear force T .These structures can be found for a Corresponding author: z toufik270169@yahoo.frexample at the blades of a compressor of an engine, the propeller blades and other structures.
The shearing force T is an internal effort that will be determined from the external forces using the method of sections.For example, the flow around the propeller blades, causes a pressure distribution which will generate pressure forces and friction around the blade propeller.The projection of these forces under the coordinate axes gives the forces responsible for these shear forces [10,20,21].
The calculation of the direct stress and deflection of a beam under the isotropic bending moment effect is the object of reference [9].The author neglected in this case the calculation of the shear stress due to the effect of shear force, since it was assumed that the beam length is very large relative to the transverse dimension of the section.The calculation done in this case is purely analytical.
The energy method for the determination of the shear stress in circular section beams is presented in reference [22].Given the simplicity of the section, the calculation is purely analytical.The goal is to calculate the Article published by EDP Sciences The calculation of the shear stress only for a halfelliptical cross-section of a beam is disclosed in reference [12].The author made an application to the elliptical section which is also considered as a very simple section, where the analytical calculation is discussed.
Reference [23] is to treat the shear coefficient calculation problem in thin-wall beams.He has given a general formula of calculation and application is made to a circular section.This study aims to calculate the movement of the shearing effect of the beams by the energetic method.
Analysis of the shear stress in the short length of beams is processed first in reference [19].The application is also made for a simple circular section.The method of calculating opted in this reference is the direct analytical integration of static moment.
In all the references mentioned previously [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]22,23], the sections chosen by the authors are simple sections that are circular or elliptical.The calculation is performed purely analytical, as there is no complex section studied to merit the development of a numerical calculation program and applications in all domains of interest.The calculation method is based on the direct integration of the static moment, as the function defining the boundary circle or ellipse is given analytically.
Similarly these references provide the advantage in the calculation of direct stress due to the bending moment that the calculation of the shear stresses due to the effect of shear force.This case has of practical value when the application requests a longitudinal dimension much greater relative to the transversal dimension of the section.We find these structures, for example, in civil engineering application; but the application in aeronautics and aerospace industry, we can find beams with fairly moderate size, which can not be applied the hypothesis opted in another discipline.Must be taken in this case, the calculation of the shear stress due to the effect of shear effort, which becomes an importance value.The τ calculation negligence can cause material rupture, as the total stress balance is not well-determined [1-19, 22, 23].
The aim of this work is to develop a new numerical computing program to assess the distribution of the stress τ due to the effect of shear force T by a new technique in view of the determination of the position and the value of the maximum stress for any complex section, in order to calculate the maximum shear stress within the bending strength test for the beams of small or moderate lengths compared to transverse dimension of the section.
In fact, two programs were performed in this context.The first program aims to calculate the shear stress due to the effect of the horizontal shear force and the second program is designed to calculate the stress due to the vertical shear force, as in the general case one can have both directions along the coordinate axes.
The computer computation time depends essentially on the chosen discretization problem by the number of triangles selected in the upper area, of the section, and on the selected step on the direction of the coordinate axis as well as on the characteristics of the selected computer, including the CPU Time speed.However, comparing between the execution times of two programs implemented, the calculation of shear stress for T x is quite superior to that of the calculation given by the stress of shear force T y due to there are other intermediate procedures performed as the inverse cubic spline interpolation.
The application will be made for airfoils, wich are very interested in the fields of aeronautics and aerospace industries.The developed technique in this contest is to make a discretization by an important number of adjacent triangular meshes; intersected all in one point (one observer) and other points will be placed on the boundary of the upper part.The static moment of that part is obtained as the sum of static moments of all triangles that form this region.Here the triangles are defined by the positions of their three nodes.The major difficulty for the selected section is to ascertain the positions of the nodes in the boundary area of each evaluation of the stress position.
The comparison and validation of our numerical results will be made with implementing the program for the circular section where the exact solution exists.

Mathematical formulation
The shear stress τ due to the shear force T for sections having one symmetrical axis is given by the following relationship [1-19, 22, 23]: In formula (1) shear force T and the stress τ are directed along the same axis.The calculation of the static moment S * , the length of the band width b and the moment of inertia I are relative to the perpendicular axis of T .
The plan of the selected section is Gxy.In relation ( 1), we can have two possible cases in the direction of T as the principal axes of inertia.If T is directed along the axis Gy, then in Equation ( 1) we have T = T y , τ = τ (y), The right member of the relation (1) depends only on the geometry of the section.
In most references [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]22,23], the authors are interested in simple sections such as circle, square and triangle.In this case, the calculation of the stress distribution becomes very simple and can be found analytically.But if the section becomes complex, as the case of airfoils, interested in industrial mechanics, numerical calculation is necessarily given the complexity of the shape of the section.As a result one can determine the position and value of the stress τ Max and the corresponding shape factor by the following equation [1][2][3][4][5][6][7]10]: As seen two shear forces T x and T y , then we will have two shape factors k x and k y , respectively.

Calculation procedure
First, the calculation of (x G , y G ), A, I x and I y of the complete section are necessary.The method of calculation of these parameters is done by the technique of subdivision section by adjacent triangles like this Figure 1 [24].The accuracy of the calculation depends on the used mesh.This figure was presented with a large mesh and the other with a fine mesh.
According to formula (1), we are interested in determining two parameters b and S * who are to be calculated for each τ value and depend on the distance of the bandwidth relative to the reference section.Consider the case of T = T y directed vertically as Figure 2 shows.
The static moment S * and the bandwidth b of the upper part, following discrete in NT adjacent triangles, which are remoted from there over the Gx axis, are given by: b With The triangle number (i) is limited by the nodes number P , i and i+1.
The number of nodes in the upper part is equal to NT + 1.The position of point P must be made on the bandwidth of the upper part.One can take the point P on the main vertical axis.Then x P = 0.0 and y P = y.
To determine the value of b by the relation (3), must be determined for each top section the coordinates of points numbers 1 and NT + 1 in Figure 2. The point number 1 is located between the right end and the y Max of upper surface of the section and point number NT + 1 is located between the left end and y Max end of the function.Hence the interest to determine the maximum of the function and points numbers 1 and NT + 1 by the bisection method [25,26].All triangles in the upper part are shared with a point like this Figure 2.
Replacing the values given by relations (3) and (4) into relation (1), one can find the value of τ at that position there.For each ordered y, then consider a new upper surface and in this case it is necessary to find the position of all points are training the boundary of this surface relative to the benchmark definition of the section.It should be noted that if the surface of the upper part increases, also increases the number of discretization triangle train to keep the precision of calculation.
If we want to calculate the position of the nodes compared to the benchmark xGy, in this case y G = 0.0 in formula (5).The numbering of nodes and triangles is counterclockwise as presented in Figure 2. The surface area of the triangles given by Equation ( 6) does not depend on the selected reference.
On the top and bottom ends where S * = 0.0 and b = 0.0, then the stress τ = 0.0.
The static moment of the lower part of the section is therefore equal less static moment of the upper part because the sum of the static moments must equal the static moment of the entire section and since the calculation is made in relation to the central axis, then it is necessary that the result should be zero.The accuracy of the calculation of S * depends on the NT number chosen on the upper part and the accuracy of calculation of τ distribution depends again on the step selected on the y axis.The τ Max stress can be evaluated by localized on.Finally we can calculate the corresponding shape factor using the relation (2).Now consider the case of the application of the horizontal shear force T = T x .In this case it is necessary to change x by y and y by x in Equations ( 3)-( 6) to obtain the following relations ( 7)-( 9).The relation ( 6) remains always valid.
The stress will be evaluated based on the x axis.Must vary the x-coordinate and determine the y by the inversion of the function defining the geometry using the bipartition method [25,26].This method needs to locate the solution provided in a defined interval.In this case the computation time of the computer by determination of stress becomes much higher relative to the time of calculation for the case of the application of shearing force T = T y .

Applications
The circle is a case for the validation of the results found by the developed numerical calculation program.The exact results in the calculation of the stress τ to a circle of radius R are summarized by [1-19, 22, 23]: Then Consequently Then, the shape factors are given by: According to the relationship (15), the maximum shear stress in a circle is 33% greater than the average stress.As a result, the application will be made on a choice of 10 symmetrical airfoils.For the latter, the boundary is usually given by tabulated values [27,28].To find an analytical form for the upper and lower surfaces was used a cubic spline interpolation [25,26].Airfoils chord is taken as C = 1.00.

Error of numerical calculation
To evaluate the error of numerical calculation compared with exact results of a circle using Equation ( 16) below: The word parameter in Equation ( 16) can be τ at a point, τ Max or k.

Results and comments
In Table 1, there is presented the discretization effect on the convergence of the maximum shear stress.The selected section is a circle of radius R. The number of triangles has been varied in the upper half of the circle, since the maximum stress is located at the axis of symmetry Gx.It is noticed the convergence of τ Max to the exact solution that is shown in Table 1.We also note the rating decrease the error with increasing of number of triangles NT .For NT = 3000 an error was found ε = 10 −4 .
The circular shape is discretized with the same manner as the airfoil to keep the same calculation method of the stress distribution.
The airfoils selected for calculation are summarized in Table 3.In Tables 5, 6, 9 and 10 are present only the section number and must refer to Table 3 to see the airfoil name.
Only for the NACA 0012 airfoil, it was found in reference [29], the function defining the upper and lower surfaces by the following equation.With: For this airfoil, the maximum thickness is t = 12%, or equal to 0.12.It is located at a distance of x = 0.2991 from the leading edge.To the lower surface, it is considered the sign (-) for reasons of symmetry.Here, the chord of the airfoils is equal to unity.Table 2 shows the minimum number of triangles considered in the discretization of the upper part of the section to find a number of significant digit after the decimal point in the calculation of the stress τ Max circle of Table 1.For example, to have 5 exact decimal numbers, we must need at least 515 triangles in the upper part of the section.
Figure 3 shows the variation of the obtained relative error of the discretization according to the number of selected triangles in the upper part of the section to obtain the value of τ Max of a circle.We clearly noticed that the more triangles increase, the relative error decreases gradually to an accepted error.Then the relative error getting mainly depends on the number of triangles of the discretization.For example, to obtain an error ε = 10 −5 %, we need about NT = 237 triangles.The points of definition of airfoils, RAF 30 Modified, NACA M1 and NACA 62 are presented respectively by 25 points each as the watches Table 4 [27,28].
According to reference [24], the geometrical characteristics of airfoils of Table 3 are shown in Table 5.The calculation of x G is made with respect to the leading edge of the airfoil and I x and I y are calculated relative to the central axes.
The y Max values required for the change in the ordinate y to calculate the shear stress τ shown in Figure 4 for the airfoils of Table 1 are summarized in Table 6.This value represents half of the maximum thickness of the airfoil.
In Figures 4a and 4b there is shown the variation of the length of the bandwidth and the static moment of the upper surface of the NACA 0012 airfoil when T = T y .The values of these curves are shown in Table 7.It is clearly noted that the bandwidth length and the static moment of the upper and lower ends are equal to zero.
In Figures 4c and 5c there is shown the variation of the shear stress on the NACA 0012 airfoil surface respectively when applying a shearing force T = T y and T = T x .Note that the shear stress is zero at the left and right ends when T = T x and zero at top and bottom ends when T = T y since in these ends the static moment is zero.Then the maximum stress has a value of an internal point of the surface of the section.τ Max position is not located at the section center of gravity.
In Figure 4, the representation of the airfoil is not in an orthonormal system for a good presentation of the change in shear stress.
In Figure 4c, the variation of τ is symmetric since the airfoil is itself symmetric.It only varies with the thickness  For the evaluation of the stress τ in Figure 5 for other airfoils, the x-coordinate must vary between x G and C − x G , such as x G shown in Table 5.
The numerical results in some points of the variation of the stress shown in Figure 4 are shown in Table 7. Then the maximum stress is equal to 14.919 and is located at a distance equal to 0.0066 with respect to the horizontal axis.
Figures 5a and 5b  in b in Figure 5a shows the variation of thickness across the length of the airfoil.
In Figure 5c, the variation of τ is a function of the abscissa and is constant along the ordinate.
The numerical results in some points of the change in the shear stress shown in Figure 5c are shown in Table 8.Then the maximum stress is equal to 16.268 and is located at a distance of 0.512 from the loading edge.
Table 9 shows the value and position of τ Max and the shape factor value of airfoils presented in Table 3 when applied a shear force T = T y .The τ Max position is given relative to the horizontal axis of symmetry.
Similarly in Table 10 are shown the value and position of τ Max and the shape factor values of airfoils presented in Table 3 when applying a shear force T = T x .The τ Max position is given relative to the leading edge of the airfoil.
We note from Tables 9 and 10 that the values found are of the same order of magnitude for the airfoils.
The selected number NT for evaluation of the values presented in Tables 9 and 10 is equal to NT = 10 000 to see 5 exact decimal digits.
Finally for moderate lengths of beams subjected to bending we can apply a criterion of resistance to avoid breakage of the beam [1-7, 10]: In relation (18) we see clearly the influence of τ Max in the balance of applied stresses, which is accompanied with σ Max stress due to the bending moment.Applying strength criterion (18) we will determine the dimension of the section with a high precision without breaking.That is to say the dimension of the section with consideration of τ Max is less than the dimension of the section without effect τ Max .
The shape factors presented in Tables 9 and 10 show the ratio of the maximum stress compared with the average stress applied to the section.From the resulting of the shape factors presented in Tables 9 and 10, we note that the maximum shear stress in an airfoil is between 16% and 77% greater than the average stress for the selected airfoils.
Whatever the shape of the beam section, the stresses σ Max and τ Max are proportional to the dimension of the section by the following relation [1][2][3][4][5][6][7]10]: So we see from the relation ( 19) that we cannot neglect the shear stress τ Max due to the effect of shear force in relation to the direct stress of bending moment if the length of the beam is the same order of magnitude as the transverse dimension of the section, which is the case for a variety of domain structure of the aerospace and industrial mechanics.

Conclusion
This work allows us to calculate the shear stress due to the effect of shear force in complex sections.We can draw the following conclusions: For non-symmetrical sections, we must seek the principal axes of inertia and consider the calculations along these two axes.
The numerical program can make any simply connected symmetric section, including airfoils, the shape of the boundary is given by tabulated Each section has two shape factors k x and k y following the two options of two shear forces T x and T y .
For T = T y , the stress τ will be forced to zero at the top and bottom ends and is zero at the left and right ends if T = T x .
When T is directed arbitrarily in the plane of the section, must be the projection of T as the principal axes of inertia and consider the calculation of the stress τ along the two axes and make the vector summation.
The stress τ Max is at a point which is not the section center of gravity.
The interpolation of tabulated values is necessary to find an analytical form the upper and lower surfaces.In our study we opted cubic spline interpolation.
The calculation must be made to the two principal axes of inertia.
The calculation of the geometric characteristics of the section is necessary.
The value of T in the program does not matter, however, calculating the value of τ /T as a function of the geometry of the section.
The convergence of the results depends on the discretization, including the number of triangles in the upper part and the step used in the direction of T .
Whatever the shape of the section, the maximum shear stress in this section is always greater than the maximum stress of circumcircle for this section.
As perspective and in the same line of research, we can consider the calculation of the shear stress due to shear forces for thin-walled hollow sections.The geometrical characteristics of this type of section applied to airfoils, necessary to this study are presented in reference [30].

A
. Kirad et al.: Mechanics & Industry 17, 608 (2016) Nomenclature T Shear force along the principal axis of inertia S Static moment with respect to the principal axis b Band width of the upper part I Principal central moments of inertia τ Shear stress σ Normal stress due to bending moment x, y Coordinates of a point A Area of the section of the upper part NT Number of triangles in the discretization k Shape factor of the section C Chord of the airfoil ε Computation relative error R Radius of the circle P Point on the band width of the section upper part L Length of the beam h Transverse dimension of the section ai Coefficients for NACA 0012 airfoil equation m Number of significant decimal digit of the results Superscripts * Upper part of the section Subscripts Max Maximum G Center of gravity (i) Counter on the triangles i Counter on the nodes x Along the main axis Gx y Along the main axis Gy el Eligible weighted coefficient in the energy equation for the shear stress.

Fig. 1 .Fig. 2 .
Fig. 1.Discretization of the full field of the section in view of determination of the geometrical characteristics.(a) Large mesh.(b) Fine mesh.

Fig. 4 .
Fig. 4. Variation of the bandwidth, static moment and the shear stress τ due to the shear force Ty on the airfoil surface.(a) Variation of bandwidth b.(b) Variation of the static moment S * .(c) Variation of the shear stress τ .
and constant along the horizontal axis.In this figure we have not presented the airfoil full scale for aim of a good visualization of the variation of the stress.The calculation of the stress exactly to the high and low ends for T = T y , and the left and right ends for T = T x , the program developed gives a problem of division by zero.Exactly there will be a kind of indeterminacy (0/0) since in this case S * = 0.0 and b = 0.0 then S * /b = 0/0.To remove this uncertainty we consider that ends in (y Max − ε) and in (−y Min + ε) for T = T y and in points (x Max − ε) and (−x Min + ε) for T = T x with ε = 10 −8 .In this case it is found that S * = 0.0 and b = 0.0 wherein the ratio of S * /b = 0.0 but close to zero and therefore the value of τ ≈ 0.0.So it follows that the value τ at the ends is exactly zero.

Fig. 5 .
Fig. 5. Variation of the bandwidth, static moment and the shear stress τ due to the shear force Tx on the airfoil surface.(a) Variation of bandwidth b.(b) Variation of the static moment S * .(c) Variation of the shear stress τ .

Table 1 .
Discretization effect on the convergence of the results for the circle of radius R.

Table 2 .
Minimum number of triangles NT giving the number of significant digit m of the results of τMax for the circle.
NT Fig.3.Relative error according to the number of triangles of the discretization of the upper part of the section for a circle.

Table 3 .
Names of studied airfoils.

Table 4 .
Defining points of the surface of the airfoils RAF 30 M, NACA M1 and NACA 62.

Table 5 .
Geometric characteristics of the airfoils of

Table 6 .
Values of yMax for airfoils of Table3.
A i r f o i l

Table 7 .
Discretization for the airfoil NACA 0012 when T = Ty.

Table 8 .
Discretization for the airfoil NACA 0012 when T = Tx.

Table 9 .
Numerical results for the airfoils of Table3when T = Ty.

Table 10 .
Numerical results for the airfoils of Table3when T = Tx.