# 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.
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]:  1. 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.

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 r1 at the start and radius 2r1 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 r1 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.

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.
(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.

1. 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.

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’);