xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

MATLAB 代写-BMET2960/BMET9960-Assignment 2

时间：2021-05-05

AMME2000/BMET2960/BMET9960 - Assignment 2

Due: 11:59 pm Friday 4th June 2021

Introduction

A team is using drones to capture high resolution images of the terrain after recent volcanic activity

in Iceland. However, the images are distorted due to variations in air pressure and density caused

by the heating from the lava. The team has provided you with some data you can use to obtain a

model of the pressure and density fields in the environment - which will then allow them to apply

corrections to their drone images.

In this assignment you will use the analytical and numerical techniques developed in this course

to obtain a steady-state temperature field in a 2D rectangular slice of the environment using the

Laplace Equation. Using this information, you will then solve the Poisson equation for pressure.

A schematic of the 2D environment is shown in Figure 1.

x

z

D

=

2

5

m

L = 25m

T(0, z) = T

0

T(x,0) = T0

T(L, z) = T

0

T(x,D) = f(x)

Figure 1: Schematic of the 2D domain.

1

1 Analytic Solution to the Steady-State Temperature Field (20%)

Taking x to be the horizontal position and z to be the depth (positive downwards) for the 2D

domain, the boundary temperature distributions are defined as (Figure 1):

T (0, z) = T (L, z) = T (x, 0) = T0 (1)

T (x,D) = f(x) = T0 + T1 sin(pix/L) + T2 sin(5pix/L) (2)

where L = 25m is the width of the environment in the x-direction, D = 25m is the depth of the

environment in the z-direction, and the temperatures are T0 = 20

◦C, T1 = 380◦C and T2 = 205◦C.

Given that the temperature field is at steady-state we can use the Laplace equation:

∂2T

∂x2

+

∂2T

∂z2

= 0 (3)

1. Assuming that the temperature distribution has the form T (x, z) = F (x)G(z), use separation

of variables to show that:

F (x) = A cos(px) +B sin(px) (4)

G(z) = Cepz +De−pz (5)

(5%)

2. By recognising that T0 appears on all four boundaries, derive the steady-state solution as:

T (x, z) = T0 +

∞∑

n=1

Bn sin

(npix

L

)

sinh

(npiz

L

)

(6)

Show your working and justify any assumptions you make. (10%)

3. Explain why Fourier coefficient solutions only exist for n = 1 and n = 5, and hence state the

exact analytic solution for the steady-state temperature field. (5%)

2 Numerical Solution to the Temperature Field (40%)

In this section you will solve the temperature field numerically using two different methods: (i)

Gauss-Seidel and (ii) Successive Over-Relaxation (SOR). We can discretise the spatial domain into

nodes along x and z according to:

xi = (i− 1)∆x where i ∈ [1, nx] (7)

zj = (j − 1)∆z where j ∈ [1, nz] (8)

and denote the temperature at position (xi, zj) by Ti,j . Because we are considering a steady-state

problem, time is not a variable, however we still need to iterate our numerical solution and so will

use the superscript n to index the current solution iteration. In other words, Tni,j represents the

numerical temperature estimate at location (xi, zj) for iteration n.

The Gauss-Seidel (G-S) method computes successive grid point values according to:

Tn+1i,j =

1

4

(

Tni+1,j + T

n+1

i−1,j + T

n

i,j+1 + T

n+1

i,j−1

)

(9)

2

Here, provided that you maintain an array Tni,j for the current iteration and a second array T

n+1

i,j

for the next iteration, you can use this scheme to sequentially compute all the component values

Tn+1i,j . This concept will be illustrated in the lecture.

The Successive Over-Relaxation (SOR) method is a modified form of the G-S method which

computes successive grid point values according to:

Tn+1i,j = T

n

i,j + 1.85

(

1

4

(

Tni+1,j + T

n+1

i−1,j + T

n

i,j+1 + T

n+1

i,j−1

)

− Tni,j

)

(10)

Using this information, complete the following tasks:

1. Using the information above, explain how you will implement these numerical methods in

MATLAB. In your explanation, you should either include a control flow diagram or steps

with pseudo code. Make sure to outline how you will:

• Implement the boundary conditions

• Initialise the solution domain

• Structure your code to iterate and update the solution

• Check for convergence

(10%)

2. Solve the steady-state temperature field for the given boundary conditions using both the G-S

and SOR methods. Choose the number of divisions in the domain to be ndiv = 2

n for n = 2 : 9

in both the x and z directions. (Note: this means the number of points is npts = ndiv + 1.)

For each method and each domain discretisation, compute the number of iterations taken to

reach the steady-state solution within a tolerance of = 1× 10−6. Present a figure showing

the number of iterations vs grid size for the two methods. (15%)

3. Using the SOR method, produce three heat-maps showing the solution resolution improve-

ment at 3 different grid sizes tested above (choose a coarse, medium and fine grid to illustrate

the changes). (5%)

Tip: In order to produce a heat-map in MATLAB, use the surf() function with edge

colour set to zero. The syntax for this is: surf(X,Y,T,’EdgeColor’,’none’) where: X is a

matrix of x positions; Y is a matrix of y positions; and T is your 2D temperature solution.

Use the command view(0,-90) to project the surface into the x, y plane, with the colour

gradient representing the temperature in that cell. The command colormap() allows you to

change the colour map, see https://au.mathworks.com/help/matlab/ref/colormap.html

for more options.

4. Compute the L1, L2 and L∞ error norms and the order of the error for each respective norm

for the SOR method (you don’t need to do this for the G-S method). For this error norm

analysis, use the same number of points tested in Part 2.2. The error you compute should be

based on a comparison with the analytic solution obtained in Section 1. Present your results

in a table, as in the lecture slides. State the number of divisions in the domain required to

achieve asymptotic convergence. (10%)

The Lp norm for a two-dimensional solution can be computed using the expression:

Lp =

∆x∆z nx∑

i=1

nz∑

j=1

|Ti,j − T (xi, zj)|p

1/p (11)

3

Note when computing the norm in MATLAB, we need to supply the norm() function with

a vector input. You can use the following code to compute the p-norms for your solution:

% Create a variable for the difference between solutions:

diff = T numeric(2:end-1,2:end-1) - T analytic(2:end-1,2:end-1);

% "Unwrap" the 2D array to a 1D row vector;

diffVec = reshape(diff',1,numel(diff));

% Compute Lp norms as usual;

L1norm = norm(diffVec,1)*(dx*dz);

L2norm = norm(diffVec,2)*sqrt(dx*dz);

LInfnorm = norm(diffVec,Inf);

3 Poisson Equation for Pressure-Density Coupling (20%)

In this section you will extend the applicability of your SOR solver from Part 2 to the Poisson

Equation. Ultimately, we want to compute the density field, however the density field depends on

both the temperature (which you computed in Part 2) and hydrostatic pressure. The hydrostatic

pressure field can be calculated by solving a Poisson equation:

∂2P

∂2x

+

∂2P

∂2z

=

1

ρ

(

∂P

∂z

∂ρ

∂z

+

∂P

∂x

∂ρ

∂x

)

= φ(x, z) (12)

where ρ is the density and P is the pressure. Given the pressure field from the above equation,

and temperature field from the previous section, the density field can be computed using the Ideal

Gas law:

ρ =

P

RspecificT

(13)

where Rspecific = 287.058Jkg

−1K−1 is the specific gas constant, and T is the temperature (specified

in Kelvin for this equation).

To solve the Poisson equation, the SOR method may be adapted to include the source term φ(x, z)

as follows:

Pn+1i,j = P

n

i,j + 1.85

(

φi,j +

1

4

(

Pni+1,j + P

n+1

i−1,j + P

n

i,j+1 + P

n+1

i,j−1

)

− Pni,j

)

(14)

You can take the boundary conditions of the pressure field to be:

P (x, z = 0) = 101325Pa (15)

∂P (x = 0, z)

∂x

=

∂P (x = L, z)

∂x

= 0 (16)

∂P (x, z = D)

∂z

= −ρ(x, z = D)g (17)

where gravitational acceleration g = 9.81m/s2. These boundary conditions can be discretised as

follows:

Pni,j=1 = 101325 (18)

Pni=1,j = P

n

i=2,j (19)

Pni=nx,j = P

n

i=nx−1,j (20)

Pni,j=nz = P

n

i,j=nz−1 + ρ

n

i,j=nz−1g∆z (21)

4

Note that the boundary conditions for the density field must be computed using the Ideal Gas

Law.

Conceptually, the sequence of steps you will follow is:

i To start your simulation, assume the pressure distribution is P = 101325Pa everywhere.

ii Use the initial pressure field together with the computed temperature field to obtain an initial

guess for the density field.

iii With the initial guess for pressure and density, you can now solve for the new value of pressure

after application of the SOR method:

• Compute the derivatives for the source term.

• Using your modified SOR solver, compute the new pressure field Pn+1i,j .

iv Use the new pressure field to update the density field at the computed points (Equation 13).

v Update the boundary point values of pressure based on the new pressure field.

vi Update the density at the boundary points based on the new pressure at the boundary.

vii Repeat steps (iii) to (vi) until the residuals are below a specified tolerance.

Note that a useful check that your solver is giving sensible results is to consider that for a constant

density field, the total pressure difference from top to bottom of the domain should be equal to

ρgD.

Using this information, complete the following tasks:

1. State the discretisation you have used to calculate the source term φ. (5%)

2. Using the information gained in Section 2, and your own grid convergence study for this

problem, determine an appropriate grid size and residual convergence criteria which gives a

grid converged solution. Present two line plots of density and pressure across the mid-height

of the domain (z=12.5 m) for at least three grid resolutions to support your choice. (10%)

3. Using your solver, compute the density and pressure fields for the domain. Produce heat-

maps showing:

(a) Initial guess for the density field;

(b) Final estimate for the density field (i.e. once pressure field has reached convergence);

(c) Final estimate for the pressure field.

(5%)

5

Assignment Submission Information

Your assignment is to be presented as a concise report in PDF format. Note that 10% of the

assignment marks are awarded based on the presentation, clarity and functionality of your MAT-

LAB code. A further 10% of the assignment marks are allocated based on overall structure, clarity

and presentation of the report. An electronic copy of the report should be submitted to

Turnitin by the due date. Late submissions will incur a penalty of 5% per day late.

IMPORTANT:

• Marks will be deducted for handwritten answers and screenshots of equations and/or figures.

See Tutorial 2 for detailed notes on how to insert figures into your assignment.

• The report should not exceed 10 pages; additional pages will not be marked, so aim to be

concise with your findings.

• Any figure/table you include in your report must be numbered and must be referred to and

discussed in the main text of your report.

• You are strongly encouraged to read through the ”Example Assignment” we have made

available on Canvas here. This demonstrates the expectations with respect to assignment

layout, explanations and figures/tables.

MATLAB Submission

• You will submit your MATLAB Code independently from your assignment submission, not as

an appendix. Instructions on how to do this will be included in the Assignment 2 Information

page on Canvas.

• You are also required to submit a script, which should solve Section 2.2 for ndiv = 64 using

the SOR method, and return a heat-map of the temperature solution. This code will be run

by the marker and checked as a part of your MATLAB mark.

6

学霸联盟

Due: 11:59 pm Friday 4th June 2021

Introduction

A team is using drones to capture high resolution images of the terrain after recent volcanic activity

in Iceland. However, the images are distorted due to variations in air pressure and density caused

by the heating from the lava. The team has provided you with some data you can use to obtain a

model of the pressure and density fields in the environment - which will then allow them to apply

corrections to their drone images.

In this assignment you will use the analytical and numerical techniques developed in this course

to obtain a steady-state temperature field in a 2D rectangular slice of the environment using the

Laplace Equation. Using this information, you will then solve the Poisson equation for pressure.

A schematic of the 2D environment is shown in Figure 1.

x

z

D

=

2

5

m

L = 25m

T(0, z) = T

0

T(x,0) = T0

T(L, z) = T

0

T(x,D) = f(x)

Figure 1: Schematic of the 2D domain.

1

1 Analytic Solution to the Steady-State Temperature Field (20%)

Taking x to be the horizontal position and z to be the depth (positive downwards) for the 2D

domain, the boundary temperature distributions are defined as (Figure 1):

T (0, z) = T (L, z) = T (x, 0) = T0 (1)

T (x,D) = f(x) = T0 + T1 sin(pix/L) + T2 sin(5pix/L) (2)

where L = 25m is the width of the environment in the x-direction, D = 25m is the depth of the

environment in the z-direction, and the temperatures are T0 = 20

◦C, T1 = 380◦C and T2 = 205◦C.

Given that the temperature field is at steady-state we can use the Laplace equation:

∂2T

∂x2

+

∂2T

∂z2

= 0 (3)

1. Assuming that the temperature distribution has the form T (x, z) = F (x)G(z), use separation

of variables to show that:

F (x) = A cos(px) +B sin(px) (4)

G(z) = Cepz +De−pz (5)

(5%)

2. By recognising that T0 appears on all four boundaries, derive the steady-state solution as:

T (x, z) = T0 +

∞∑

n=1

Bn sin

(npix

L

)

sinh

(npiz

L

)

(6)

Show your working and justify any assumptions you make. (10%)

3. Explain why Fourier coefficient solutions only exist for n = 1 and n = 5, and hence state the

exact analytic solution for the steady-state temperature field. (5%)

2 Numerical Solution to the Temperature Field (40%)

In this section you will solve the temperature field numerically using two different methods: (i)

Gauss-Seidel and (ii) Successive Over-Relaxation (SOR). We can discretise the spatial domain into

nodes along x and z according to:

xi = (i− 1)∆x where i ∈ [1, nx] (7)

zj = (j − 1)∆z where j ∈ [1, nz] (8)

and denote the temperature at position (xi, zj) by Ti,j . Because we are considering a steady-state

problem, time is not a variable, however we still need to iterate our numerical solution and so will

use the superscript n to index the current solution iteration. In other words, Tni,j represents the

numerical temperature estimate at location (xi, zj) for iteration n.

The Gauss-Seidel (G-S) method computes successive grid point values according to:

Tn+1i,j =

1

4

(

Tni+1,j + T

n+1

i−1,j + T

n

i,j+1 + T

n+1

i,j−1

)

(9)

2

Here, provided that you maintain an array Tni,j for the current iteration and a second array T

n+1

i,j

for the next iteration, you can use this scheme to sequentially compute all the component values

Tn+1i,j . This concept will be illustrated in the lecture.

The Successive Over-Relaxation (SOR) method is a modified form of the G-S method which

computes successive grid point values according to:

Tn+1i,j = T

n

i,j + 1.85

(

1

4

(

Tni+1,j + T

n+1

i−1,j + T

n

i,j+1 + T

n+1

i,j−1

)

− Tni,j

)

(10)

Using this information, complete the following tasks:

1. Using the information above, explain how you will implement these numerical methods in

MATLAB. In your explanation, you should either include a control flow diagram or steps

with pseudo code. Make sure to outline how you will:

• Implement the boundary conditions

• Initialise the solution domain

• Structure your code to iterate and update the solution

• Check for convergence

(10%)

2. Solve the steady-state temperature field for the given boundary conditions using both the G-S

and SOR methods. Choose the number of divisions in the domain to be ndiv = 2

n for n = 2 : 9

in both the x and z directions. (Note: this means the number of points is npts = ndiv + 1.)

For each method and each domain discretisation, compute the number of iterations taken to

reach the steady-state solution within a tolerance of = 1× 10−6. Present a figure showing

the number of iterations vs grid size for the two methods. (15%)

3. Using the SOR method, produce three heat-maps showing the solution resolution improve-

ment at 3 different grid sizes tested above (choose a coarse, medium and fine grid to illustrate

the changes). (5%)

Tip: In order to produce a heat-map in MATLAB, use the surf() function with edge

colour set to zero. The syntax for this is: surf(X,Y,T,’EdgeColor’,’none’) where: X is a

matrix of x positions; Y is a matrix of y positions; and T is your 2D temperature solution.

Use the command view(0,-90) to project the surface into the x, y plane, with the colour

gradient representing the temperature in that cell. The command colormap() allows you to

change the colour map, see https://au.mathworks.com/help/matlab/ref/colormap.html

for more options.

4. Compute the L1, L2 and L∞ error norms and the order of the error for each respective norm

for the SOR method (you don’t need to do this for the G-S method). For this error norm

analysis, use the same number of points tested in Part 2.2. The error you compute should be

based on a comparison with the analytic solution obtained in Section 1. Present your results

in a table, as in the lecture slides. State the number of divisions in the domain required to

achieve asymptotic convergence. (10%)

The Lp norm for a two-dimensional solution can be computed using the expression:

Lp =

∆x∆z nx∑

i=1

nz∑

j=1

|Ti,j − T (xi, zj)|p

1/p (11)

3

Note when computing the norm in MATLAB, we need to supply the norm() function with

a vector input. You can use the following code to compute the p-norms for your solution:

% Create a variable for the difference between solutions:

diff = T numeric(2:end-1,2:end-1) - T analytic(2:end-1,2:end-1);

% "Unwrap" the 2D array to a 1D row vector;

diffVec = reshape(diff',1,numel(diff));

% Compute Lp norms as usual;

L1norm = norm(diffVec,1)*(dx*dz);

L2norm = norm(diffVec,2)*sqrt(dx*dz);

LInfnorm = norm(diffVec,Inf);

3 Poisson Equation for Pressure-Density Coupling (20%)

In this section you will extend the applicability of your SOR solver from Part 2 to the Poisson

Equation. Ultimately, we want to compute the density field, however the density field depends on

both the temperature (which you computed in Part 2) and hydrostatic pressure. The hydrostatic

pressure field can be calculated by solving a Poisson equation:

∂2P

∂2x

+

∂2P

∂2z

=

1

ρ

(

∂P

∂z

∂ρ

∂z

+

∂P

∂x

∂ρ

∂x

)

= φ(x, z) (12)

where ρ is the density and P is the pressure. Given the pressure field from the above equation,

and temperature field from the previous section, the density field can be computed using the Ideal

Gas law:

ρ =

P

RspecificT

(13)

where Rspecific = 287.058Jkg

−1K−1 is the specific gas constant, and T is the temperature (specified

in Kelvin for this equation).

To solve the Poisson equation, the SOR method may be adapted to include the source term φ(x, z)

as follows:

Pn+1i,j = P

n

i,j + 1.85

(

φi,j +

1

4

(

Pni+1,j + P

n+1

i−1,j + P

n

i,j+1 + P

n+1

i,j−1

)

− Pni,j

)

(14)

You can take the boundary conditions of the pressure field to be:

P (x, z = 0) = 101325Pa (15)

∂P (x = 0, z)

∂x

=

∂P (x = L, z)

∂x

= 0 (16)

∂P (x, z = D)

∂z

= −ρ(x, z = D)g (17)

where gravitational acceleration g = 9.81m/s2. These boundary conditions can be discretised as

follows:

Pni,j=1 = 101325 (18)

Pni=1,j = P

n

i=2,j (19)

Pni=nx,j = P

n

i=nx−1,j (20)

Pni,j=nz = P

n

i,j=nz−1 + ρ

n

i,j=nz−1g∆z (21)

4

Note that the boundary conditions for the density field must be computed using the Ideal Gas

Law.

Conceptually, the sequence of steps you will follow is:

i To start your simulation, assume the pressure distribution is P = 101325Pa everywhere.

ii Use the initial pressure field together with the computed temperature field to obtain an initial

guess for the density field.

iii With the initial guess for pressure and density, you can now solve for the new value of pressure

after application of the SOR method:

• Compute the derivatives for the source term.

• Using your modified SOR solver, compute the new pressure field Pn+1i,j .

iv Use the new pressure field to update the density field at the computed points (Equation 13).

v Update the boundary point values of pressure based on the new pressure field.

vi Update the density at the boundary points based on the new pressure at the boundary.

vii Repeat steps (iii) to (vi) until the residuals are below a specified tolerance.

Note that a useful check that your solver is giving sensible results is to consider that for a constant

density field, the total pressure difference from top to bottom of the domain should be equal to

ρgD.

Using this information, complete the following tasks:

1. State the discretisation you have used to calculate the source term φ. (5%)

2. Using the information gained in Section 2, and your own grid convergence study for this

problem, determine an appropriate grid size and residual convergence criteria which gives a

grid converged solution. Present two line plots of density and pressure across the mid-height

of the domain (z=12.5 m) for at least three grid resolutions to support your choice. (10%)

3. Using your solver, compute the density and pressure fields for the domain. Produce heat-

maps showing:

(a) Initial guess for the density field;

(b) Final estimate for the density field (i.e. once pressure field has reached convergence);

(c) Final estimate for the pressure field.

(5%)

5

Assignment Submission Information

Your assignment is to be presented as a concise report in PDF format. Note that 10% of the

assignment marks are awarded based on the presentation, clarity and functionality of your MAT-

LAB code. A further 10% of the assignment marks are allocated based on overall structure, clarity

and presentation of the report. An electronic copy of the report should be submitted to

Turnitin by the due date. Late submissions will incur a penalty of 5% per day late.

IMPORTANT:

• Marks will be deducted for handwritten answers and screenshots of equations and/or figures.

See Tutorial 2 for detailed notes on how to insert figures into your assignment.

• The report should not exceed 10 pages; additional pages will not be marked, so aim to be

concise with your findings.

• Any figure/table you include in your report must be numbered and must be referred to and

discussed in the main text of your report.

• You are strongly encouraged to read through the ”Example Assignment” we have made

available on Canvas here. This demonstrates the expectations with respect to assignment

layout, explanations and figures/tables.

MATLAB Submission

• You will submit your MATLAB Code independently from your assignment submission, not as

an appendix. Instructions on how to do this will be included in the Assignment 2 Information

page on Canvas.

• You are also required to submit a script, which should solve Section 2.2 for ndiv = 64 using

the SOR method, and return a heat-map of the temperature solution. This code will be run

by the marker and checked as a part of your MATLAB mark.

6

学霸联盟