Continuum Problems

The Rayleigh-Ritz method is a variational method. Instead of solving the differential equation, it finds an approximate solution that makes a functional (the total potential energy $\Pi$) stationary. For a stable elastic system, the stationary point is a minimum.

Steps

1. Write the functional. For elastic bodies, $\Pi = U - W$ (strain energy minus work of external loads). For a bar, $\Pi = \dfrac{1}{2}\displaystyle\int_0^L AE\left(\dfrac{du}{dx}\right)^2dx - \int_0^L f u\,dx - Pu(L)$.
2. Choose a trial function with unknown coefficients:

$$u(x) \approx \sum_{i=1}^{n} a_i\,\phi_i(x)$$

The functions $\phi_i$ must be continuous, independent and must satisfy the essential (geometric) boundary conditions. 3. Substitute into $\Pi$ so that it becomes a function of $a_1, a_2, \dots, a_n$. 4. Make $\Pi$ stationary: $\dfrac{\partial \Pi}{\partial a_i} = 0$ for $i = 1, \dots, n$. This gives $n$ algebraic equations $[K]\{a\} = \{F\}$. 5. Solve for $a_i$ and compute displacement, strain and stress.

Remarks

• Natural boundary conditions need not be satisfied by the trial function; they come out of the minimisation.
• More terms give a better answer. The computed stiffness is never less than the true stiffness (the displacement is underestimated), so the method gives a bound.
• If the trial function includes the exact solution, the result is exact.

Limitations

• Choosing one trial function for the whole domain is hard for complex shapes.
• It needs a functional, which exists only for self-adjoint problems.
• FEM removes this difficulty by applying the Ritz idea element by element with simple local trial functions.

Choice of trial function

$w = a\sin(\pi x/L)$ gives $w(0) = w(L) = 0$, which satisfies the essential boundary conditions of a simply supported beam.

Potential energy

$$\Pi = \frac{EI}{2}\int_0^L \left(\frac{d^2w}{dx^2}\right)^2dx - q\int_0^L w\,dx$$

With $w'' = -a\left(\dfrac{\pi}{L}\right)^2\sin\dfrac{\pi x}{L}$, $\displaystyle\int_0^L\sin^2\dfrac{\pi x}{L}\,dx = \dfrac{L}{2}$ and $\displaystyle\int_0^L\sin\dfrac{\pi x}{L}\,dx = \dfrac{2L}{\pi}$:

$$\Pi = \frac{EI}{2}\,a^2\left(\frac{\pi}{L}\right)^4\frac{L}{2} - qa\,\frac{2L}{\pi}$$

Stationary condition

$$\frac{d\Pi}{da} = EI\,a\left(\frac{\pi}{L}\right)^4\frac{L}{2} - \frac{2qL}{\pi} = 0 \Rightarrow a = \frac{4qL^4}{\pi^5EI}$$

Substituting values:

$$a = \frac{4(10\,000)(4)^4}{\pi^5(1.6\times10^6)} = \frac{1.024\times10^7}{306.02\times1.6\times10^6} = 0.02091\ \text{m}$$

Mid-span values

• Deflection: $w(L/2) = a = 20.91$ mm.
• Bending moment: $M = EI\,w''$ in magnitude at $x = L/2$: $M = EI\,a\left(\dfrac{\pi}{L}\right)^2 = 1.6\times10^6(0.02091)(0.6169) = 20\,641$ N m.

Comparison

Deflection (the primary variable) is accurate; the moment (a second derivative) is less accurate because differentiation reduces accuracy.

Answer: $a = 20.9$ mm; $w_{max} = 20.9$ mm; $M_{max} = 20.6$ kN m.

Residual

With $u = a\,x(1-x)$: $u'' = -2a$. The trial function satisfies $u(0) = u(1) = 0$.

$$R(x) = u'' + u + x = -2a + a\,x(1-x) + x$$

(a) Galerkin method

The weight function equals the trial function shape, $W = x(1-x)$, and $\int_0^1 W R\,dx = 0$.

Needed integrals: $\int_0^1 x(1-x)dx = \dfrac{1}{6}$, $\int_0^1 x^2(1-x)^2dx = \dfrac{1}{30}$, $\int_0^1 x^2(1-x)dx = \dfrac{1}{12}$.

$$\begin{aligned} \int_0^1 x(1-x)\,R\,dx &= -2a\left(\frac{1}{6}\right) + a\left(\frac{1}{30}\right) + \frac{1}{12} \\ &= -\frac{3a}{10} + \frac{1}{12} = 0 \end{aligned}$$ $$a = \frac{10}{36} = \frac{5}{18} = 0.2778$$

$u(0.5) = 0.2778(0.5)(0.5) = 0.06944$.

(b) Point collocation at $x = 0.5$

$$R(0.5) = -2a + a(0.25) + 0.5 = -1.75a + 0.5 = 0 \Rightarrow a = \frac{2}{7} = 0.2857$$

$u(0.5) = 0.2857(0.25) = 0.07143$.

Comparison

Exact: $u(0.5) = \dfrac{\sin 0.5}{\sin 1} - 0.5 = \dfrac{0.47943}{0.84147} - 0.5 = 0.06975$.

Galerkin averages the residual over the whole domain, so it is more accurate than collocation, which forces zero residual only at one point.

Answer: Galerkin $u = 0.2778\,x(1-x)$, $u(0.5) = 0.0694$; collocation $u = 0.2857\,x(1-x)$, $u(0.5) = 0.0714$; exact $0.0698$.

Method of weighted residuals (MWR)

For a differential equation $\mathcal{D}(u) + f = 0$ on a domain $\Omega$, an approximate solution $\tilde{u} = \sum a_i\phi_i(x)$ does not satisfy the equation exactly, leaving a residual

$$R(x, a_i) = \mathcal{D}(\tilde{u}) + f \neq 0$$

The MWR makes the residual small by forcing a weighted integral to be zero:

$$\int_\Omega W_i(x)\,R(x, a_i)\,d\Omega = 0,\qquad i = 1, \dots, n$$

This gives $n$ equations for the $n$ unknowns $a_i$. The choice of the weight function $W_i$ defines the method.

Comparison

Remarks

• Collocation is simplest but depends on the choice of points.
• Least squares always gives a symmetric matrix but needs higher-order derivatives, so the algebra is heavy.
• Galerkin gives a symmetric matrix for self-adjoint problems and equals the Ritz method when a functional exists. It is the basis of most finite element formulations.
• All the methods become exact if the residual is zero everywhere.

(a) Strong form

Find $u(x)$ such that

$$\frac{d}{dx}\left(AE\frac{du}{dx}\right) + f = 0,\quad 0 < x < L$$

with $u(0) = 0$ (essential BC) and $AE\dfrac{du}{dx}\Big|_{x=L} = P$ (natural BC).

It needs $u$ to be twice differentiable, and it must hold at every point.

(b) Weak form (Galerkin)

Multiply the equation by an arbitrary weight (test) function $w(x)$ with $w(0) = 0$, and integrate over the length:

$$\int_0^L w\left[\frac{d}{dx}\left(AE\frac{du}{dx}\right) + f\right]dx = 0$$

Integrate the first term by parts:

$$\int_0^L w\,\frac{d}{dx}\left(AE\frac{du}{dx}\right)dx = \left[w\,AE\frac{du}{dx}\right]_0^L - \int_0^L \frac{dw}{dx}AE\frac{du}{dx}\,dx$$

At $x = 0$, $w = 0$. At $x = L$, $AE\,u' = P$. So the boundary term is $w(L)P$. The weak form is:

$$\int_0^L \frac{dw}{dx}\,AE\,\frac{du}{dx}\,dx = \int_0^L w f\,dx + w(L)\,P$$

for all admissible $w$.

Note the natural BC entered through the boundary term, and the essential BC is imposed on $u$ and $w$.

(c) Why weak form is preferred

• Only first derivatives appear, so the trial and test functions need to be only $C^0$ continuous; simple linear shape functions can be used.
• Natural boundary conditions are included automatically.
• It leads to a symmetric stiffness matrix, $K_{ij} = \int AE\,\phi_i'\phi_j'\,dx$.
• It is the physical statement of the principle of virtual work, and it works when $A$, $E$ or the load are discontinuous.

Data (units: N, mm)

$AE = 300 \times 70\times10^3 = 2.1\times10^7$ N, $L = 2000$ mm, $q = 6$ N/mm, $P = 12\,000$ N.

The trial function satisfies $u(0) = 0$.

Potential energy

$$\Pi = \frac{AE}{2}\int_0^L\left(\frac{du}{dx}\right)^2dx - q\int_0^L u\,dx - P\,u(L)$$

With $u' = a_1 + 2a_2x$:

$$\begin{aligned} \int_0^L u'^2dx &= a_1^2L + 2a_1a_2L^2 + \tfrac{4}{3}a_2^2L^3 \\ \int_0^L u\,dx &= \tfrac{1}{2}a_1L^2 + \tfrac{1}{3}a_2L^3 \\ u(L) &= a_1L + a_2L^2 \end{aligned}$$

Stationary conditions

$$\begin{aligned} \frac{\partial\Pi}{\partial a_1} &= AE\,(a_1L + a_2L^2) - \tfrac{1}{2}qL^2 - PL = 0 \\ \frac{\partial\Pi}{\partial a_2} &= AE\,(a_1L^2 + \tfrac{4}{3}a_2L^3) - \tfrac{1}{3}qL^3 - PL^2 = 0 \end{aligned}$$

Divide the first by $L$ and the second by $L^2$:

$$\begin{aligned} AE\,(a_1 + a_2L) &= \tfrac{1}{2}qL + P = 6000 + 12\,000 = 18\,000 \\ AE\,(a_1 + \tfrac{4}{3}a_2L) &= \tfrac{1}{3}qL + P = 4000 + 12\,000 = 16\,000 \end{aligned}$$

Subtract: $AE\,\tfrac{1}{3}a_2L = -2000$, so $a_2 = \dfrac{-6000}{AE\,L} = \dfrac{-6000}{4.2\times10^{10}} = -1.4286\times10^{-7}\ \text{mm}^{-1}$.

Then $a_1 = \dfrac{18\,000}{AE} - a_2L = 8.5714\times10^{-4} + 2.8571\times10^{-4} = 1.1429\times10^{-3}$.

End displacement

$$u(L) = a_1L + a_2L^2 = 2.2857 - 0.5714 = 1.714\ \text{mm}$$

Exactness

The exact solution of $AE\,u'' + q = 0$ with $u(0) = 0$ and $AE\,u'(L) = P$ is

$$u = \frac{(P + qL)}{AE}x - \frac{q}{2AE}x^2$$

which is a quadratic. Our coefficients match it ($\frac{24\,000}{2.1\times10^7} = 1.1429\times10^{-3}$ and $\frac{6}{4.2\times10^7} = 1.4286\times10^{-7}$). So the result is exact because the trial function contains the exact solution.

Answer: $a_1 = 1.143\times10^{-3}$, $a_2 = -1.429\times10^{-7}\ \text{mm}^{-1}$; end displacement $= 1.714$ mm; yes, exact.