Mcd 4140 computing for engineers - you have been asked to

Question

QUESTION 1

Background

You have been asked to manage a project to install a pipeline from an offshore gas platform to an onshore gas processing plant, and to find the cheapest design that is possible. The platform is Q km offshore, and D km along the coastline from the processing plant. The pipeline will follow a straight line from the platform to the shore, and hit land at a point that is a distance ‘x' from that point on the shore that is closest to the platform. We assume that the coastline is a straight line. The layout is shown schematically in Fig 1.

Fig 1 A schematic of the gas pipeline from an offshore platform to an onshore processing facility.

Q1a

Company A has quoted the cost of laying the subsea part of the pipeline at $PSA per km (regardless of water depth) and the cost of laying onshore pipeline at $POA per km. Write a MATLAB function called CostA that calculates the total cost of laying the pipeline as a function of ‘x', the distances D and Q and the cost factors $ PSA and $ POA. Note: write the function so that ‘x' can be a vector - assume all other input parameters are scalars.

If Q = 50km, D = 100, PS_A = $3.0M per kilometer and POA = $1.85M, plot the cost of the pipeline (in units of $100M) for 0 <= x <= 100km for 101 equally spaced points.

Q1b

Write a MATLAB function called dCostAdx that calculates the derivative of the cost with respect to the distance x. On a new figure, plot the derivative of cost w.r.t. x for 0 <= x <= 100km for 101 equally spaced points. Assume the same parameters from Q1a.

Use the Modified secant method to determine the distance ‘x' that gives the lowest cost pipeline. Print the lowest cost (in $M) to the MATLAB command window to 3 decimal places with a suitable statement attached (i.e. DO NOT just write the cost). On a new line print the distance ‘x' (in km - it might not be a whole number) where the cost is lowest (again to 3 decimal places with a suitable statement). In your fprintf statements, include units.

(NOTES:
1. This is a minimization problem.
2. You should write your own modified secant code.
3. You should use a starting guess in modified secant of 20.0
4. You should use an absolute precision of 0.001 (i.e., not a % precision)
5. You CANNOT easily pass a function that has multiple input parameters into another function as an input parameter. However, if you know all of the input values except x (i.e. you know Q, D, PSA , POA) you can define an anonymous function that you can pass as an input parameter

As an example, if x is a vector of values, and MyFunc a function with 4 parameters (x,A,B,C) then you cannot do the following
plot(x,MyFunc(x,A,B,C))
However, if you set values for A, B, C, then you CAN do this A=val1
B=val2; C=val3;
f=@(x) MyFunc(x,A,B,C)
plot(x,f(x))
NOTE: Every time you change one of the values (A, B or C), then you must redefine the anonymous function

Q1c

Keeping D, PS_A and PO_A fixed as in Q1a (i.e. D = 100, PS_A = $3.0M and PO_A = $1.85M), vary Q over the range 1 <= Q <= 75 km (in increments of 1 km). In a new figure, with two vertically stacked sub-plots, plot the optimal value of ‘x' for each value of Q in the top subplot. In the second subplot below the first, plot the total cost (in units of $100M) as a function of Q.

Q1d

Company B has also given a quote, where the cost of laying the onshore pipeline is fixed at $POB per km and the cost per kilometre of the subsea pipeline depends on the water depth η (also specified in kilometres). The cost per kilometer is

PS = PSB (1 + αη).

Consider the case where the depth of the water in kilometres (η) increases linearly with distance (in kilometres) from the shore (y) like
η (y) = εy

In this case, it is possible to calculate that the TOTAL cost of the subsea portion of this pipeline is

Cost_ss = PSB(1 + 0.5αεQ)(√(x² + (1 + ε²)Q²)

Write a function (called CostB) that calculates the TOTAL cost of the pipeline from Company B as a function of ‘x', the distances D and Q, cost factors PSB and POB, the slope factor ε and the depth-cost factor α.

If you assume the following parameters
• D=100km
• Q=50km
• PS_B = $2.1M/km
• PO_B = $1.5M/km
• α = 0.5
• ε = 0.05

In a new figure, plot the cost of the pipeline for 0 <= x <= 100km for 101 equally spaced points. ADD your plot from part 1a for comparison and add a legend to the figure. You will see the two curves cross at one value of x.

Q1e

Find the distance ‘x' at which company A and company B's quotes are equal. Check your answer is correct by calculating the cost using both Pipe Cost models and make sure they are the same. Print the distance (in km) and cost from both Companies (in $M) to the command window, accurate to 3 decimal places only (new line for each number/value you print).

NOTE: You must do this numerically, NOT visually. You will need to re-frame this question as a root finding problem. You can use fzero to find the root, or any other root finding method you like.

Q1f

You have just received an environmental report that shows that a reef is situated between the platform and the shore. This reef is the home of a rare species of red-lipped guppy and you are told that the only place you are allowed to bring the pipeline ashore is at a location given by x=65km. At this location the quote from Company B is more expensive.

You begin negotiating with Company B to drop their price. They refuse to change the cost per kilometer for the onshore section (i.e. POB = $1.5M/km) because they will subcontract this to another party. However you believe they might be prepared to change the cost per kilometer of the subsea section. The depth factor, α, is not negotiable, but the factor PSB might be. Currently they have quoted PSB = $2.1M/km. Find the value they need to drop PSB to in order to be competitive with Company A's quote for pipe landfall at x = 65km.

YOU MUST find this value automatically in MATLAB, not by trial and error. Again, you should frame this question as a root finding problem, and again you may use any root finding method you like (including fzero).

Print the competitive PSB value to the MATLAB command window (in $M/km to 3 decimal places).

Question 2

Background

A residence time distribution (or RTD) is a probability density function that describes how long fluid stays in a continuous flow chemical reactor. One way of measuring an RTD is to inject a small pulse of a chemical tracer (e.g. salt, radioactive material or coloured dye) at the inflow of the rector and then measure the signal at the outflow of the reactor.

At one extreme is a so-called "plug-flow" reactor which has no mixing and in which all of the input pulse of tracer exits the reactor at the same time. Provided there is no short-circuiting, this time is given (in seconds) as

τ= V/Q

where V is the volume of the reactor (in m3) and Q is the flow rate (in m³s^-1).

At the other extreme, a Continuously stirred tank reactor (CSTR) is one in which each element of fluid that is injected into the reactor is instantly uniformly mixed with everything else inside the reactor. The RTD for a CSTR is a negative exponential function.

Both of these extremes are idealised concepts and can never be realised in practice. In the real world, the RTD (usually) rises quickly, and then decays slowly, and in this question we investigate some different RTD's.

Q2a

You have been asked to determine if several different reactors are operating with similar behaviour, but you have not been given any information on their size, or design and have only been given a set of concentration measurements as a function of time, one for each reactor.

First you must open the rtd data file (‘rtd.dat') and read it into MATLAB using importdata. The first column of data is the time of the measurement (in seconds), and the other columns are the concentration measurements (C(t)) of the tracer at the exit of the reactor for an unknown number of reactors. Determine how many different reactors have been included in the file and print this number to the command window.

Plot each of the concentration versus time curves on the same figure using a different coloured line (in order, use as many as needed of black, red, green, blue, magenta, yellow). Ensure your plot has a legend using the text headers contained in the file.

Q2b

In order to compare the curves, they must be normalized. The normalized RTD curve (often called E(t)) is defined as

E(t) = C(t)/₀∫C^∞(t)dt

Write a function called CompTrap that calculates the integral of a function using the Composite Trapezoidal rule. The input parameters are a vector of (evenly spaced) times over the time range [t₁, t₂] and a vector of function values that corresponds to the time vector (you can assume that spacing of the data is uniform - you do not need to confirm this).

Normalise each of your concentration curves to give E(t) and in a new figure, plot these for each reactor using the same colours from part a. (NOTE: Instead of integrating to t = ∞ you should integrate to the last point in the data that you read in). Write the normalising value of ₀∫C^∞(t)dt for each reactor to the command window, one to a line.

Q2c

The mean residence time in the reactor can be determined from the following integral

τ_M = ₀∫^∞tE(t)dt

The mean residence time is also known as the "first moment" of the RTD. Higher order moments are also important and can be used to categorise and compare different reactors.

The second order moment is called the variance of the distribution and is

σ² = ₀∫^∞(t - τ_M)E(t)dt

It can be normalized by the mean residence time to provide a value that indicates how big the standard deviation is compared to the mean:

σ_N = √1/τ²_M0∫^∞(t - τ_M)²E(t)dt = (1/τ_M)√σ²

The third order moment is called the skewness of the distribution and is

s³ = 1/σ^3/2₀∫^∞(t-τ_M)³E(t)dt

Calculate τ_M, σN and s³ (i.e. equation 1, 2, 3) for each of the reactors using your CompTrap function. Print them to the command window in tabular form using fprintf (DO NOT USE THE MATLAB table function). The output should looks like the table below

Reactor	Mean RT	StDev_N	Skewness
1
2
etc.

If two reactors have similar values of both σN and s3 they are operating in a similar way, even if their mean residence times are quite different (for example, they might be the same design but have different sizes and different flow rates). For the purposes of this question, we will assume that two reactors are similar if their values of σN do not differ from each other by more than 10% AND if their s3 also vary from each other by less than 10%.

Based on the results in your table, automatically calculate if any two reactors are similar and print a sentence to the command window for each pair that are, saying which pair. If reactors 1 and 2 are similar, make sure you do not also print that reactors 2 and 1 are similar.