Numerical Techniques in Engineering

Numerical Techniques in Engineering
This is a course work on Numerical Techniques in Engineering (ENGT5140),
which is the part of the module assessment. This course work determines 50%
of your grade for the module.
Please note that there is a deadline for handing in this work. All coursework
will be marked against the following criteria:
• Correctness of results (50%), including mathematical derivations and ac-
curacy of simulation
• Elegance and style of computer programming work (30%)
• Elegance and clarity of presentation (20%). Your work should be presented
as if writting for an audience of fellow students, and explain everything
they would need to reproduce your work.
• Any course work must be conducted individually, cheating of any kind
will not be tolerated. You may be asked for a presentation or an oral
examination based on your course work, if necessary.
The submission deadline is Friday 11th January 2013, before 2:00pm and
handed in to the Faculty of Technology Support Centre on the ground floor
of the Gateway House with a coursework cover sheet. Late submission will be
penalised.
In your submitted coursework, you should include:
• Details of your coursework report.
• Programming codes in Matlab or other computer languages.
• Any necessary plots, tables and diagrams.
1
Coursework Questions
(You need to choose THREE QUESTIONS ONLY as your coursework)
1. A cylindrical pipe has a hot fluid flowing through it. Because the pressure is very high,
the walls of the pipe are thick. For such a situation, the temperature distribution in
the metal wall related to radial distance can be modelled by a second order differential
equation
r
d2
dr2 +
d
dr
= 0
where  is temperature, and r is the radial distance from the centreline.
(a) Solve for the temperature (r) within a pipe whose radius is 1 cm and whose
outer radius is 2 cm if the fluid is at 540oC and the temperature of the outer
circumference is 20oC.
(b) The pipe is insulated to reduce the heat loss. The insulation used has proper-
ties such that the gradient
d
dr
at the outer circumference is proportional to the
difference in temperature from the outer wall to the surroundings, which can be
modelled by
d
dr

r=2
= 0.083((2) – 20)
Solve the equation with this boundary condition.
(c) In this question, you should use at least two numerical methods to solve the
equation, and compare the solutions by each method. Plot solutions and make
your judgement for each method.
2. A cable of constant length passing frictionless rings is used as shown in Fig.1. Ini-
tially, both masses m1 and m2 are held at rest at positions shown. When they are
released, their instantaneous positions can be denoted as y(t) and z(t), respectively.
The instantaneous angle and the length of the inclined portion of the cable can be
expressed in term of y as:
 = tan-1

y + h
b

L = [(y + h)2 + b2]1/2
(1)
If we denote the cable tension as T, then the Newtons second law applied to the two
masses leads to
m1g – 2T sin  = m1
d2y
dt2
T – m2g = m2
d2z
dt2
(2)
By eliminating T, we obtain:
d2y
dt2 +
2m2
m1
sin 
d2z
dt2 = g

1 –
2m2
m1
sin 

(3)
2
Figure 1: Diagram for two connected masses m1 and m2.
The displacements y and z are restricted by the condition that the cable length must
remain unchanged. That is
z = 2
p
[(y + h)2 + b2]
p
[h2 + b2]
This relationship can be differentiated with respect to t to obtain another equation
relating
d2y
dt2 and
d2z
dt2 which is
d2z
dt2 =
-2
L2 (y + h)
dL
dt
dy
dt
+
2
L

dy
dt
2
+ (y + h)
d2y
dt2
!
(4)
Equation (4) can be substituted into Equation (3) to obtain
d2y
dt2 =
sin 
L(2y + h + L)

gL2 –
2m2
m1
”
gL2 – 2(y + h)
dL
dt
dy
dt
+ 2L

dy
dt
2
#!
(5)
(a) By letting x1 = y, x2 =
dy
dt
, x3 = z and x4 =
dz
dt
, show that the final form of
the system from Equations (4) and (5) can be written as a set of four first order
differential equations
dx1
dt
= x2
dx2
dt
=
GL2(1 – 2rm sin ) + 4rmv2 sin 
L(1 + 4r sin2 )
dx3
dt
= x4
dx4
dt
=
-2v2
L
+ 2
dx2
dt
sin 
where
rm =
m2
m1
, v2 = x2

dL
dt
sin  – x2

(b) Investigate a numerical case where b = 3.6m, h = 0.15m,
m2
m1
= 0.8, and initial
conditions are
y = z =
dy
dt
=
dz
dt
= 0
3
Plot the numerical solution for 0  t  12 seconds, and find their maximum
speeds, |dy/dt| and |dz/dt|.
3. Fig.2 is the circuit diagram for a small network of resistors. There are five nodes,
Figure 2: Diagram for a resistor network.
eight resistors, and one constant voltage source. We want to compute the voltage
drops between the nodes and the currents around each of the loops.
Several different linear systems of equations can be formed to describe this circuit.
Let vk, k = 1, · · · , 4, denote the voltage difference between each of the first four nodes
and node number 5 and let ik, k = 1, · · · , 4, denote the clockwise current around each
of the loops in the diagram. Ohm’s law says that the voltage drop across a resistor
is the resistance times the current. For example, the branch between nodes 1 and 2
gives
v1 – v2 = r12(i2 – i1)
Using the conductance, which is the reciprocal of the resistance, gkj =
l
rkj
, Ohm’s
law becomes
i2 – i1 = g12(v1 – v2)
The voltage source is included in the equation
v3 – vs = r35i4
Kirchhoff’s voltage law says that the sum of the voltage differences around each loop
is zero. For example, around loop 1,
(v1 – v4) + (v4 – v5) + (v5 – v2) + (v2 – v1) = 0
Combining the voltage law with Ohm’s law leads to the loop equations for the currents:
Ri = b
Here i is the current vector, and b is the source voltage vector,
i =
.
…
i1
i2
i3
i4
.
…
, b =
.
…
0
0
0
vs
.
…
4
and R is the resistance matrix,
R =
.
…
r25 + r12 + r14 + r45 -r12 -r14 -r45
-r12 r23 + r12 + r13 -r13 0
-r14 -r13 r14 + r13 + r34 -r34
-r45 0 -r34 r35 + r45 + r34
.
…
Kirchhoff’s current law says that the sum of the currents at each node is zero. For
example, at node 1,
(i1 – i2) + (i2 – i3) + (i3 – i1) = 0
Combining the current law with the conductance version of Ohm’s law leads to the
nodal equations for the voltages:
Gv = c
Here v is the voltage vector and c is the source current vector:
v =
.
…
v1
v2
v3
v4
.
…
, c =
.
…
0
0
g35vs
0
.
…
and G is the conductance matrix,
G =
.
…
g12 + g13 + g14 -g12 -g13 -g14
-g12 g12 + g23 + g25 -g23 0
-g13 -g23 g13 + g23 + g34 + g35 -g34
-g14 0 -g34 g14 + g34 + g45
.
…
You can solve the linear system obtained from the loop equations to compute the
currents and then use Ohm’s law to recover the voltages. Or you can solve the linear
system obtained from the node equations to compute the voltages and then use Ohm’s
law to recover the currents. Your assignment is to verify that these two approaches
produce the same results for this circuit. You can choose your own numerical values
for the resistances and the voltage source.
4. A cable hung at its two ends as shown in Figure 3 by its own weight will have a
catenary shape described by the equation
y =
Tx
w

cosh
wx
Tx
– 1

(6)
where w is the weight per unit length and Tx is the horizontal, x-component of the
tension of the cable. Equation (6) is for the case when both w and Tx are constant
throughout the cable. In fact, Equation (6) is the solution of the differential equation
d2y
dx2 =
w
Tx
ds
dx
=
w
Tx
”
1 +

dy
dx
2
#1/2
(7)
where s is a variable along the length of the cable.
(a) Applying the Runge-Kutta method to solve the equation (6), where w = 0.12
kN/m, xA = yA = 0, xB = 200 m, and yB = 50 m. Let the initial conditions be
y = 0 at x = xA = 0, iterate Tx value until yB is within 99.9% of 50 m.
5
Figure 3: A diagram for a cable.
(b) How could the problem be solved if xA = 100 m and yA = 25 m by application
of the Runge-Kutta method (w remains equal to 0.12 kN/m)?
5. Consider the problem of an orbit of a satellite, whose position and velocity are obtained
as the solution of the following equations
x’
1(t) = x3(t)
x’
2(t) = x4(t)
x’
3(t) =
-GMEx1(t)
(x2
1(t) + x2
2(t))3/2
x’
4(t) =
-GMEx2(t)
(x2
1(t) + x2
2(t))3/2
(8)
where G = 6.672×10-11 N m2/kg2 is the gravitational constant, andME = 5.97×1024
kg is the mass of the earth. Note that (x1, x2) and (x3, x4) denote the position and
velocity, respectively, of the satellite on the plane having the earth at its origin. This
state equation is defined in the M-file orbit.m below (only for reference)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dx = df_sat(t,x)
global G Me Re
dx = zeros(size(x));
r = sqrt(sum(x(1:2).^2));
if r <= Re, return; end % when colliding against the earth surface
GMr3 = G*Me/r^3;
dx(1) = x(3); dx(2) = x(4); dx(3) = ??????; dx(4) = ??????????????;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%orbit.m to solve a nonlinear d.e. on the orbit of a satellite
clear, clf
global G Me Re
G = 6.67e-11; Me = 5.97e24; Re = 64e5;
t0 = 0; T = 24*60*60; tf = T; N = 2000;
R = 4.223e7;
v20s = [3071 3500 2000];
for iter = 1:length(v20s)
6
x10 = R; x20 = 0; v10 = 0; v20 = v20s(iter);
x0 = [x10 x20 v10 v20]; tol = 1e-6;
[tR,xR] = ode_RK4(f,[t0 tf],x0,N);
[t45,x45] = ode45(????????????);
[t23s,x23s] = ode23s(???????????????????);
plot(xR(:,1),xR(:,2),b, x45(:,1),x45(:,2),k., ????????????)
[t45,x45] = ode45(??????????????????????????));
[t23s,x23s] = ode23s(?????????????????????????????????);
plot(xR(:,1),xR(:,2),b, x45(:,1),x45(:,2),k., ????????????)
end
(a) Supplement the following program orbit.m which uses the three routines ode_RK4(),
ode45(), and ode23s() to find the paths of the satellite with the following initial
positions/velocities for one day.
i. (x10, x20) = (4.223 × 107, 0) [m] and (x30, x40) = (v10, v20) = (0, 3071) [m/s]
ii. (x10, x20) = (4.223 × 107, 0) [m] and (x30, x40) = (v10, v20) = (0, 3500) [m/s]
iii. (x10, x20) = (4.223 × 107, 0) [m] and (x30, x40) = (v10, v20) = (0, 2000) [m/s]
Rewrite the program and check if the plotting results are as depicted in Figure
4.
(b) In Fig. 4, we see that the ode23s() solution path differs from the others for
case (ii) and the ode45() and ode23s() paths differ from the ode_RK4() path
for case (iii). But, we do not know which one is more accurate. In order to find
which one is the closest to the true solution, apply the two routines ode45() and
ode23s() with smaller relative error tolerance of tol = 1e – 6 to find the paths
for the three cases. Which one do you think is the closest to the true solution
among the paths obtained by these methods? Give a detail of your discussion.
(c) The purpose of this problem is not to compare the several MATLAB routines,
but to warn the users of the danger of abusing them. With smaller number of
steps (N) (i.e., larger step size), the routine ode_RK4() will also deviate much
from the true solution. The MATLAB built-in routines have too many good
features to be mentioned here. Note that setting the parameters such as the
relative error tolerance (RelTol) is sometimes very important for obtaining a
reasonably accurate solution.
6. A typical steady-state heat flow problem is the following: A thin steel plate is a 10×20
(cm)2 rectangle. If one side of the 10 cm edge is held at 1000C and the other three
edges are held at 00C, what are the steady-state temperature at interior points? We
can state the problem mathematically in this way if we assume that heat flows only
in the x and y directions:
Find u(x, y) (temperature) such that
@2u
@x2 +
@2u
@y2 = 0 (9)
with boundary conditions
u(x, 0) = 0
u(x, 10) = 0
u(0, y) = 0
u(20, y) = 100
7
Figure 4: The paths of a satellite with the same initial position and different initial velocities.
We replace the differential equation by a difference equation
1
h2 [ui+1,j + ui-1,j + ui,j+1 + ui,j-1 – 4ui,j ] = 0 (10)
which relates the temperature at the point (xi, yj) to the temperature at four neigh-
boring points, each the distance h away from (xi, yj). An approximation of Equation
(9) results when we select a set of such points (these are often called as nodes) and
find the solution to the set of difference equations that result.
(a) If we choose h = 5 cm , find the temperature at interior points.
(b) Write a program to calculate the temperature distribution on interior points with
h = 2.5, h = 0.25, h = 0.025 and h = 0.0025 cm. Discuss your solutions and
examine the effect of grid size h.
7. Modified the difference equation (10) so that it permits to solve the equation
@2u
@x2 +
@2u
@y2 = xy(x – 2)(y – 2)
on the region
0  x  2, 0  y  2
with boundary condition u = 0 on all boundaries except for y = 0, where u = 1.0.
Write and run the program with different grid sizes h and discuss your numerical
results.
8. If a antilever beam of length L, which bends due to a uniform load of w lb/ft, is also
subject to an axial force P at its free end (see Figure 5). The equation of its elastic
8
curve is
EI
d2y
dx2 = Py –
wx2
2
(11)
For this equation, the origin O has been taken at the free end. I is the moment of
inertia; here I =
bh3
12
. At the point x = L,
dy
dx
= 0 and at x = 0, y = 0.
PO X
YLFigure 5: Diagram for a cantilever beam.
(a) Solve the boundary value problem by any numerical method for a wooden beam,
2 in × 4 in × 10 ft, where E = 12 × 105 lb/in2. Find y versus x when the beam
has the 4 in dimension vertical with w = 25 lb/ft and a tension force of P = 500
lb.
(b) And also solve for the deflections if the beam is turned so that the 4 in dimension
is horizontal.
(c) If the factor of the radius of curvature has been taken into account, Equation
(11) becomes
EI
[1 + (y’)2]3/2
d2y
dx2 = Py –
wx2
2
(12)
If the defection of the beam is small, the difference is negligible, but in some cases
this is not true. Furthermore, if there is considerable bending of the beam, the
horizontal distance from the origin to the wall is less than L, the original length
of the beam. Solve the equation (12) with the same conditions, and determine
by how much the deflections differ from those previous calculations.
9. Condon and Odishaw (1967) discuss Duffing equation for the flux  in a transformer.
This nonlinear differential equation is
d2
dt2 + !2
0 + b3 =
!
N
E cos(!t), (0) = '(0) = 0
In this equation , E cos !t is the sinusoidal source voltage and N is the number of
turns in the primary winding, where !0 and b are parameters of the transformer design.
Make a plot of  versus t (and compare to the source voltage) if E = 165, ! = 120,
N = 600, !2
0 = 83 and b = 0.14. For approximate calculations, the nonlinear term
b3 is sometimes neglected. Evaluate your result to determine whether this makes a
significant error in the results.
(Condon E. and Odishaw H., (1967) Handbook of Physics, New York, McGraw-Hill)
9
10. The Lorenz equations
The Lorenz equations are a simplification of a miniature atmosphere model that he
designed to study Rayleigh-Benard convection, the movement of heat in a fluid, such
as air from a lower warm medium (such as the ground) to a higher cool medium (like
the upper atmosphere). In this model of a two-dimensional atmosphere, a circulation
of air develops that can be described by the following system of three equations
x’ = -sx + sy
y’ = -xz + rx – y
z’ = xy – bz
(13)
The variable x denotes the clockwise circulation velocity, y measures the temperature
difference between the ascending and descending columns of air, and z measures the
deviation from a strictly linear temperature profile in the vertical direction. The
Prandtl number s, the Rayleigh number r and b are parameters of the system. The
most common setting for the parameters is s = 10, r = 28 and b = 8/3. These setting
would be used in your calculations.
The Lorenz equations are an important example because the trajectories show great
complexity, despite the fact that the equations are deterministic and fairly simple
(almost linear). The explanation for the complexity is similar to the double pendulum
or three body problem: sensitive dependence on initial conditions. Consider the initial
condition (x, y, z) = (5, 5, 5, ) and another initial condition at a distance
= 10-5
from the first. Compute both trajectories by fourth-order Runge-Kutta with step
sizes h = 0.001 and h = 0.0001, respectively. Explore the sensitive dependence of this
so-called chaotic attractor.
10

ORDER THIS ESSAY HERE NOW AND GET A DISCOUNT !!!

You can place an order similar to this with us. You are assured of an authentic custom paper delivered within the given deadline besides our 24/7 customer support all through.