# Droplet Simulation

1 **
**1.1 Introduction

An important reason for the construction of the Range Crossing is to reduce the gradient of the road for trucks. You will analyse the data in ass4q1in.csv to calculate gradients of the road with and without

the viaduct. For this question, we are assuming that the viaduct follows a straight line between

its end-points and the x-coordinate has been reoriented to act along the length of the bridge.

1.2 Requirements**
**For this assessment item, you must perform hand calculations:

1. Assuming that there is no viaduct, estimate the value of the gradient at the second

measurement point using backward, forward and central differences.

2. Discuss which is the best result out of the backward, forward and central differences and

why.

3. Estimate the overall gradient of the viaduct.

You must also produce MATLAB code which:

4. Calculates backward, forward and central differences for for the length of the road,

assuming there is no viaduct. The calculations should be for each possible section.

5. Verifies the results from Requirement 4 using the results of Requirement 1.

6. Estimates the numerical error in the calculations of Requirement 4. Perform this using

numerical derivatives for as many data points as possible; at the ends of the data, use the

value closest to the end (e.g. for the 1st datum, use the second derivative calculated at

the 2nd datum, and the third derivative calculated at the 3rd datum). For data spaced

evenly in *x*, the necessary formulae are [the methods are accurate]:

- Plots the results from Requirements 4 and 6. Also plots together with the central difference from Requirement 4. On the plots of the gradient, also show the viaduct average gradient. Discusses the results, including which is the best method.

8. Has appropriate comments throughout.

2

2.1 Introduction

The resistance of a resistor which has a constant cross-sectional area A can be calculated using

the formula

Thermoelectricity is a phenomenon that relies on changing material properties: electric

currents produced by an applied voltage difference generate heat (which could replace chemical

means), or temperature differences generate potential differences (most commonly applied in

thermocouples). It is, however a very inefficient means of converting between heat and electricity

which is fundamentally due to the relevant material properties.

For your assignment, the following values are to be used for resistivities, radii and lengths.

The properties for resistor 1 are:

In this question, you will calculate the resistances of a number of different resistor configurations:

1. Resistor 1 with uniform cross-section.

2. Resistor 1 with radius *r*1 at the start and radius 2*r*1 at the end (the radius changes at a

constant rate between the two ends).

3. A resistor constructed by joining resistor 2 end-to-end with resistor 3 (the component

resistors each have a constant cross-section).

4. A resistor constructed by joining resistor 2 end-to-end with resistor 3 (the radius changes

at a constant rate along each of the two component resistors so that the radius is *r*1 where

they join).

**2.2 Requirements****
**For this assessment item, you must perform hand calculations:

1. Find the resistances of configurations 1- 4 analytically.

2. Use the trapezoidal method to determine the resistance of configurations 2 and 3 using

four intervals along the total length of each configuration. Validate the results with

Requirement 1.

You must also produce MATLAB code which:

3. Repeats Requirement 2 and verifies the code.

4. Uses the trapezoidal method to solve configurations 1- 4 using appropriate intervals.

5. Reports the values from Requirements 1 and 4.

6. Estimates, reports and discusses the overall numerical error in the calculations of each

configuration in Requirement 4.

7. Plots

*R*(

*x*) from Requirements 1 and 4.

8. Has appropriate comments throughout.

**3
**

**3.1 Introduction**

Your task is to simulate the droplet acceleration using the formulae and parameters supplied in Assignment 1, and starting with the initial condition supplied in ass1in.csv.

To obtain consistent results from your random number generation, you should initialise the

seed to a fixed value using rng(seed). For your assignment, the value is:

seed = 6*:*721

**3.2 Requirements**

For this assessment item, you must perform hand calculations:

1. Use Euler’s method to estimate the velocity for the 2nd data value in ass1in.csv using 5

equal steps; validate against that data value.

You must also produce MATLAB code which:

2. Repeats Requirement 1 to verify the method.

3. Simulates the entire duration of ass1in.csv using appropriate parameters for the following

methods:

(a) Euler’s method

(b) ode23

(c) ode45

Reports the timestep used for Euler’s method; if you use the default settings for the

MATLAB ode solvers, report this, otherwise report what values you used for any changed

settings. Graphically validates the simulation results with the data in ass1in.csv and

reports the outcome of the validation.

4. Uses Simulink to:

(a) Simulate up to the 5th data value.

(b) Simulate the entire duration of ass1in.csv.

Plots the results and reports the final velocity within Simulink. Uses MATLAB to validate the results with the data in ass1in.csv.

- Simulates droplet collisions and accelerations.

(a) Loads the data in ass2in.csv and multiplies the velocities by 50.

(b) Assigns each droplet an \age”. The age grows so that, on average, for each 0.5 mm

the droplet travels, it collides with another droplet. Initialise the age of each droplet

to be a uniformly-distributed random number between 0 mm and 1 mm.

(c) At each timestep, accelerates every droplet according to Euler’s method as in Requirement 3.

(d) For each timestep, increments every droplet’s age. If the age is greater than a

uniformly-distributed random number between 0 and 1 mm, then the droplet is

to collide after the new velocity has been calculated for the timestep.

(e) Pairs up droplets that are to collide and collide them according to the algorithm of

Assignment 2. Randomly assign droplets to be paired with each other. If there is

an odd number of droplets marked for collision, choose one to not collide. Reset the

age of the new, combined droplet according to the initialisation in Requirement 5b.

(f) Simulates up to the entire duration of ass1in.csv.

(g) Plots the mean velocity, standard deviation of velocity and total number of droplets

as functions of time.

(h) Reports to the Command Window the total number of droplets remaining at the end

of the simulation. If there is only a single droplet remaining, reports when the last

droplet collision occurred.

6. Has appropriate comments throughout.

**Solution**** **

**f1.m**** **

function y=f1(x)

L=0.01;

%function defined

y=1/((1+(x/L))^2);

end

**f2.m**

function y=f2(x)

r2=0.001379;

y=1/(pi*r2^2);

end** **

**f3.m**** **

function y=f3(x)

r3=0.0014714;

y=1/(pi*r3^2);

end** **

**f4.m**** **

function y=f4(x)

r1=0.001;

r3=0.0014714;

y=1/(pi*(((r3-r1)*x/r1)+r1)^2);

end** **

**f5.m**** **

function y=f5(x)

r1=0.001;

r2=0.001379;

y=1/(pi*(((r2-r1)*x/r1)+r1)^2);

end** **

**f11.m**** **

function y=f11(x,L)

%function defined

y=1/((1+(x/L))^2);

end

**req_1.m**** **

x=[0,227,427,624,800]; %values given in provided excel sheet

z=[0,-51,-10,-46,-5];

%backward difference

for i=2:length(z)

dz_dx_back(i-1)=(z(i)-z(i-1))/(x(i)-x(i-1));

end

%forward difference

for i=1:length(z)-1

dz_dx_for(i)=(z(i+1)-z(i))/(x(i+1)-x(i));

end

%central difference

for i=2:length(z)-1

dz_dx_cent(i-1)=(z(i+1)-z(i-1))/(x(i+1)-x(i-1));

end

%second derivative

for i=2:length(z)-1

d2z_dx2(i-1)=2*((x(i+1)-x(i))*z(i-1)-(x(i+1)-x(i-1)*z(i))+(x(i)-x(i-1))*z(i+1))/((x(i)-x(i-1))*(x(i+1)-x(i))*(x(i+1)-x(i-1)));

end

%third derivative

for i=3:length(z)-2

d3z_dx3(i-2)=6*((x(i+2)-x(i))*(x(i)-x(i-2))*(x(i+2)-x(i-2))*(((x(i+1)-x(i))^2)*z(i-1)-(x(i+1)-x(i-1))*(x(i+1)-2*x(i)+x(i-1))*z(i)-((x(i)-x(i-1))^2)*z(i+1))-(x(i+1)-x(i))*(x(i)-x(i-1))*(x(i+1)-x(i-1))*(((x(i+2)-x(i))^2)*z(i-2)-(x(i+2)-x(i-2))*(x(i+2)-2*x(i)+x(i+2))*z(i)-((x(i)-x(i-2))^2)*z(i+2)))/(((x(i+2)-x(i))*(x(i+2)-x(i-2))*(x(i)-x(i-2))*(x(i+1)-x(i))*(x(i+1)-x(i-1))*(x(i)-x(i-1))*((x(i+2)-x(i))*(x(i)-x(i-2))-(x(i+1)-x(i))*(x(i)-x(i-1)))));

end

%for verification, u can check both the values that are the hand calculated

%and MATLAB calculated

%error in forward diff

for i=3:length(z)-2

error_for(i-2)=((x(i+2)-x(i-2))/2)*d2z_dx2(i)+(((x(i+2)-x(i-2))^2)/6)*d3z_dx3(i-2);

end

%error in backward diff

for i=3:length(z)-2

error_back(i-2)=((x(i+2)-x(i-2))/2)*d2z_dx2(i)-(((x(i+2)-x(i-2))^2)/6)*d3z_dx3(i-2);

end

%error in forward diff

for i=3:length(z)-2

error_cent(i-2)=(((x(i+2)-x(i-2))^2)/6)*d3z_dx3(i-2);

end

%plots of derivatives

plot(dz_dx_back);

hold on;

plot(dz_dx_for);

plot(dz_dx_cent);

xlabel(‘index’);

ylabel(‘derivative value’);

title(‘Plots of derivatives’);

legend(‘Backward difference’, ‘Forward difference’,’Central difference’);

%plots of errors

figure();

plot(error_back);

hold on;

plot(error_for);

plot(error_cent);

xlabel(‘index’);

ylabel(‘error value’);

title(‘Plots of errors’);

legend(‘Backward difference’, ‘Forward difference’,’Central difference’);** **

**req_2_3.m**** **

rho=1.68*10^-8;

r1=0.001;

L=0.01;

x=linspace(0,L,4);

sum=0;

for i=1:length(x)-1

sum=sum+((x(i+1)-x(i))/2)*(f1(x(i+1))+f1(x(i)));

end

int1=((rho)/(pi*r1^2))*sum; %see in req 2 solution 2

%for configuration 3

rho2=51.207*10^-8;

r2=0.001379;

L2=0.0126;

rho3=88.366*10^-8;

r3=0.0014714;

L3=0.01728;

x2=linspace(0,L2,4);

x3=linspace(0,L3,4);

sum2=0;

for i=1:length(x2)-1

sum2=sum2+(rho2*L2)*((x2(i+1)-x2(i))/2)*(f2(x2(i+1))+f2(x2(i)));

end

sum3=0;

for i=1:length(x3)-1

sum3=sum3+(rho3*L3)*((x3(i+1)-x3(i))/2)*(f3(x3(i+1))+f3(x3(i)));

end

int=sum2+sum3;

**req2_4.m**** **

n=100;

%config 1

rho1=1.68*10^-8;

r1=0.001;

L1=0.01;

x1=linspace(0,L1,n);

sum1=0;

for i=1:length(x1)-1

sum1=sum1+rho1*L1*((x1(i+1)-x1(i))*(1/(pi*r1^2)));

end

R1=sum1;

%config 2, using same values of rho, L and r

sum2=0;

for i=1:length(x1)-1

sum2=sum2+((rho1)/(pi*r1^2))*((x1(i+1)-x1(i))/2)*(f1(x1(i+1))+f1(x1(i)));

end

R2=sum2; %see in req 2 solution 2

%for configuration 3

rho2=51.207*10^-8;

r2=0.001379;

L2=0.0126;

rho3=88.366*10^-8;

r3=0.0014714;

L3=0.01728;

x2=linspace(0,L2,n);

x3=linspace(0,L3,n);

sum3=0;

for i=1:length(x2)-1

sum3=sum3+(rho2*L2)*((x2(i+1)-x2(i))/2)*(f2(x2(i+1))+f2(x2(i)));

end

sum4=0;

for i=1:length(x3)-1

sum4=sum4+(rho3*L3)*((x3(i+1)-x3(i))/2)*(f3(x3(i+1))+f3(x3(i)));

end

R3=sum3+sum4;

%config 4

sum5=0;

for i=1:length(x2)-1

sum5=sum5+(rho2*L2)*((x2(i+1)-x2(i))/2)*(f4(x2(i+1))+f4(x2(i)));

end

sum6=0;

for i=1:length(x3)-1

sum6=sum6+(rho3*L3)*((x3(i+1)-x3(i))/2)*(f5(x3(i+1))+f5(x3(i)));

end

R4=sum6+sum5;

%in workspace see the values from R1 to R4

%for error

%absolute value of error is h^2 * (b-a)/12 *f”(x)

error1=(((x1(2)-x1(1))^2)*L1/12)*0; %since double derivative of constant is 0

error2=(((x1(2)-x1(1))^2)*L1/12)*(f1(3)-2*f1(2)+f1(1))/((x1(2)-x1(1))^2);

error3=(((x2(2)-x2(1))^2)*L2/12)*(f2(3)-2*f2(2)+f2(1))/((x2(2)-x2(1))^2)+(((x3(2)-x3(1))^2)*L3/12)*(f3(3)-2*f3(2)+f3(1))/((x3(2)-x3(1))^2);

error3=(((x2(2)-x2(1))^2)*L2/12)*(f4(3)-2*f4(2)+f4(1))/((x2(2)-x2(1))^2)+(((x3(2)-x3(1))^2)*L3/12)*(f5(3)-2*f5(2)+f5(1))/((x3(2)-x3(1))^2);;

%to compare the values of error, refer to the error values in the workspace** **

**req2_7.m**** **

%this is for resistance change with change in length, I’ll just use the

%req2_4.m script and just vary length values

%x has been replaced by z and L have been replaced by x

n=100;

j=1;

for x=0.005:0.001:0.03

%config 1

rho1=1.68*10^-8;

r1=0.001;

z1=linspace(0,x,n);

sum1=0;

for i=1:length(z1)-1

sum1=sum1+rho1*x*((z1(i+1)-z1(i))*(1/(pi*r1^2)));

end

R1(j)=sum1;

%config 2, using same values of rho, L and r

sum2=0;

for i=1:length(z1)-1

sum2=sum2+((rho1)/(pi*r1^2))*((z1(i+1)-z1(i))/2)*(f11(z1(i+1),x)+f11(z1(i),x));

end

R2(j)=sum2; %see in req 2 solution 2

%for configuration 3

rho2=51.207*10^-8;

r2=0.001379;

rho3=88.366*10^-8;

r3=0.0014714;

z2=linspace(0,x,n);

z3=linspace(0,x,n);

sum3=0;

for i=1:length(z2)-1

sum3=sum3+(rho2*x)*((z2(i+1)-z2(i))/2)*(f2(z2(i+1))+f2(z2(i)));

end

sum4=0;

for i=1:length(z3)-1

sum4=sum4+(rho3*x)*((z3(i+1)-z3(i))/2)*(f3(z3(i+1))+f3(z3(i)));

end

R3(j)=sum3+sum4;

%config 4

sum5=0;

for i=1:length(z2)-1

sum5=sum5+(rho2*x)*((z2(i+1)-z2(i))/2)*(f4(z2(i+1))+f4(z2(i)));

end

sum6=0;

for i=1:length(z3)-1

sum6=sum6+(rho3*x)*((z3(i+1)-z3(i))/2)*(f5(z3(i+1))+f5(z3(i)));

end

R4(j)=sum6+sum5;

j=j+1;

end

x=linspace(0.005,0.03,26);

plot(x,R1);

hold on;grid on;

plot(x,R2);

plot(x,R3);

plot(x,R4);

xlabel(‘length(in m)’);

ylabel(‘Resistance (in ohm)’);

title(‘Plot of R(x)’);

legend(‘R1’, ‘R2’, ‘R3’, ‘R4’);