Axisymmetric Twist and Bending

When using a 2D axisymmetric Solid Mechanics interface, the primary dependent variables are, by default, the displacements in the radial and axial directions, u and w. It is, however, possible to also include the displacement component in the circumferential direction, v, which opens up the possibility to model both twisting and bending deformations. The extension of a 2D model to solve for a true 3D displacement field is sometimes referred to as a 2.5D formulation.

The primary objective of this model is to show how to model twist and bending in 2D axisymmetry. A hollow axisymmetric shaft is subjected to different load cases, and the stress concentration factors are derived for each case. The analyses are performed both in 2D axisymmetry and 3D to demonstrate their equivalence. However, if the geometry allows it, modeling in 2D axisymmetry is generally preferable being computationally significantly cheaper as well as more convenient to set up.

Figure 1: Shaft geometry and relevant geometric properties.

Model Definition

The geometry of the hollow shaft is shown in Figure 1. The transition between the thicker and thinner sections is slightly smoothened by a fillet. The shaft is made of steel (see Table 1) and behaves linear elastically.

Table 1: Material Properties.
Property	Cable
Young’s modulus	200 GPa
Poisson's ratio	0.3
Density	7850 kg/m3

The wider end is fixed, while a load is applied at the narrow end. Three different load cases are examined:

Axial extension

Torsion

Bending

In the transition region between the small and large shaft diameters, local stress concentrations are expected due to the disturbance in the stress field. For the three load cases, stress concentration factors are evaluated numerically in this model. In general, the geometrical stress concentration factor is defined as the ratio between the actual maximum stress, σmax, and the nominal stress, σnom, in the region where the stress field is disturbed:

Here, the nominal stress, σnom, is the stress computed from elementary theories. The stress concentration factors are defined as follows for the three load cases:

Axial extension:

where

Torsion:

where

Bending:

where

Here, Fz, Mt, and Mb are the applied force in the z direction, the torsional moment, and the bending moment, respectively. Furthermore, A, J, and I are the cross-section area, the torsional constant, and the area moment of inertia for the thinner section of the hollow shaft.

In order to solve for torsion and bending in 2D axisymmetry, the standard formulation has to be augmented by including the displacement in the circumferential direction, v, as a new dependent variable. By default, this circumferential displacement is assumed to be zero, and the extended 2D axisymmetric formulation must be turned on explicitly in COMSOL Multiphysics. Adding a new dependent variable increases the size of the system matrices slightly.

In addition to adding the displacement field component in the φ direction, the displacement gradient, ∇u, is augmented to include the spatial variation of the circumferential displacement. The displacement gradient is a fundamental quantity used to define important variables such as strains and subsequently stresses. There can be alternative formulations of the displacement gradient depending on the chosen study type. In general, the gradient of the displacement field in a cylindrical coordinate system is defined as

Assuming that the displacement is constant in the circumferential direction, the derivatives with respect to φ become zero, and ∇u simplifies to

This formulation, which allows a constant displacement around the circumference, is available for all study types.

For frequency-domain and eigenfrequency studies, the displacement field is allowed to have a periodic variation in the φ direction. This approach is sometimes referred to as a circumferential mode extension. The following complex-valued ansatz for the displacement field is assumed:

Here, m is the mode number. This ansatz allows periodic displacements in the circumferential direction, and it has two noteworthy special cases. The case of m = 0 corresponds to a constant circumferential displacement, and m = 1 allows a bending-type displacement. With the new ansatz, the displacement gradient becomes

This formulation is only available for frequency-domain and eigenfrequency studies, and it is only valid for small strains.

The assumed complex-valued displacement field has implications for how an excitation load should be applied when working with a cylindrical coordinate system since the force vector is to be entered according to the base vectors of the cylindrical system. In the corresponding 3D case, the bending load is applied as a total force acting in the positive X direction of the global Cartesian coordinate system. From a parallelogram of forces it can be shown that the force FX can be written in terms of force components expressed in the cylindrical system.

The force entered in the load feature acts directly on the displacement field, hence the force vector is implicitly multiplied by

. Therefore, the load vector should be defined only as

which corresponds to a real-valued load FX acting in the X direction for m = 1.

Results and Discussion

The von Mises stress distribution is shown in Figure 2–Figure 4 for the three load cases, that is, axial extension, torsion, and bending. The 2D axisymmetric solution closely matches the 3D solution for all cases.

To find the highest stress value in the region of interest, the stress intensification factors are evaluated using a maximum operator (see the evaluation group Stress Concentration Factors under the node of the ). Again, the computationally leaner 2D axisymmetric solution shows good agreement with the full 3D result. Axial extension produces the highest stress concentration factor, followed by bending and torsion.

Figure 2: von Mises stress for axial load as computed in 2D axisymmetry (left) and 3D.

Figure 3: von Mises stress for torsional load as computed in 2D axisymmetry (left) and 3D.

Figure 4: von Mises stress for bending load as computed in 2D axisymmetry (left) and 3D.

It is worth highlighting that the number of degrees of freedom (DOFs) is roughly 25 times higher for the 3D model, which directly affects assembly and solver time. When investigating stress concentrations, it is often of interest to examine effects of geometrical changes such as varying fillet or diameter sizes in between different shaft sections. Such geometrical investigations can easily be performed in COMSOL Multiphysics using parametric sweeps. The geometry needs to be rebuilt and remeshed, and the model needs to be re-solved. This can be done more efficiently with only a 2D geometry.

Notes About the COMSOL Implementation

The model is set up and solved for twice in order to show the equivalence between the 2D axisymmetric and full 3D formulations.

Bending deformations in 2D axisymmetry can only be computed in frequency-domain or in eigenfrequency studies. Therefore, a frequency-domain study is added with an excitation frequency of 0 Hz. This corresponds to a stationary analysis since the applied loads will have no time-harmonic part.

Note also that all three load cases could have been solved for with a single frequency-domain study. In this case, the azimuthal mode number would have had to be defined as

m = if(group.lg3, 1, 0)

Here, group.lg3 is the weight variable of load group 3 (bending). The variable is defined in the section when adding load cases. The above expression for m sets the azimuthal mode number to 0 for the first and second load cases (that is, axial extension and torsion) and to 1 for the third load case (that is, bending).

Structural_Mechanics_Module/Verification_Examples/axisymmetric_twist_and_bending

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
r0	2[mm]	0.002 m	Inner radius
R1	12[mm]	0.012 m	Outer radius, thick section
R2	5[mm]	0.005 m	Outer radius thin section
L1	18[mm]	0.018 m	Length, thick section
L2	50[mm]	0.05 m	Length, thin section
I	(pi/4)*(R2^4 - r0^4)	4.7831E-10 m4	Area moment of inertia
J	2*I	9.5661E-10 m4	Torsion constant
A_load	pi*(R2^2 - r0^2)	6.5973E-5 m²	Area of applied load
F_load	10[MPa]*A_load	659.73 N	Total applied force
Mt	(2pi/3)(F_load/A_load)*(R2^3 - r0^3)	2.4504 N·m	Equivalent twisting moment
Mb	F_load*(L2 - L1)	21.112 N·m	Bending moment at transition