AMME2000/BMET2960/BMET9960 - Assignment 2, 2023

•   Your report must be your own work. Please see theFAQ on Academic Integrityfor informaon on what is appropriate and what is not appropriate for Assignments.

•   Present your assignment as a concise report in PDF format;

•   The report must be typed, including all mathemacal working; handwriten assignments will receive a mark of 0 and marks will be deducted for screenshots of equaons and/or figures;

•   Your MATLAB code must be included as an appendix to the assignment report. 10% of the assignment marks are allocated for the readability and quality of your MATLAB code;

•   The report should not exceed 10 pages (this limit does NOT include the code in the appendix);

•   Structure the report using the appropriate Secon and Queson numbers (e.g. 1.1, 2.2). Ensure all relevant working and final answers are in that secon;

•   All figures and tables in your report must be numbered and must be referred to and discussed;

•   Submit your report to Turnin by the due date. Late submissions incur a penalty of 5%/day. We have prepared anexemplary assignmentto guide you in wring your own.

QUESTION 1: WAVE EQUATION [40 marks]

Relevant parts of the course for this queson:

•   Week 6 Lecture 1

•   Week 6 Lecture 2

•   Week 7 Pre-work video: ‘The starting problem’

•   Week 7 tutorialandWeek 8 tutorial

A torsional wave can propagate along a stiff rod if the rod experiences a torque. The resulting wave can be modelled using the wave equation:

a2e     a2e

2

where e = e(x, t) is the torsion angle of the wave in radians at location x (along the length of the rod) and time t, and c is given by:

c = l                                          (1.2)

where G is the shear modulus of the rod in units of Pa, and p is the density of the rod in units of kg.m-3 . If the boundary conditions are e(0, t) = e(L, t) = 0, the general solution to Eqn (1. 1) is given by:

e(x, t) = x[Ancos(入nt) + Bnsin(入nt)] sin (nL(冗)x)

n=0                    (1.3)

where 入n = nL(几)c are the eigenvalues, L is the length of the rod, and An and Bn are coefficients which depend on the initial conditions. Consider a thin metallic filament of length L being used as a sensor in an engineering application (Figure 1.1):

The filament is prone to torsional waves, which in turn impact the sensor function. To fully understand the impact, you want to characterise the torsion e = e(x, t) in the sensor filament. You may assume the boundary conditions  are e(0, t) = e(L, t) = 0 so that the solution is given by Eqn (1.3).

1.1.   If the initial torsion angle along the length of the rod is 0 everywhere and the initial torsional velocity along

the rod is given by g(x), shown in Figure 1.2, write down the displacement initial condition and the velocity initial condition mathematically. [4 marks]

1.2.   Starting from Eqn (1.3), show why for this particular problem the solution reduces to:

e(x, t) =x(∞) Bnsin(入nt) sin (nL(几)x)

=1

[2 marks]

1.3.    Given the initial velocity profile (Figure 1.2), solve for the Bn coefficients in Eqn (1.4). [4 marks]

Note: You may use the fact that the Fourier sine series for a triangular function of the form:

is given by:

f(x) = n几K几2) sin (nL(几)x)

Therefore, you do NOT need to perform any integrations yourself for this question.

(1.4)

(1.5)

1.4.   Implement your analytical solution for the torsional wave in MATLAB. For this implementation use the

following parameter values:

•    L = 0.1 m

•    G = 30 GPa

•    p = 9 × 103 kg. m − 3

•    g0  = 50 rad. s − 1

•    nx  = 101 (number of points in spatial domain)

In your report, provide a plot of your solution at 4 different time points, all overlaid on the same plot. Make  sure you indicate the time corresponding to each curve and provide a brief discussion of your results. Include your MATLAB code in the appendix of your report. [10 marks]

1.5.   Write down the form of the CTCS scheme to solve this problem numerically. (Note: You do not need to

derive the CTCS scheme.) Briefly justify why the CTCS scheme is a suitable choice to solve this problem. (Hint: Consider the properties of the scheme and the expected physics of the problem. You do not need to  include any derivations for this part.) [2 marks]

1.6.   Briefly describe what is meant by the ‘starting problem’ in the CTCS scheme and explain how you will

handle it for this problem. You should outline the specific steps and how you plan to implement those steps in MATLAB. [4 marks]

1.7.   Now implement the CTCS numerical scheme in MATLAB to solve the problem numerically. For the

numerical implementation you should use the following parameters:

•    ∆t = the maximum stable time step

•    nx  = 101 (number of points in spatial domain)

In your report you should include the following:

•    Explain how you computed the maximum stable time step. (Note that you are NOT required to perform a von Neumann stability analysis.) [2 marks]

•    Provide a plot showing the analytical and numerical solutions overlaid at the same 4 time points used in 1.4. Briefly discuss your results. [6 marks]

Note: Make sure your implementation compares the analytical and numerical solutions at the same (not different) time at each step.

1.8.   Plot the maximum absolute percentage error between the analytical and numerical solution as a function of

time over one complete cycle of the wave up and down the rod. On the same plot overlay the mean absolute  percentage error as a function of time. Use 51 evenly spaced time points for your horizontal axis. (Note: In   your computation of the mean and maximum absolute percentage error, make sure you avoid any time points or spatial locations that will introduce a divide-by-0 error.) Briefly explain how you generated the plot and    discuss the results. [6 marks]

QUESTION 2: LAPLACE EQUATION [40 marks]

Relevant parts of the course for this queson:

•   Week 7 Lecture 1

•   Week 7 Lecture 2

•   Week 9 tutorial

An engineering control experiment involves exposing a small rectangular wafer to several lasers. The laser light results in a ‘photon pressure’ distribution on the wafer which we can model using the Laplace equation:

 +  = 0                                     (2.1)

where S = S(x,y) is the 2-dimensional pressure distribution on the wafer in units of Pa. Figure 2.1 shows the        wafer dimensions and the pressure distribution at the boundary corresponding to x = a. The other boundaries have no pressure applied.

2.1    For the scenario shown in Figure 2. 1, write down the 4 boundary conditions for the pressure S

mathematically. [4 marks]

2.2    Use separation of variables to show that the steady-state pressure distribution for the scenario in Figure 2.1 is

given by:

S(x, y) = x An(∗) sinh (nb(几)x) sin (nb(几)y)                                                (2.2)

=1

where An(∗) are coefficients yet to be determined. [4 marks]

2.3    Apply the non-zero boundary condition to solve for the An(∗) coefficients and hence write down the complete

solution for S(x, y). [4 marks]

Note: For this question you may use Eqn (1.5) to help you solve for An(∗) – i.e. you should NOT perform any integrations for this part.

2.4    Suppose the final arrangement of the lasers results in the boundary conditions shown in Figure 2.2:

Briefly explain how you would obtain the steady-state solution for S(x, y) for this new scenario, and write down this solution. Note: You do NOT need to perform any lengthy derivations, you should simply explain what mathematical principles allow you to adapt your solution from 2.3 to this new solution. [4 marks]

2.5     Implement your analytical solution from 2.4 in MATLAB. For your implementation use the following

parameters:

•    a = b = 0.1 m

•    nx = ny = 101 (number of points in each dimension of the spatial domain)

•    No. of Fourier terms = 100

•    S1 = S2 = 2 × 10−3 Pa

In your report, provide a clearly labelled contour plot of the steady-state solution. Note that you can use the pcolor function in MATLAB to generate this plot. Be sure to provide a brief discussion of your plot. Include your MATLAB code in the appendix of your report. [8 marks]

2.6    Now implement the Gauss-Jacobi numerical scheme from Week 7 Lecture 2 to solve the same problem

numerically. In your report:

•    Briefly outline the steps you used to implement the Gauss-Jacobi method, including your stopping criteria. [4 marks]

•    Present 2 plots side-by-side – the left plot should be your numerical solution to the steady-state         pressure distribution; the right plot should be the signed difference image between the numerical and analytical solutions. Discuss your results. [8 marks]

2.7     The scientists want to know how much they can shorten the wafer before the pressure at the dead-centre

becomes too great. Using your analytical solution, determine at what wafer width, ath , the pressure at the centre first exceeds 75% of S1 . You should determine ath using a tolerance of 1 µm. Repeat your estimate using your numerical solution. In your report:

•    Briefly outline the method you used to compute ath . [2 marks]

•    Present your results using the analytical and numerical solutions and briefly discuss any differences.

[2 marks]

QUESTION 3: Fourier Integral [10 marks]

Relevant parts of the course for this queson:

•   Week 8 Lecture 2

•   Week 9 Lecture 1

•   Week 10 tutorial

A breach in a factory pipeline has caused a large amount of the insecticide permethrin to be released into a long narrow waterway (Figure 3. 1). The concentration C(x, t) of permethrin at location x and time t along the waterway can be modelled using Fick’s second law, as the heat equation:

where D is the mass diffusivity. Because the waterway is very long, there are no boundary conditions.

3.1    State the boundary conditions mathematically and briefly explain why you cannot solve this problem analytically using Fourier series. [2 marks]

3.2    C(x, t) can be obtained using the Fourier Integral form of the heat equation solution, as shown in Eqn (3.2):

C(x, t) = 1 j∞ [j∞ f(v) cos(px − pv)e−Dp2 t dv| dp                                 (3.2)

where p is a continuous variable and f(x) = f(v) is the initial condition. Using the known integral on the Data Sheet:

j0∞- −sz  cos(2bs) . ds =2(√几) - −bz

show how Eqn (3.2) can be simplified to:

C(x, t) = 2√D(1)几t j−f(v)- − ( (z  )dv

[2 marks]

(3.3)

(3.4)