Singularity functions are a mathematical tool used to describe systems where discontinuity occurs. These functions are written using a special notation using singularity brackets. The general form is written as $\langle x-a \rangle^n$, where $a$ is the point of discontinuity and $n$ is the function degree. A detailed explanation of the properties of these functions can be found on Wikipedia. The singularity functions method is sometimes referred to as “Macaulay’s method” and the brackets as “Macaulay brackets”.

The main advantage of using singularity functions for solving the continuous beam problem is the ability to describe and solve the system analytically with any kind of load. This can be done for complex loading systems such as having the load described as a high order function. The symbolic solution can be reached without any approximation which gives us a pure analytical solution to any continuous beam.

To demonstrate this, the following continuous beam is used as an example:

Continuous beam problem (lengths in meter) Continuous beam problem (lengths in meter)

The load is described as a function of the length starting at $x=15.75$ m to $x=24.75$ m from left to right. Equation (1) describes the load function; where: $q_0 = 60\ \text{kN}$ and $a = 9\ \text{m}$.

\[q(x) = -\frac{4q_0}{a^2} x^2 + \frac{4q_0}{a} x \tag{1}\]

First we should determine the reaction forces $R_a$, $R_b$ and $R_c$.

Reaction forces Reaction forces

This is done by simplifying the load as follows:

Load simplification Load simplification

\[\begin{array}{l} q(x) = -\frac{4q_0}{a^2} x_0^2 + \frac{4q_0}{a} x_0 \\ F = \int_0^9 q(x)\, dx = 6\, q_0 = 360\ \text{kN} \\ l_0 = \frac{\int_0^9 q(x)\, x\, dx}{F} = 4.5\ \text{m} \end{array} \tag{2}\]

Using simple statics we can determine the system equations (3) and (4):

\[\sum F_y = 0 \rightarrow R_a + R_b + R_c = 6\, q_0 \tag{3}\] \[\sum M_c = 0 \rightarrow R_b + 2R_c = 13.5\, q_0 \tag{4}\]

We start the solution by writing the equation of normal forces as shown in (5). To get the equation, a couple of rules should be known: for single loads we use a singularity function with a degree of $n=-1$, constant line loads get a degree of $n=0$, and so on. Rules for applying singularity functions can be seen in this link. The challenge occurs when having loads with higher order (such as the one in our example). In this case the load should be described as a Taylor expansion. The reason for this is explained in this discussion.

\[\begin{aligned} q(x) &= R_a\langle x-0\rangle^{-1} + R_b\langle x-9\rangle^{-1} + R_c\langle x-18\rangle^{-1} \\ &+ \frac{4q_0}{9}\langle x-15.75\rangle^{1} - \frac{4q_0}{81}\langle x-15.75\rangle^{2} \end{aligned} \tag{5}\]

Now we have finished the most important part of the solution and the rest is simple mathematics! To get the shear force equation, we integrate the force equation (5) as seen in equation (6):

\[\begin{aligned} V(x) &= \int q(x)\,dx \\ V(x) &= R_a\langle x-0\rangle^{0} + R_b\langle x-9\rangle^{0} + R_c\langle x-18\rangle^{0} \\ &+ \frac{4q_0}{18}\langle x-15.75\rangle^{2} - \frac{4q_0}{243}\langle x-15.75\rangle^{3} \end{aligned} \tag{6}\]

Similarly we get the bending moment equation (7):

\[\begin{aligned} M(x) &= \int V(x)\,dx \\ M(x) &= R_a\langle x-0\rangle^{1} + R_b\langle x-9\rangle^{1} + R_c\langle x-18\rangle^{1} \\ &+ \frac{4q_0}{54}\langle x-15.75\rangle^{3} - \frac{4q_0}{972}\langle x-15.75\rangle^{4} \end{aligned} \tag{7}\]

Rotation is calculated in (8). Here we should take an integration constant $C_1$ into account. This is due to the fact that the rotation at $x=0$ does not equal 0:

\[\begin{aligned} EI\theta(x) &= \int M(x)\,dx \\ EI\theta(x) &= \frac{R_a}{2}\langle x-0\rangle^{2} + \frac{R_b}{2}\langle x-9\rangle^{2} + \frac{R_c}{2}\langle x-18\rangle^{2} \\ &+ \frac{4q_0}{216}\langle x-15.75\rangle^{4} - \frac{4q_0}{4860}\langle x-15.75\rangle^{5} + C_1 \end{aligned} \tag{8}\]

Finally we can derive the deflection equation as in (9):

\[\begin{aligned} EIy(x) &= \int EI\theta(x)\,dx \\ EIy(x) &= \frac{R_a}{6}\langle x-0\rangle^{3} + \frac{R_b}{6}\langle x-9\rangle^{3} + \frac{R_c}{6}\langle x-18\rangle^{3} \\ &+ \frac{4q_0}{1080}\langle x-15.75\rangle^{5} - \frac{4q_0}{4860}\langle x-15.75\rangle^{6} + C_1 x \end{aligned} \tag{9}\]

Determining $C_1$: when $x=9$ then $EIy(9) = 0$, therefore by substituting in the deflection equation we determine:

\[0 = \frac{R_a}{6}\langle 9-0\rangle^{3} + 9C_1 \rightarrow C_1 = -13.5\, R_a\]

When $x=18$ then $EIy(18) = 0$, therefore by substituting in the deflection equation we determine equation (10):

\[729\, R_a + 121.5\, R_b = \left(\frac{8019}{40960}\right) q_0 \tag{10}\]

Now by solving equations (3), (4) and (10) we can determine the reaction forces:

\[\begin{aligned} R_a &= \left(\frac{30753}{81920}\right) q_0 \\ R_b &= -\left(\frac{92193}{40960}\right) q_0 \\ R_c &= \left(\frac{645153}{40960}\right) q_0 \\ q_0 &= 60\ \text{kN/m} \rightarrow \\ R_a &= 22.524\ \text{kN} \\ R_b &= -135.048\ \text{kN} \\ R_c &= 472.524\ \text{kN} \end{aligned} \tag{11}\]

Reaction forces result Final reaction forces

Since we already determined the equations of the internal forces and calculated the reactions, we can simply plot the equations by varying the value of $x$.

Interactive Diagrams

The following is the VBA code used to plot the solution. It can be used as a reference on the proper way to implement singularity functions in VBA/Excel:

Public Const Ra = 22.524
Public Const Rb = -135.048
Public Const Rc = 472.524
Public Const q0 = 60
Public Const E = 2.1 * 10 ^ 8
Public Const I = 3945001.5 * 10 ^ -8
Public Const EI = E * I

Function Shear(x)
    y1 = Abs(x >= 0) * Ra * (x - 0) ^ 0
    y2 = Abs(x >= 9) * Rb * (x - 9) ^ 0
    y3 = Abs(x >= 18) * Rc * (x - 18) ^ 0
    y4 = Abs(x >= 15.75) * (-4 * q0 / 18) * (x - 15.75) ^ 2
    y5 = Abs(x >= 15.75) * (4 * q0 / 243) * (x - 15.75) ^ 3

    Shear = Abs(x <= 24.75) * (y1 + y2 + y3 + y4 + y5)
End Function

Function Moment(x)
    y1 = Abs(x >= 0) * Ra * (x - 0) ^ 1
    y2 = Abs(x >= 9) * Rb * (x - 9) ^ 1
    y3 = Abs(x >= 18) * Rc * (x - 18) ^ 1
    y4 = Abs(x >= 15.75) * (-4 * q0 / 54) * (x - 15.75) ^ 3
    y5 = Abs(x >= 15.75) * (4 * q0 / 972) * (x - 15.75) ^ 4

    Moment = -Abs(x <= 24.75) * (y1 + y2 + y3 + y4 + y5)
End Function

Function Rotation(x)
    If x > 24.75 Then
        Rotation = Rotation(24.75)
        Exit Function
    End If

    y1 = Abs(x >= 0) * (1 / 2) * Ra * (x - 0) ^ 2
    y2 = Abs(x >= 9) * (1 / 2) * Rb * (x - 9) ^ 2
    y3 = Abs(x >= 18) * (1 / 2) * Rc * (x - 18) ^ 2
    y4 = Abs(x >= 15.75) * (-4 * q0 / 216) * (x - 15.75) ^ 4
    y5 = Abs(x >= 15.75) * (4 * q0 / 4860) * (x - 15.75) ^ 5
    y6 = -13.5 * Ra

    EItheta = Abs(x <= 24.75) * (y1 + y2 + y3 + y4 + y5 + y6)
    Rotation = EItheta / EI
End Function

Function Deflection(x)
    If x > 24.75 Then
        Deflection = Deflection(24.75) + Rotation(24.75) * (x - 24.75)
        Exit Function
    End If

    y1 = Abs(x >= 0) * (1 / 6) * Ra * (x - 0) ^ 3
    y2 = Abs(x >= 9) * (1 / 6) * Rb * (x - 9) ^ 3
    y3 = Abs(x >= 18) * (1 / 6) * Rc * (x - 18) ^ 3
    y4 = Abs(x >= 15.75) * (-4 * q0 / 1080) * (x - 15.75) ^ 5
    y5 = Abs(x >= 15.75) * (4 * q0 / 29160) * (x - 15.75) ^ 6
    y6 = -13.5 * Ra * x

    EIy = Abs(x <= 24.75) * (y1 + y2 + y3 + y4 + y5 + y6)
    Deflection = EIy / EI
End Function