Worksheet 4

Workspace Setup

Ensure your MATLAB path is set to

<INSTALLDIR>/i2sc_worksheets/Worksheet4

where <INSTALLDIR> is the directory where you extracted the course content (see the Getting Started page)

Shear force and bending moment diagrams

A cantilevered beam under two, positive point loads is shown in Figure.

Cantilevered beam under multiple point loads

The two point loads, which have size \(F_1\) and \(F_2\) and are at a distance from the free-end of the beam of \(x_1\) and \(x_2\) respectively, can also be generalised to \(N\) point loads. The \(n\)^th point load has a strength \(F_n\), and is at \(x_n\) from the free-end of the beam. Given this generalisation, and \(\forall x\in[0,1]\) (where \(x\) is the location anywhere along the beam and the wall is located at \(x=1\)), the shear force, \(V(x)\), and bending moment, \(M(x)\), are given by:

\[\begin{aligned} V(x) &=-\sum_{n=1}^{N} \overline{F_n} \\ M(x) &=-\sum_{n=1}^{N}(x-x_n) \overline{F_n} \end{aligned}\]

where the force term, \(\overline{F_n}\), is given by:

\[\begin{aligned} \overline{F_n} = \left\{ \begin{array}{l l} F_n & \quad x \ge x_n\\ 0 & \quad x < x_n \end{array} \right. \\ \end{aligned}\]

The Task

Complete the following tasks in the supplied file problem_1.m

Write a function with the signature [V,M] = GetLoads(x,forces) to create shear force and bending moment diagrams, where x is a 1xN array of values to evaluate the loads at, forces is a Nx2 array where the first column denotes the \(x\) location of a force, and the second column denotes its magnitude. For example:
```
  forces = [0 2;0.5 1];
```
denotes two forces, a 2N force at \(x=0\) and a 1N force at \(x=0.5\). Add argument validation to your function to ensure all values in the array x are between 0 and 1, throw an error if they are not.
Plot V and M for the loads forces = [0 -2;0.5 -1]
Create a MATLAB function with the signature forces = DiscretiseUniformForce(N,m), which discretises a uniformly distributed load of m Newtons into N point loads. The output forces should be of a format to be used in the function GetLoads(x,forces).
Hint: discretise the load into N point loads, so that the sum of the point forces is \(N\) Newtons.
Plot V and M for a uniformly distributed load of \(-3N\). How does the figure change as you vary the number of discretised point forces?
Imagine the cantilever beam is a wing which is elliptically loaded, e.g. the force per unit length is described as
\[1 = \frac{F^2}{b^2} + (x-1)^2\]
rearrange this equation to solve for F, then find the value of b such that \(\int^1_0{F}dx = m\). Only consider positive values of m.
Hint: Is there a canonical equation for the area of an elipse?
Create a MATLAB function with the signature forces = DiscretiseEllipticalForce(N,m) which discretise a N Newton elliptical load into a series of m point loads and returns them as an array which can be used in the function [V,M] = GetLoads(x,forces). To calculate each point load use the trapezium method
Plot V and M for both a \(3N\) uniform, and \(3N\) elliptical load distribution on the same figures. E.g. one figure for V and another for M.

Solving the heat equation

In this task you will be modelling heat transfer in a one-dimensional rod. This problem deals with a one-dimensional rod of length 1.0 units, where the distance along the bar is denoted \(x\), as shown in Figure.

The temperature at a location along this rod, \(x\), at a specific time \(t\) is given as \(u(x,t)\). The system is governed by the equation

\[\frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} = 0\]

Physically, this equation says that the rate of change of temperature with time (first gradient) is proportional to the curvature of that temperature (second gradient). Hence, the sharper/peakier the temperature in the rod, the quicker it will cool down.

The heat transfer equation is not solvable for all but the simplest of problems. Hence, given a generic problem, the rod must be discretised into a number of nodes and then a numerical solution applied. In this exercise, you will deal with the rod split into evenly spaced nodes. Hence, given a spacing between consecutive nodes, \(\Delta x\), and in a total of \(N=1+(1/\Delta x)\) nodes, the problem is given in Figure.

Schematic of discretised one-dimensional rod

To solve this problem, a numerical method is provided. For this, it is convenient to write the temperature at a given time \(t\) at the \(i\)^th node as \(u_i(t)\). Given the temperature at all nodes at the current time, the temperature after some small increment in time, \(\Delta t\), is given by:

\[u_i(t+\Delta t)=u_i(t)+\frac{1}{2}\Big[ u_{i+1}(t)-2u_{i}(t)+u_{i-1}(t) \Big]\]

where \(\Delta t=\frac{1}{2}(\Delta x)^2\).

Therefore, in pseudo-code form, the general problem structure for simulating this system is:

Determine constants and calculate location of nodes \(x\)
\(t = 0\), specify initial conditions \(u_i(0)\) to all nodes
WHILE suitable stopping criteria not met DO

Increment time level, \(t = t + \Delta t\)

Specify temperature at boundaries: \(u_1(t)\), \(u_N(t)\)

Calculate temperature at next time level for remaining nodes \(u_i(t+\Delta t)\) using numerical scheme

END WHILE

The Task

Write a code to simulate how temperature in the rod varies with time. The code should:
- define the number of elements to discretise the rod into
- define the initial condition of the rod (e.g. the temperature at each node) as an array
- define the boundary conditions for each end of the rod
Using your code, simulate a system in which the rod is initially at a uniform temperature of 2 units. But the ends of the rod are held at a constant lower temperature of 1 unit.
- Ensure you choose a suitable discretisation of time and space
- How long does the rod take to settle to a constant temperate (i.e. when the average temperature stops changing to within 0.01 of a unit)?
Simulate these other systems
- Linear initial temperature, then bringing the cool end to a high fixed temperature
- Constant initial temperature, then having one end hot and the other cool
- Constant initial temperature of zero, one boundary is constant at zero, the other varies according to the function
  \[u_b(t) = sin(2 \pi f t)\]
  where f is the frequency of the oscillation in Hertz. plot the evolution of the temperature with time for different frequencies.