The moments of inertia in the XY coordinate system are

The angle needed to rotate the x-axis to the axis of the largest principal moment of inertia (x1) is denoted α. From the definition of Mohr’s circle, the angle is

When implemented using the atan2 function, the angle can be correctly evaluated for all rotations, and returns in the interval −π < α < π.

The following notation is used: . This is a shear stress in the x2 direction (acting on the plane with z as normal) caused by a unit shear force acting in the x1 direction.

The max shear stress factor is the ratio between the maximum shear stress in the cross section and the average shear stress. For a shear force in the x1 direction, the definition is

is the resulting shear stress from a unit load in the x1 direction. Similarly:

The shear correction factor is also computed. The shear correction factor is a multiplier which makes the strain energy from the average shear stress and shear strain in the cross section equal to the true shear energy in the cross section. The shear correction factor can be introduced through the concept shear area. The shear area is the reduced area which should replace the true area when computing the shear deformation of a beam. In terms of the shear correction factor it can be written as

where κ1 is the shear correction factor for a shear force in the x1 direction. Thus, for a shear flexible beam, the constitutive relation for the average shear is

Since T1 is a unit shear force, the shear correction factor can be computed as

The shear center (or, equivalently, the center of rotation) is the point around which the shear stresses from bending has no torque. In COMSOL it is represented as the distance from the center of gravity of the cross section in the principal axes coordinate system. The torque can be computed as

The two first equations simply state that the shear stresses are independent of z, whereas the third equation is the one on which to focus the interest. Assume that the shear stresses can be derived from a scalar stress potential, ψ, through:

Given the stress state from Equation 9-1, the only two nontrivial equations are

Integration of the first equation with respect to x1 and the second equation with respect to x2 gives

This is the same equation as Equation 9-2. It is thus possible to fulfill equilibrium, compatibility, and the constitutive relation with a single equation of Poisson type.

The proofs for the two components are analogous, so it is shown only for the x1 direction:

To compute the integral of the x1-derivative of ψ, a term containing the differential equation itself is added. This is a zero contribution, but it makes further simplifications possible.

Since the only part of the displacement that is relevant for the bending shear stresses is independent of the z-coordinate, the functions g1 and g2 can be considered as independent of z:

In the last transformation Green’s theorem is used. The uniqueness of the u2 displacement can be shown in the same way.

The uniqueness of the out-of-plane component of the displacement is shown in Equation 9-3:

In the last step of Equation 9-3, all integrals are zero since the coordinate system is located at the center of gravity of the section. This proves that all displacement components are unique.

When solving the problem, the shear stresses caused by a unit force in each of the two principal directions must be separated, so two separate problems are solved. For the force in the x1 direction it is formulated as

The corresponding problem for the x2 direction is

where ϕ is the stress function. For a singly connected region the boundary condition is ϕ = 0 along the whole boundary. Having solved this problem, the torsional rigidity can be computed as

In the case that there are internal holes in the section the situation is slightly more complex. The ϕ = 0 condition is now applicable only to the external boundary, whereas each boundary of an internal hole i needs a Dirichlet boundary condition:

where Hi is a constant to be determined. The constant value of the stress function fulfills the stress-free boundary conditions. There is also a compatible condition that must be fulfilled: the displacements must be single valued when going around each hole along its boundary Γi. This is trivially fulfilled for the in-plane displacements, but the out-of-plane displacement, w, generates the necessary equations to determine Hi:

is employed. This assumption implies that the origin of the coordinate system is at the center of rotation. This is true only for doubly symmetric sections, but adding a constant offset to the x and y coordinates does not contribute to the integral.

The gradient of ϕ depends linearly on the yet unknown variables Hi, the values of which can be solved by adding one equation,

where Ai is the area of hole i.

The warping function ω(x, y) describes the out-of-plane deformation related to torsion. It fulfills Laplace’s equation Δω = 0.

The offset by the shear center coordinates (ex, ey) is introduced because the torsion theory assumes that the coordinate system has its origin in the center of rotation (which is the same as the shear center).

The warping displacement, u, for a beam with unrestrained warping can be computed by providing the twist of the beam, ϕxl as input:

where B is the bimoment. The maximum axial stress is

Given the warping constant, it is possible to compute a nondimensional number that can be used to characterize the influence of nonuniform torsion in a beam with a certain length L. This number is