When I was a kid, long before Spotify, iTunes and even Winamp, people used to go to physical stores to buy compact discs (a.k.a. CDs™). These devices proved to be a significant improvement over the previous ways of distributing music, namely cassettes and vinyl discs in a wide variety of ways, all of them being far outside the scope of this article. In any case, the CD player would lift the disc and make it spin around its center at a variable angular speed (between 200 and 500 RPM) so as to keep an approximately constant tangential speed at the laser location. This way, a uniform stream of bits was sent from the data stored optically in a spirally-shaped track containing pits and lands to the music player.

If the electrical motor in charge of spinning the disc failed and its angular speed was too high—and I have witnessed it first hand—the CD would literally shatter and shred up in little shiny pieces as illustrated in the video shown in fig. 1. Of course the “centrifugal force” is to be blamed here and is a source of deformation or possible failure for other parts subject to rotational motion. Is there a way we can simulate and study these effects in Onscale Solve? Sure there is.

1 Equivalent static case

A first naïve approach to tackle the problem would be to set up a transient calculation and write a time-dependent boundary condition at the face of the part being driven by the shaft. We will leave the details of this way of solving the problem for another post. In this one, we focus on the smarter way: to view the problem from a steady non-inertial frame of reference and adding the appropriate pseudo-forces so as to comply with Newton’s mechanical laws.

Let’s assume that our disc lies on the x-y plane, as depicted in fig. 2. The vector representing the angular speed \omega is aligned with the z axis. If we as observers also rotated with the same angular speed \omega around the z axis, we would see the disc as a static part so we would not need to solve a transient problem. Yet there is a catch: Newton’s laws will not hold unless we add pseudo-forces—in particular, the centrifugal forces. The detailed mathematics can be checked in any classical mechanics book (or even in Wikipedia, but the idea is that we have to add an acceleration field

\mathbf{a}(\mathbf{r}) = - \mathbf{\omega} \times \mathbf{\omega} \times \mathbf{r} \qquad{(1)}

Note that the disc has a finite width so the 3D position vector \mathbf{r} of a point on the disc is given by [x,y,z].. Therefore, \mathbf{\omega} \times \mathbf{r} is

\mathbf{\omega} \times \mathbf{r} = \begin{bmatrix}0 \\ 0 \\ \omega\end{bmatrix} \times \begin{bmatrix}x \\ y \\ z\end{bmatrix} = \begin{bmatrix}-\omega \cdot y \\ +\omega \cdot x \\ 0\end{bmatrix}

and \mathbf{\omega} \times \mathbf{\omega} \times \mathbf{r} is

\mathbf{\omega} \times \left( \mathbf{\omega} \times \mathbf{r} \right) = \begin{bmatrix}0 \\ 0 \\ \omega\end{bmatrix} \times \begin{bmatrix}-\omega \cdot y \\ +\omega \cdot x \\ 0\end{bmatrix} = \begin{bmatrix}-\omega^2 \cdot x \\ -\omega^2 \cdot y \\ 0\end{bmatrix}

Hence, minding the negative sign in eq. 1, the acceleration field is

\mathbf{a}(\mathbf{r}) = \omega^2 \begin{bmatrix} x \\ y \\ 0\end{bmatrix} \qquad{(2)}

Not only does the analysis above apply to planar discs but also to any generic part rotating with angular speed \omega around the z axis. Therefore, to model a rotating body as a static problem, we add a distributed acceleration in the outward radial direction proportional to the square of the angular speed. In effect, if r_c is the radial coordinate of a cylindrical coordinate system centered at the original aligned with the z axis, the acceleration is

\mathbf{a}(\mathbf{r}) = \omega^2 \cdot \mathbf{r_c} \qquad{(3)}

2 A solid spinning cylinder

The stresses of a thin uniform solid disc (i.e. a CD without the hole) of radius R rotating with a constant angular speed in the linear elastic regime can be computed analytically for the plane-stress case. Under plane-stress conditions a disc rotating with constant angular speed \omega experiences hoop and radial stress fields respectively equal to [1]:

\sigma_h(r) = \frac{\rho \omega^2}{8} \cdot \Big[ (3+\nu) \cdot R^2 - (1+ 3\nu) \cdot r^2 \Big] \qquad{(4)}

and

\sigma_r(r) = \frac{\rho \omega^2}{8} \cdot (3+\nu) \cdot (R^2-r^2) \qquad{(5)}

where \rho is the density and \nu is the material’s Poisson ratio.

We can now check these results with Onscale. We first open Onshape and draw an angular section of a thin disc, like this one (fig. 3).

Why an “angular section” instead of a full 360º disc? Two reasons:

1. Since the problem is axisymmetric, once we know the result along any radial direction we have the full solution.
2. It is easier to set up the boundary conditions of the problem this way because we have three non-coplanar faces where we can set symmetry conditions. More on this topic on a future post.

After loading the CAD geometry from Onshape into Onscale Solve, let us set up the problem. First, we have to make sure the density is what we need. In this case, let us assume the disc is made of structural steel that in Onscale’s material library has a uniform density \rho=7850~\text{kg}\cdot\text{m}^{-3}.

In this example we set the centrifugal acceleration field as a space-dependent gravity field (in sec. 3 we use a body acceleration field instead). In the Physics step, we “turn on” gravity by opening the Global Settings dialog in Physics as shown in fig. 4. The default values appear. Do not worry about them now, we are going to modify them in the SimAPI code console below.

Next thing we pay attention to is the boundary conditions. Since we are modeling only 30º of the disc, we set symmetry conditions on the sectioned faces of the CAD geometry. We also restrain the disc in the z-direction by adding a symmetry condition on the bottom face to arrive at fig. 5.

Before proceeding, allow me to stress some comments here:

1. The stress fields do not depend on the magnitude of the Young modulus E (the displacement field does).
2. Since the bottom face has a symmetry condition, we are effectively simulating a disc of twice the height in the CAD, i.e. two millimeters instead of one. If we added another symmetry condition of the top face we would be modeling an infinitely long cylinder in the vertical direction, i.e. a plane-strain condition.
3. The choice of 30º was arbitrary. We could have chosen either a smaller or larger angle, with pros and cons. The smaller the angle, the less accurate the solution. The larger the angle, the bigger (and more expensive) the problem we have to solve. For the particular cases of 90º, 180º and 270º the symmetry condition avoids triggering a Lagrange multiplier in Onscale’s solver which might result in a better efficiency, but 30º is accurate and small enough for this case.
4. The disc lies on the horizontal x-y plane and rotates around the vertical  z plane, so we can use the equations we derived in the previous section. In sec. 3 we introduce a more general approach for other directions.

Next we add the centrifugal acceleration we discussed in sec. 1 to set up a static case equivalent to the rotating one. So far we have been using \omega for the angular speed that has units of radians per second. Let us say we want to see what happens to our disc if it rotates at 1,000 revolutions per minute. To obtain the angular speed in radians per seconds we multiply 1,000 revolutions per minutes by 2\pi/60:

\omega~\left[\frac{\text{rad}}{\text{s}}\right] = 1{,}000~\frac{\text{rev}}{\text{min}} \cdot \frac{2\pi~\text{rad}\cdot\text{rev}^{-1}}{60~\text{s} \cdot \text{min}^{-1}} \approx 104.72~\text{rad}\cdot{\text{s}}^{-1}

Since Onscale Solve understands algebraic expressions of space (and time) we can then open the code console and edit the default gravity values to match eq. 2. Find the line

on.settings.Gravity([0,0,-9.81])

and replace it with:

omega = 1000
on.settings.Gravity([f"({omega} * 2*pi/60)^2 * x", f"({omega} * 2*pi/60)^2 * y","0"])

Click “Build” and then “Save.” If you now go back to the “Global Settings” dialog, you should see something like fig. 6.¹

After launching the simulation, we can get the results back. Displacement sensors (fig. 7) are placed at equal intervals along radial line from the center of the disc to compare against the analytical solution given by eqns. 4, 5 as shown in fig. 8.

3 An actual part of rotating machinery

Instead of theoretical solid cylinders, let us now shift our attention to simulating a mechanical part which is subject to centrifugal forces. From the Onscale CAD library, we are going to choose the propeller as shown in fig. 10.

We fully fix the front and back annuli as if the shaft was attached to these faces (fig. 11). Now, in the previous example we used eq. 2 to define the centrifugal acceleration field in the global Cartesian coordinate systems. For the sake of completeness, in this case we are going to use the cylindrical version from eq. 3 by adding an acceleration field to the body (keeping gravity disabled).

c = on.CoordinateSystem(
        vector_1=[1, 0, 0],
        vector_2=[0, 0, 1],
        center=[0, 0, 0],
        system_type="cylindrical"
    ) 
omega = 5000 
centrifugal = on.loads.Acceleration(1, [f'({omega} * 2*pi/60)^2 * {c.r}', "0", "0"], coord_system=c) centrifugal >> geometry.parts[0]

References