Shear and Bending Stresses in Prismatic, Non-Circular-Profile Shafts with Epitrochoidal Contours Under Shear Force Loading

1. Introduction

As a general rule, the shear stress is neglected in comparison with the bending stresses in a beam or shaft loaded with shear force. However, this is only permissible for large length ratios ($l/d$) depending on the bearing type of the shaft, where l is the length and d is the diameter of the shaft. Based on classical engineering mechanics (see, e.g., [1]), an aspect ratio limit of $l/d>6.67$ for double-sided beams (Figure 1; case 1) and an aspect ratio limit of $l/d>3.34$ for single-sided beams (Figure 1; case 2) can be determined for the stress ration limit of $\tau /\sigma <5\%$.

In many industrial applications in drive technology, bearing supports similar to case 1 are generally used, although the length ratios can also be in the range of 2.0 to 3.0. In such applications, the magnitude and effect of the transverse shear stresses should, therefore, be checked. The classical approach to calculating the shear stress, according to textbooks, assumes that the transverse shear stress also remains parallel to the shear force everywhere on the edge of the cross-section (see Figure 2). However, this assumption violates the boundary condition of the free lateral surface (e.g., [1]) and, depending on the shear contraction number, can lead to deviations of up to approx. $10\%$ for round cross-sections. Furthermore, there are no general solutions for determining the shear stresses for trochoidal profile shafts using the classical formulation.

In [2], the torsionally loaded non-circular profile cross-sections were treated in detail with the aid of conformal mapping. Suitable mapping functions were found using the successive approximation for hypotrochoidal profiles and torsional stresses, and deformations were determined on this basis.

The pure bending loads for both stationary and rotating load cases were also discussed in [3], and the corresponding solutions for the bending stresses were derived. This article, therefore, also deals with the case of bending with shear force.

2. Geometrical Description

The non-circular contours used in technical applications are often developed on the basis of general cycloids or trochoids. These curves are also referred to and classified in the technical literature as trochoids [4] and wheel lines [5]. Cycloids can be generated by rolling a circle on a path curve (usually also a circle). Depending on the rolling side, epicycloids (Fig. 1a), hypocycloids, and hybrid cycloids are produced. If the resulting curve is used as the base curve for a new trochoid, the secondary curve is referred to as a higher-level trochoid. This, in turn, can be used as a new base for further trochoids, whereby the level of the trochoid is increased accordingly.

The general representation and classification of higher-level trochoids is described in detail in [6].

Simple epicycloids

Figure 3 shows the generation of a simple epicycloid or a second-stage wheel line. Here, a circle rolls from the outside on a base circle, and point $A_2$ generates a simple epitrochoid with the eccentricity $A_1{A}_2$.

The geometry of the movement of the generating point $A_2$ can be described as follows:

$$x=r\cos \left(\mathrm{t}\right)+e_1\cos \left(\left(n+1\right)t\right)$$

$$y=r\sin \left(\mathrm{t}\right)+e_1\sin \left(\left(n+1\right)t\right)$$

(1)

Here, $r$ denotes the nominal radius, $e_1$ the eccentricity, and $n$ the periodicity (number of sides) of the cycloids. Equation (1) can be represented as follows by using the complex function:

$$z=r\zeta +e_1{\zeta }^{n+1}$$

(2)

where $z=x+iy$ and $\zeta ={\rho e}^{it}$ with $i=\sqrt{-1}$. If the generation point $A_2$ lies on the circumference of the rolling circle $K_2$, the path describes a limiting case where the curve begins to cut itself. The following applies to the eccentricity:

$$e={\displaystyle \frac{r}{n+1}}$$

(3)

For an epitrochoid with a flat point (on the flank), the following applies according to ([7]):

$$e={\displaystyle \frac{r}{{\left(n+1\right)}^2}}$$

(4)

The epitrochoid with a flat point with $n=3$ was used for the first time in industry as a non-circular polygonal shaft–hub connection. These profiles became known as K-profiles [8] and were the basis of the later standardized P3G profiles [9].

Equation (2) can "conformally" map a unit circle to the profile cross-section. This property plays an important role in determining the stress state in a profile shaft [10].

Figure 4 shows simple epitrochoids with different numbers of drivers.

If the rolling circle $K_2$ is on the inside of the base circle, the so-called hypotrochoid results.

Simple epitrochoids can be produced economically as inner and outer contours using the twin-spindle turning process [11]. Figure 5 shows some examples produced by [12] using this process.

Higher epitrochoids

The higher rolling curves were systematically treated for the first time by Wunderlich in [5] with the aid of complex functions. This allows the two parameter equations to be combined into one complex equation and reduces the mathematical effort. This procedure was also used in [4] for the representation of plane curves. According to Figure 6, the higher epitrochoid $K_m$ is defined by the planetary motion of several levers with the corresponding angular velocities ${\omega }_1,{\omega }_2,\cdot \cdot \cdot {\omega }_m$. The position of the creation point $A_m$ is determined by the sum of the vectors $A_0{A}_1$, $A_1{A}_2$, … $A_{m-1}{A}_m$ and can be determined as follows.

The following conditions apply to the higher epitrochoids [3]:

Angular velocities, ${\omega }_j$, all act in the same direction and in the positive direction.

$$r>0,{\omega }_1=n+1,{\omega }_1=2n+1,\cdot \cdot \cdot {\omega }_m=m\cdot n+1,$$

where $n\ge 1$ is arbitrary and an integer.

The higher trochoids are described in detail in [3]. The profiles generally have the following parameter descriptions:

$$x=r\mathrm{c}\mathrm{o}\mathrm{s}\left(t\right)+e_1\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(n+1\right)t\right)+e_2\mathrm{c}\mathrm{o}\mathrm{s}\left(\left(2n+1\right)t\right)+\dots +e_m\cos \left(\left(m·n+1\right)t\right)$$

$$y=r\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right)+e_1\mathrm{s}\mathrm{i}\mathrm{n}\left(\left(n+1\right)t\right)+e_2\mathrm{s}\mathrm{i}\mathrm{n}\left(\left(2n+1\right)t\right)+\dots +e_m\sin \left(\left(m·n+1\right)t\right)$$

(5)

In a complex form, Equation (3) can be combined into a single equation:

$$z=r\cdot {\mathrm{e}}^{i{\omega }_{0}t}+e_1\cdot {\mathrm{e}}^{i{\omega }_{1}t}+e_2\cdot {\mathrm{e}}^{i{\omega }_{2}t}+\cdot \cdot \cdot =r\cdot {\mathrm{e}}^{i{\omega }_{0}t}+{\sum }_{j=1}^m{e}_j\cdot {\mathrm{e}}^{i{\omega }_{j}t}$$

(6)

If $\zeta ={\mathrm{e}}^{it}$ and the associated ${\omega }_j$ are inserted into Equation (6), the following equation can be obtained to describe the higher epitrochoid:

$$\omega \left(\zeta \right)=r\zeta +e_1{\zeta }^{n+1}+e_2{\zeta }^{2n+1}+e_3{\zeta }^{3n+1}+e_4{\zeta }^{4n+1}\cdot \cdot \cdot .$$

(7)

Figure 7 shows examples of a class of higher epitrochoids with four eccentricities from the generalized Equation (7) for different numbers of flanks, using the following values for nominal radii and eccentricities:

$$r=20,e_0={\displaystyle \frac{r}{{\left(n+1\right)}^2}},$$

$$e_1={\displaystyle \frac{3n}{4}}{e}_0,e_2={\displaystyle \frac{1}{3}}{e}_1,e_3={\displaystyle \frac{1}{3}}{e}_2\mathrm{and}{e}_4={\displaystyle \frac{1}{3}}{e}_3.$$

(8)

Figure 8 compares a simple epitrochoid with one eccentricity and a higher epitrochoid with four eccentricities, according to Equations (7) and (8). The figure shows that the higher epitrochoid (on the right) has a larger radius of curvature when the eccentricity $e_0$ is increased. This results in advantages from a manufacturing point of view, as well as a smaller shape number (stress factor) in terms of the profile when compared to a simple epitrochoid (left in Figure 8).

2.1. State of the Art

Torsional loading

A torsional function was found using a suitable mapping that transforms the area enclosed by the unit circle to conform to the profile base area, and the torsional stresses and deformations were determined based on this.

$$f\left(\zeta \right)=\phi \left(\zeta \right)+i\psi \left(\zeta \right),$$

(9)

Here, $\phi$ and $\psi$ are conjugate to each other. The real part of the torsion function $\phi$ can be used to determine the torsional moment of inertia and the stress components. The method is discussed in detail in [2].

Pure bending load

For the case of pure bending, it is assumed that the plane cross-sections remain plane after loading (see, e.g., [13]). This means the following for the stresses:

$${\sigma }_z=-{\displaystyle \frac{M_b}{I_y}}x,{\sigma }_x={\sigma }_y={\tau }_{xy}={\tau }_{zx}={\tau }_{zy}=0.$$

(10)

Based on this, the task is reduced to determining the bending moment of inertia, $I_y$, of the profile cross-section. The procedure is described in detail in [3], and solutions for the M-sections were derived.

The method was used in [14], specifically and in detail, for the hypotrochoidal profiles standardized according to DIN 3689, and the so-called stress factors were determined and presented in tabular form.

2.2. Shear force Bending and Formulation

The general formulation for shear force bending is presented in [13,15] based on the elasticity behavior of the material. The method is described in more detail in [15], but a corresponding practical application is complex, as four stress functions must first be determined as a function of the profile contour. The method is, therefore, no longer suitable for practical applications. In contrast, the procedure described in [13] is presented in a more compact form. However, this also requires a complex bending stress as a function of the profile geometry. Based on the latter formulation, a general procedure was developed in the present work from a parametric profile description, according to Equations (3) and (4). The application was demonstrated in more detail using a number of test examples.

In order to analyze the transverse shear stresses, the origin of the co-ordinate is placed in the center of gravity of the profile surface, as shown in Figure 9. Force F is applied on the right side at the center of the cross-section and in the $x$-direction.

At any axial point at $z$ distance from the co-ordinate origin $O$, the resulting bending moment $M_b$ in relation to the $y$-axis is as follows:

$$M_b=F\left(l-z\right),$$

(11)

where $l$ is the length of the beam. However, in the case of shear force bending, ${\tau }_{xz}$ and ${\tau }_{yz}$ no longer disappear, as the shear force causes a shear stress on the cross-section. The following equations are, hence, valid:

$${\sigma }_z=-{\displaystyle \frac{F\left(l-z\right)}{I_y}}x,{\sigma }_x={\sigma }_y={\tau }_{xy}=0.$$

(12)

If we insert these values into the equilibrium conditions, we obtain

$${\displaystyle \frac{{\partial \tau }_{xz}}{\partial z}}={\displaystyle \frac{{\partial \tau }_{yz}}{\partial z}}={\displaystyle \frac{{\partial \tau }_{xz}}{\partial x}}+{\displaystyle \frac{\partial }{\partial y}}\left({\partial \tau }_{yz}+{\displaystyle \frac{Fx}{I_y}}\right)=0$$

(13)

The first two parts of the equation mean that ${\tau }_{xz}$ and ${\tau }_{yz}$ are not dependent on $z$.

For the further procedure, a complex stress function for bending $G\left(z\right)$ is required, which satisfies the corresponding compatibility conditions of the deformations, as well as the boundary conditions (load-free profile contour). It has the following general form:

$$G\left(z\right)={\Phi }_1\left(x,y\right)+i{\Psi }_1\left(x,y\right),$$

(14)

where ${\Phi }_1\left(x,y\right)$ describes the real and ${\Psi }_1\left(x,y\right)$ the imaginary part of the $G$-function. Both functions are dependent on the contour geometry of the profile (e.g., Equations (1) and (3)). If the variable of $z$ is transformed into ζ, Equation (5) can be presented as follows:

$$g\left(\zeta \right)={\phi }_1\left(\zeta \right)+i{\psi }_1\left(\zeta \right),$$

(15)

Here, ${\phi }_1\left(\zeta \right)$ and ${\psi }_1\left(\zeta \right)$ are conjugate to each other, i.e.,

$${\displaystyle \frac{{\partial \phi }_1}{\partial x}}={\displaystyle \frac{{\partial \psi }_1}{\partial y}},{\displaystyle \frac{{\partial \phi }_1}{\partial y}}=-{\displaystyle \frac{{\partial \psi }_1}{\partial x}}$$

(16)

${\phi }_1$ is labeled as a bending stress function. If a conformal mapping for transforming the unit circle on the cross-sectional area is available (Figure 10), a solution of ${\phi }_1$ for the entire profile cross-section can be determined from the boundary condition for ${\psi }_1$ on the profile contour. The corresponding procedure is described in more detail below using an example.

The boundary condition on the profile contour, ${\psi }_1$ can be derived as follows [13]:

$${\psi }_1=-\left(1-{\displaystyle \frac{\upsilon }{2}}\right){\displaystyle \frac{y^3}{3}}-\upsilon {\displaystyle \frac{x^{2}y}{2}}+2\left(1+\upsilon \right)\int yxdx+C_1$$

(17)

Due to the closed profile contour, the constant value $C_1$ on the right side of Equation (8) disappears. If the ${\psi }_1$ function is known, $g\left(\zeta \right)$ or ${\phi }_1$ can be determined using Cauchy's integral formula.

Procedure for determining ${\mathit{\phi }}_1$

If the contour geometry has been described by parametric equations with “cos” and “sin” terms (e.g., Equation (1) or (3)), the cos and sin functions are first rewritten as expressions of $\lambda ={\mathrm{e}}^{it}$:

$$\mathrm{c}\mathrm{o}\mathrm{s}\left(t\right)={\displaystyle \frac{1}{2}}\left({\mathrm{e}}^{it}+{\mathrm{e}}^{-it}\right)={\displaystyle \frac{1}{2}}\left(\lambda +{\displaystyle \frac{1}{\lambda }}\right)$$

$$\mathrm{s}\mathrm{i}\mathrm{n}\left(t\right)={\displaystyle \frac{1}{2i}}\left({\mathrm{e}}^{it}-{\mathrm{e}}^{-it}\right)={\displaystyle \frac{1}{2i}}\left(\lambda -{\displaystyle \frac{1}{\lambda }}\right)$$

(18)

where e is the Euler number. If these expressions are inserted into the parametric equations, the variables $x$ and $y$ can be represented as functions of $\lambda$. $\lambda$ is the edge of the unit circle and corresponds to the profile contour in the real plane. Based on this, Equation (8) is represented as follows:

$${\psi }_1\left(\lambda \right)=-\left(1-{\displaystyle \frac{\upsilon }{2}}\right){\displaystyle \frac{{y\left(\lambda \right)}^3}{3}}-\upsilon {\displaystyle \frac{{x\left(\lambda \right)}^{2}y\left(\lambda \right)}{2}}+2\left(1+\upsilon \right)\int x\left(\lambda \right)y\left(\lambda \right)x´\left(\lambda \right)d\lambda +C_1$$

(19)

This can be used to determine $g\left(\zeta \right)$ or ${\phi }_1$using Cauchy's integral formula:

$$g\left(\zeta \right)={\displaystyle \frac{1}{\pi i}}\oint {\displaystyle \frac{{\psi }_1\left(\lambda \right)}{\lambda -\zeta }}d\lambda +const$$

(20)

${\phi }_1$and ${\psi }_1$ are the real and imaginary parts of the complex bending stress function in the so-called model or ζ plane. The stress components can be determined using Equations (4) and (5) as follows:

$${\tau }_{xz}=-{\displaystyle \frac{F}{2\left(1+\upsilon \right)I_y}}\left[{\displaystyle \frac{\partial {\phi }_1}{\partial x}}+{\displaystyle \frac{1}{2}}\upsilon x^2+\left(1-{\displaystyle \frac{\upsilon }{2}}\right)y^2\right]$$

$${\tau }_{yz}=-{\displaystyle \frac{F}{2\left(1+\upsilon \right)I_y}}\left[{\displaystyle \frac{\partial {\phi }_1}{\partial y}}+\left(2+\upsilon \right)xy\right]$$

(21)

$${\sigma }_z=-{\displaystyle \frac{F\left(l-z\right)}{I_y}}x$$

For the partial derivatives of ${\phi }_1$occurring in Equation (21), ${\psi }_1$ is also required:

$${\displaystyle \frac{\partial {\phi }_1}{\partial x}}={\displaystyle \frac{\dot{x}{\dot{\phi }}_1+\dot{y}{\dot{\psi }}_1}{{\dot{x}}^2+{\dot{y}}^2}},$$

$${\displaystyle \frac{\partial {\phi }_1}{\partial y}}={\displaystyle \frac{\dot{y}{\dot{\phi }}_1-\dot{x}{\dot{\psi }}_1}{{\dot{x}}^2+{\dot{y}}^2}}$$

(22)

where $\dot{x},\dot{y},{\dot{\phi }}_1$, and ${\dot{\psi }}_1$ denote the corresponding derivatives according to parameter $t$.

$I_y$ is the moment of inertia around the $y$-axis and can be determined from the following relationship according to [3]:

$$I_y=\iint x^{2}dxdy$$

(23)

The deformations can then also be determined as follows [13]:

$$u_x={\displaystyle \frac{F}{E{I}_y}}\left[{\displaystyle \frac{\upsilon }{2}}\left(l-z\right)\left(x^2-y^2\right)+{\displaystyle \frac{1}{2}}l{z}^2-{\displaystyle \frac{1}{6}}{z}^3\right]$$

$$u_y={\displaystyle \frac{F}{E{I}_y}}\upsilon \left(l-z\right)xy$$

(24)

$$u_z=-{\displaystyle \frac{F}{E{I}_y}}\left[x\left(lz-{{\displaystyle \frac{1}{2}}z}^2\right)+{\phi }_1-x{y}^2\right]$$

The displacement in the z-direction $u_z$ (the 3rd equation in (24)) consists of a $z$-dependent component, $u_{z1}$, which represents the rotation of the cross-section (Equation (25)), and a warping, $u_{z2}$, which is caused by the distribution of the shear stress. $u_{z2}$ is independent of $z$ and can be determined from Equation (26):

$$u_{z1}=-{\displaystyle \frac{F}{E{I}_y}}\left[x\left(lz-{{\displaystyle \frac{1}{2}}z}^2\right)\right],$$

(25)

$$u_{z2}=-{\displaystyle \frac{F}{E{I}_y}}\left({\phi }_1-x{y}^2\right)$$

(26)

The deflection is obtained from the first equation in (24), where $x$ and $y$ are equal to zero:

$$f_x={\displaystyle \frac{F}{2E{I}_y}}\left(l{z}^2-{\displaystyle \frac{1}{3}}{z}^3\right)$$

(27)

In the following, the procedures are first applied to the circular cross-section as a well-known classic example, and the results are compared with the solutions known from the technical literature. Then, the epitrochoidal profile cross-sections are also investigated.

3. Applications

3.1. Circular Cross-Section

The following conformal mapping applies to a circular cross-section:

$$z\left(\zeta \right)=r\zeta ,$$

(28)

This function forms the unit circle conformally and completely on the cross-section of the shaft. The contour edge can be described as parametric:

$$x=r\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{t}\right)$$

$$y=r\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{t}\right)$$

(29)

If Equation (18) is inserted into (29) and the results for $x$ and $y$ are inserted into Equation (19), the function ${\psi }_1$ for the contour can be obtained. Based on this, the complex stress function $g\left(\zeta \right)$ is determined from Equation (20). If $\zeta =\rho {\mathrm{e}}^{it}$ is inserted into $g\left(\zeta \right)$, the real part can be determined as follows:

$${\phi }_1={\displaystyle \frac{1}{4}}\rho r^3\cos \left(t\right)\left(-2\nu -{\rho }^2+2{\rho }^2\cos \left(2t\right)-3\right)$$

(30)

Here, $\rho$ assumes values between 0 and 1. Substitute ${\phi }_1$ from Equation (30) into Equation (21) and replace $x$ and $y$ with the following relationships:

$$x=r\rho \cos \left(t\right)={\displaystyle \frac{r\rho }{2}}\left(\lambda +{\displaystyle \frac{1}{\lambda }}\right)$$

$$y=r\rho \sin \left(t\right)={\displaystyle \frac{r\rho }{2i}}\left(\lambda -{\displaystyle \frac{1}{\lambda }}\right),$$

(31)

Thus, the equations for the shear stress components are obtained:

$${\tau }_{xz}=-{\displaystyle \frac{F\left({\rho }^2\left(2+\left(2\nu +1\right)\cos \left(2t\right)\right)-2\nu -3\right)}{2\pi r^2\left(\nu +1\right)}}$$

(32)

$${\tau }_{yz}=-{\displaystyle \frac{F\left(2\nu +1\right){\rho }^2\sin \left(2t\right)}{2\pi r^2\left(\nu +1\right)}}$$

(33)

The resulting shear stress can then be determined as the vectorial sum of the two components:

$$\tau =\sqrt{{\tau }_{xz}^2+{\tau }_{yz}^2}$$

(34)

With $A=\pi r^2$, the shear stress components can be determined as follows:

$${\tau }_{xz}=-{\displaystyle \frac{F}{A}}·{\displaystyle \frac{\left({\rho }^2\left(2+\left(2\nu +1\right)\cos \left(2t\right)\right)-2\nu -3\right)}{2\left(\nu +1\right)}}$$

(35)

$${\tau }_{yz}=-{\displaystyle \frac{F}{A}}·{\displaystyle \frac{\left(2\nu +1\right){\rho }^2\sin \left(2t\right)}{2\left(\nu +1\right)}}$$

(36)

The stress distribution on the cross-section is shown in Figure 11. This is the result of the shear stress component ${\tau }_{xz}$ from the numerical investigations (FEA) for a round cross-section with $r=25$ mm and $F=100kN$. The material chosen for the shaft was steel with $E=210000N/{mm}^2$ and a Poisson’s ratio of $\nu =0.3$. According to Equation (32), the analytical solutions for the two parameter angle positions $t=0$ and $t={\displaystyle \frac{\pi }{2}}$ are also shown, where the agreement of the results can be traced by observing the color scalar.

A more concrete comparison can be made using the stress distributions on the lateral surface, where $\rho =1$ applies. If $\rho =1$ is inserted into Equations (32) to (34), the following relationships are obtained:

$${\tau }_{xz}\left(1,t\right)={\displaystyle \frac{F\left(2\nu +1\right){\mathrm{s}\mathrm{i}\mathrm{n}}^2\left(t\right)}{\pi r^2\left(\nu +1\right)}}$$

(37)

$${\tau }_{yz}\left(1,t\right)=-{\displaystyle \frac{F\left(2\nu +1\right)\sin \left(t\right)\cos \left(t\right)}{\pi r^2\left(\nu +1\right)}}$$

(38)

The resulting shear stress can be determined from (34) as follows:

$$\tau \left(1,t\right)={\displaystyle \frac{F\left(2\nu +1\right)\sin \left(\mathrm{t}\right)}{\pi r^2\left(1+\nu \right)}}$$

(39)

Figure 12 shows the course of ${\tau }_{xz}$ on the contour of the cross-section. A very good agreement between the numerical results and the analytical solution—according to Equation (37)—can be seen from the figure.

The maximum stress ${\tau }_{xz,\mathrm{m}\mathrm{a}\mathrm{x}}$of the contour occurs at $t={\displaystyle \frac{\pi }{2}}$ and can be determined according to the following equation:

$${{\tau }_{xz,\mathrm{m}\mathrm{a}\mathrm{x}}=\tau }_{xz}\left(1,{\displaystyle \frac{\pi }{2}}\right)={\displaystyle \frac{F\left(2\nu +1\right)}{\pi r^2\left(\nu +1\right)}}={\displaystyle \frac{F}{A}}·{\displaystyle \frac{2\nu +1}{\nu +1}}$$

(40)

The maximum shear stress occurs at the center of the cross-section (0,0). At $\rho =0$, ${\tau }_{xz}$ is maximum and ${\tau }_{yz}$=0; therefore,

$${{\tau }_{max}=\tau }_{xz}\left(\mathrm{0,0}\right)={\displaystyle \frac{F\left(2\nu +3\right)}{2\pi \left(\nu +1\right)r^2}}={\displaystyle \frac{F}{A}}·{\displaystyle \frac{\left(2\nu +3\right)}{2\left(\nu +1\right)}}$$

(41)

For a Poisson’s ratio of $\nu =0.3$ (steel), the maximum shear stress at the center point can be determined from Equation (41), and the maximum shear stress on the profile outline from Equation (40) is as follows:

$${\tau }_{max}=1.385{\displaystyle \frac{F}{A}}\left(\mathrm{at}\mathrm{the}\mathrm{center}\right),$$

(40a)

$${\tau }_{max}=1.231{\displaystyle \frac{F}{A}}\left(\mathrm{on}\mathrm{the}\mathrm{profile}\mathrm{outline}\right).$$

(41a)

Compared with existing solutions

The solution for circular cross-sections known from the literature is based on the formulation of the classical theory of elasticity, which always depends on finding a suitable real shear function. This shear function $\varphi (x,y)$ should satisfy the corresponding equilibrium and compatibility conditions depending on the profile geometry. Therefore, solutions are only available for some specific profile geometries such as a circle, ellipse, or quadrilateral. Timoshenko [16] presents such a solution for the circular cross-section. Here, a real stress function is given in Cartesian co-ordinates, which leads to the same results from Equations (32) and (33):

$$\varphi \left(x,y\right)={\displaystyle \frac{\left(1+2\nu \right)}{8\left(1+\nu \right)}}{\displaystyle \frac{F}{I_y}}\left(x^2+y^2-{rR}^2\right)y,$$

(42)

$R$ is here the radius of the cross-section. The same results from (40a) and (41a) are obtained in [16].

3.2. Epitrochoidal Profile Cross-Sections

For the epitrochoidal cross-section, the conformal mapping, according to Equation (2), applies, where $n$ is the number of sides of the profile. This function maps the unit circle conformally and completely on the profile cross-section (see Figure 10). This guarantees the solution over the entire profile cross-section according to Equation (20).

The contour in Equation (1) can be rewritten as follows:

$$x={\displaystyle \frac{r}{2}}\left(\lambda +{\displaystyle \frac{1}{\lambda }}\right)+{\displaystyle \frac{e}{2}}\left({\lambda }^{n+1}+{\displaystyle \frac{1}{{\lambda }^{n+1}}}\right)$$

$$y={\displaystyle \frac{r}{2i}}\left(\lambda -{\displaystyle \frac{1}{\lambda }}\right)+{\displaystyle \frac{e}{2i}}\left({\lambda }^{n+1}-{\displaystyle \frac{1}{{\lambda }^{n+1}}}\right),$$

(43)

Substituting Equation (43) into Equation (19) gives the function ${\psi }_1$ for the profile contour. Based on this, the complex stress function $g\left(\zeta \right)$ can be determined from Equation (20) using Cauchy's integral formula (20). If $\zeta =\rho {\mathrm{e}}^{it}$ is inserted into $g\left(\zeta \right)$, the real part of the stress function ${\phi }_1$ can be determined. Substituting ${\phi }_1$ into Equation (21) and replacing $x$and $y$ with the following relations for the entire profile cross-section yields

$$x={\displaystyle \frac{\rho r}{2}}\left(\lambda +{\displaystyle \frac{1}{\lambda }}\right)+{\displaystyle \frac{\rho e}{2}}\left({\lambda }^{n+1}+{\displaystyle \frac{1}{{\lambda }^{n+1}}}\right)$$

$$y={\displaystyle \frac{\rho r}{2i}}\left(\lambda -{\displaystyle \frac{1}{\lambda }}\right)+{\displaystyle \frac{\rho e}{2i}}\left({\lambda }^{n+1}-{\displaystyle \frac{1}{{\lambda }^{n+1}}}\right),$$

(44)

The components of the shear stresses can be determined for the entire profile cross-section. Here, $\rho$ takes values between 0 and 1 ($0\le \rho \le 1$).

Because of the size of the equation, the general solution is not shown here. This can be found in the appendix.

In the following, two application examples for $n=3$ (Figure 5 left) and n = 5 (Figure 5 right) are presented, and the results are compared with the FE analysis.

Example 1: $\mathit{n}=3$

with flat flanks

For an epitrochoid with three sides, the following mapping function applies, which conformally transforms the unit circle on the epitrochoid:

$$z\left(\zeta \right)=r\zeta +e{\zeta }^4$$

(45)

If e=r/16 is inserted into Equation (44), the epitrochoid with three flat point flanks is obtained:

$$z\left(\zeta \right)=r\left(\zeta +{\displaystyle \frac{{\zeta }^4}{16}}\right).$$

(46)

Figure 13 shows the iso-$\rho$ and $t$ lines for this profile with $e={\displaystyle \frac{r}{16}}$ .

If $n=3$ and $e={\displaystyle \frac{r}{16}}$ are inserted into Equation (A4) in the appendix, the following relationship is obtained for the stress function:

$${\phi }_1={\displaystyle \frac{1}{114688}}\rho \left(\rho \left(\rho \left(28672r^3\cos \left(3t\right)-7\rho \left(2\nu +256\left(\nu +3\right)r^2+3\right)\cos \left(4t\right)+16{\rho }^{3}r\left(336r\cos \left(6t\right)-\left(2\nu +9\right)\rho \cos \left(7t\right)+21{\rho }^3\cos \left(9t\right)\right)+7{\rho }^9\cos \left(12t\right)\right)-1792\left(4\nu +3\right)r^2\cos \left(2t\right)\right)-224r\left(8\nu +128\left(2\nu +3\right)r^2+9\right)\cos \left(t\right)\right)$$

(47)

If $\nu ={\displaystyle \frac{3}{10}}$ is also used for steel, the following function for the distribution of the shear stress component ${\tau }_{xz}\left(\rho ,t\right)$ on the profile cross-section is obtained from (A6):

$${\tau }_{xz}\left(\rho ,t\right)=-{\displaystyle \frac{F}{221325\pi r^2({\rho }^6+8{\rho }^3\mathrm{c}\mathrm{o}\mathrm{s}(3t)+16)}}128(5{\rho }^{14}+2000{\rho }^8-2121{\rho }^6+20480{\rho }^2+4{\rho }^{14}\mathrm{c}\mathrm{o}\mathrm{s}(8t)+200{\rho }^{11}\mathrm{c}\mathrm{o}\mathrm{s}(3t)+144{\rho }^{11}\mathrm{c}\mathrm{o}\mathrm{s}(5t)+16{\rho }^{11}\mathrm{c}\mathrm{o}\mathrm{s}(11t)-96{\rho }^9\mathrm{c}\mathrm{o}\mathrm{s}(3t)+640{\rho }^8\mathrm{c}\mathrm{o}\mathrm{s}(6t)+576{\rho }^8\mathrm{c}\mathrm{o}\mathrm{s}(8t)-384{\rho }^6\mathrm{c}\mathrm{o}\mathrm{s}(6t)+12800{\rho }^5\mathrm{c}\mathrm{o}\mathrm{s}(3t)+6144{\rho }^5\mathrm{c}\mathrm{o}\mathrm{s}(5t)+256(16{\rho }^4-21)\rho \mathrm{c}\mathrm{o}\mathrm{s}(t)-17928{\rho }^3\mathrm{c}\mathrm{o}\mathrm{s}(3t)+64(24{\rho }^6-21{\rho }^2+256){\rho }^2\mathrm{c}\mathrm{o}\mathrm{s}(2t)-37776)$$

(48)

Furthermore, the following equation for the shear stress component ${\tau }_{yz}\left(\rho ,t\right)$ results from (A7):

$${\tau }_{yz}\left(\rho ,t\right)=-{\displaystyle \frac{F\rho }{221325\pi r^2\left({\rho }^6+8{\rho }^3\cos \left(3t\right)+16\right)}}1024\sin \left(t\right)\left(20{\rho }^{10}+12{\rho }^8+256{\rho }^4-120{\rho }^2+{\rho }^{13}\cos \left(3t\right)+{\rho }^{13}\cos \left(5t\right)+{\rho }^{13}\cos \left(7t\right)+40{\rho }^{10}\cos \left(4t\right)+4{\rho }^{10}\cos \left(6t\right)+4{\rho }^{10}\cos \left(8t\right)+4{\rho }^{10}\cos \left(10t\right)+144{\rho }^7\cos \left(3t\right)+144{\rho }^7\cos \left(5t\right)+144{\rho }^7\cos \left(7t\right)+96{\rho }^5\cos \left(3t\right)+96{\rho }^5\cos \left(5t\right)+1536{\rho }^4\cos \left(4t\right)+8\left(5{\rho }^8+3{\rho }^6+192{\rho }^2-30\right){\rho }^2\cos \left(2t\right)+\left({\rho }^{12}+528{\rho }^6+96{\rho }^4-336{\rho }^2+4096\right)\rho \cos \left(t\right)+672\right)$$

(49)

Figure 14 shows the distribution of the shear stress ${\tau }_{xz}$ from the numerical FE investigations for the epitrochoidal cross-section with $n=3$, $r=16$ mm, $e=1$ mm, and $F=100$ kN (for a load). The material chosen for the bar was steel with $E=210000N/m{m}^2$ and a Poisson’s ratio of $\nu ={\displaystyle \frac{3}{10}}$. Analytical solutions according to Equation (48) for the two parameter angle positions of $t=0$ and $t={\displaystyle \frac{\pi }{2}}$ are also shown, where the agreement of the results can be traced by observing the color scalar.

In order to determine the resulting shear stresses on the profile outline, $\rho =1$ is used in Equations (48) and (49). Equation (34) is used to obtain the following relationship for the total stress distribution on the contour of the shaft:

$$\tau \left(1,t\right)={\displaystyle \frac{1024F\left(1056\sin \left(t\right)+80\sin \left(2t\right)+129\sin \left(4t\right)+4\sin \left(7t\right)\right)}{221325\pi r^2\sqrt{8\cos \left(3t\right)+17}}}$$

(50)

A comparison of Equation (50) with the numerical results is shown in Figure 15, where very good agreement can be recognized.

Maximum shear stress occurs at the center point of the cross-section $\left(\mathrm{0,0}\right)$:

$${\tau }_{max}=\tau \left(\mathrm{0,0}\right)={\displaystyle \frac{64\left(88\nu +131\right)}{5675\left(\nu +1\right)}}·{\displaystyle \frac{F}{{\pi r}^2}}$$

(51)

Example 2: $\mathit{n}=5$

For an epitrochoid with five sides, the following mapping function applies, which can conformally transform the unit circle on the epitrochoid area:

$$z\left(\zeta \right)=r\zeta +e{\zeta }^6$$

(52)

Figure 16 shows this profile for this example with $e=r/20$. The iso-$t$ and $\rho$ lines are shown in the figure.

If $n=5$ and $e=r/20$ are inserted into Equation (A4), the following relationship is obtained for the stress function:

$${\phi }_1={\displaystyle \frac{1}{132}}\left({\rho }^3\left(e\rho \left(-11{\rho }^2\cos \left(6t\right)\left(e^2\left(6\nu +9\right)+2\left(\nu +4\right)r^2\right)+33e^2{\rho }^{14}\cos \left(18t\right)-6e\nu {\rho }^{7}r\cos \left(11t\right)+99e{\rho }^{9}r\cos \left(13t\right)-39e{\rho }^{7}r\cos \left(11t\right)-33-33\left(3\nu +2\right)r^2\cos \left(4t\right)+99{\rho }^4{r}^2\cos \left(8t\right)\right)+33r^3\cos \left(3t\right)\right)-33\rho r\cos \left(t\right)\left(e^2\left(24\nu +26\right)+\left(2\nu +3\right)r^2\right)\right)$$

(53)

Furthermore, Equation (A5) yields the following relationship for ${\psi }_1$:

$${\psi }_1={\displaystyle \frac{1}{132}}\left({\rho }^3\left(e\rho \left(-11{\rho }^2\sin \left(6t\right)\left(e^2\left(6\nu +9\right)+2\left(\nu +4\right)r^2\right)+33e^2{\rho }^{14}\sin \left(18t\right)-6e\nu {\rho }^{7}r\sin \left(11t\right)+99e{\rho }^{9}r\sin \left(13t\right)-39e{\rho }^{7}r\sin \left(11t\right)-33\left(3\nu +2\right)r^2\sin \left(4t\right)+99{\rho }^4{r}^2\sin \left(8t\right)\right)+33r^3\sin \left(3t\right)\right)-33\rho r\sin \left(t\right)\left(e^2\left(24\nu +26\right)+\left(2\nu +3\right)r^2\right)\right)$$

(54)

Consequently, the following equation for the shear stress component ${\tau }_{xz}(\rho ,t)$ on the profile cross-section can be obtained for $\nu =3/10$ (steel) from (A6):

$${\tau }_{xz}\left(\rho ,t\right)={\displaystyle \frac{F}{26I_y\left(9{\rho }^{10}+60{\rho }^5\cos \left(5t\right)+100\right)}}\left(-45{\rho }^{22}-24500{\rho }^{12}+25881{\rho }^{10}-200000{\rho }^2+4{\rho }^2\left(-30{\rho }^{15}\cos \left(17t\right)-\left(9{\rho }^{10}+1300\right){\rho }^{10}\cos \left(12t\right)-500\left(24{\rho }^4-29\right)\rho \cos \left(3t\right)-50\left(96{\rho }^{10}-87{\rho }^6+800\right)\cos \left(2t\right)-10\left(39{\rho }^{10}+1600\right){\rho }^5\cos \left(7t\right)+50\left(17-30{\rho }^2\right){\rho }^8\cos \left(10t\right)+5\left(-105{\rho }^{12}+51{\rho }^{10}-7000{\rho }^2+9838\right){\rho }^3\cos \left(5t\right)\right)+368300\right)$$

(55)

Furthermore, the following relation follows from (A7) for ${\tau }_{yz}\left(\rho ,t\right)$:

$${\tau }_{yz}\left(\rho ,t\right)={\displaystyle \frac{2F}{13I_y\left(9{\rho }^{10}+60{\rho }^5\cos \left(5t\right)+100\right)}}{\rho }^2\left(-30{\rho }^{15}\sin \left(17t\right)-\left(9{\rho }^{10}+1300\right){\rho }^{10}\sin \left(12t\right)-850{\rho }^8\sin \left(10t\right)+500\left(24{\rho }^4-29\right)\rho \sin \left(3t\right)-50\left(96{\rho }^{10}-87{\rho }^6+800\right)\sin \left(2t\right)-10\left(39{\rho }^{10}+1600\right){\rho }^5\sin \left(7t\right)+5\left(1211-51{\rho }^{10}\right){\rho }^3\sin \left(5t\right)\right)$$

(56)

Figure 17 shows the distribution of the shear stress components ${\tau }_{xz}$ from the numerical FE investigations for the epitrochoidal cross-section with $n=5$, $r=20$mm, and $e=1$ mm. A shear force of $F=100kN$ was applied, and the chosen material was steel with $E=210000N/{mm}^2$ and a Poisson’s ratio of $\nu =0.3$. Analytical solutions according to Equation (55) for the three parameter angle positions $t=0$, $t=\pi /2$, and $t=3\pi /5$ are also shown, where the agreement of the results can be identified using the color scalar.

The total shear stress is the vectorial sum of the two components according to Equations (55) and (56). Figure 18 shows the distribution of the resulting shear stresses over different $\rho$ -lines.

The shear stresses on the profile contour correspond to $\rho =1$. This fact leads to the following equation for the distribution of the shear stress on the profile contour with $n=5$:

$$\tau \left(1,t\right)={\displaystyle \frac{F\left(\left(2\nu +1\right)r\left(12e^2+r^2\right)\sin \left(t\right)+e\left(2\nu +1\right)\left(2\left(3e^2+r^2\right)\sin \left(6t\right)+\mathit{er}\sin \left(11t\right)\right)+2e{r}^2\sin \left(4t\right)\right)}{4I_y\left(\nu +1\right)\sqrt{36e^2+12er\cos \left(5t\right)+r^2}}}$$

(57)

The following then applies to the investigation case with $r=20$ and $e=1$:

$$\tau \left(1,t\right)={\displaystyle \frac{F\left(4120\left(2\nu +1\right)\sin \left(t\right)+\left(2\nu +1\right)\left(403\sin \left(6t\right)+10\sin \left(11t\right)\right)+400\sin \left(4t\right)\right)}{4I_y\left(\nu +1\right)\sqrt{60\cos \left(5t\right)+109}}}$$

(58)

Equation (58) represents the curve with $\rho =1$ in Figure 18.

Bending stress according to [3]

For the general case of bending stress, according to [3], the following relationship can be determined for the bending stress due to a shear force:

$${\sigma }_z=-{\displaystyle \frac{4F\rho \left(l-z\right)\left(e{\rho }^n\sin \left(nt+t+\varphi \right)+\mathrm{rsin}\left(t+\varphi \right)\right)}{\pi \left(e^4\left(n+1\right)+2e^2\left(n+2\right)r^2+r^4\right)}}$$

(59)

Here, $\varphi$ denotes the angle of rotation of the shear force, $F$, around the $z$-axis. Variable $\rho$ accepts the value 0 for the center point (neutral axis) and the value 1 for the outer profile contour. Figure 19 shows the distribution of the bending stress, according to Equation (59), for different $\rho$ -lines.

Deformations

The deformations of the profile cross-section can be determined from (24) as functions of $x$ and $y$ with the known stress function ${\phi }_1$. The results for the three displacement components in curvilinear co-ordinates $(\rho ,t)$ are summarized in the appendix. If $n=5$ is inserted into Equations (A9) to (A11), the following equations are obtained for the deformation components:

$$u_x={\displaystyle \frac{F}{6E{I}_y}}\left(3\nu {\rho }^2\left(l-z\right)\left(r^2\cos \left(2t\right)+2{\rho }^{5}r\cos \left(7t\right)+{\rho }^{10}\cos \left(12t\right)\right)+3l{z}^2-z^3\right)$$

(60)

$$u_y={\displaystyle \frac{F}{E{I}_y}}\nu {\rho }^2\left(l-z\right)\left(r\sin \left(t\right)+{\rho }^5\sin \left(6t\right)\right)\left(r\cos \left(t\right)+{\rho }^5\cos \left(6t\right)\right)$$

(61)

$$u_z={\displaystyle \frac{F}{132E{I}_y}}\rho \left(11{\rho }^5\cos \left(6t\right)\left(3\left(4lz-2\nu +{\rho }^{12}-2z^2-3\right)-2r^2\left(\nu -3{\rho }^2+4\right)\right)-33r\cos \left(t\right)\left(2z\left(z-2l\right)+24\nu -2{\rho }^{12}+r^2\left(2\nu -{\rho }^2+3\right)+26\right)+33{\rho }^3{r}^2\left(-3\nu +{\rho }^4-2\right)\cos \left(4t\right)+3{\rho }^{10}r\left(-2\nu +11{\rho }^2-13\right)\cos \left(11t\right)\right)$$

(62)

The deflection can then be determined from Equation (60) for the center of the cross-section (with $\rho =0$) as follows:

$$f_x=u_x\left(\rho =0\right)={\displaystyle \frac{F}{6E{I}_y}}\left(3l-z\right)z^2$$

(63)

Figure 20 compares Equation (63) with the results of the FEA, wherein a good match can be recognized.

4. Summary and Outlook

In this paper, first, the geometry of simple and higher epitrochoids, which are suitable for technical applications as positive shaft–hub connections, was depicted. An approach to determine the transverse shear stresses in non-circular elastic profile cross-sections was then presented. The formulation is based on the complex stress functions of Muskhelishvili [13]. The developed method was first applied to a round profile cross-section, and the results were compared to the solutions known from the literature, according to Timoshenko [16]. A complete agreement was proven. Furthermore, two epitrochoidal profiles with three and five flanks were examined. The analytical stress results show very good agreement with the numerical investigations from FEA. For the example with five flanks, the bending stresses were also determined and presented based on the approach described in Ziaei [3]. The elastic deformation components were also determined for both examples. An equation for determining the deflection was derived from $u_x$, for which very good agreement with the FEA was also observed.

For the so-called H-profiles (hypotrochoids), no suitable mapping functions can be determined directly from parametric equations. Therefore, appropriate mappings using the method of Kantorovich [18] are to be developed and used to determine the stress function in future work. The results may be of interest for the hypotrochoidal profiles standardized in DIN 3689 [17].

In addition, the shear stresses due to a shear force can be super-imposed onto the torsional stresses in the case of combined loading. This issue is to be investigated in further study.