**Chemical and Process Engineering Resources**

**Chemical and Process Engineering Resources**

## Solving Integral Equations in MS Excel

Nov 08 2010 01:10 PM | Steve Hall in Calculations and Tips Share this topic:TankVolume, as it appears in the online store, calculates the liquid volume of partially filled tanks, horizontal or vertical, with various types of heads.Â One of these, the torispherical head, requires that integral equations be solved.Â This article shows how the integrals are evaluated.

Torispherical heads are made up of two parts.Â The rounded end has a spherical shape.Â It has a radius that is described with the numerical value â€œfâ€ which is the ratio of the dish radius to the vessel shell diameter.

Torispherical heads are also called â€œF&Dâ€ or â€œflanged and dishedâ€.Â The standard F&D head has f = 1.0 and kD = 3 times the metal thickness of the head.Â ASME F&D heads require a dish radius no greater than the diameter (f <= 1.0) and knuckle radius no less than 6% of the diameter (k >= 0.06) or three times the metal thickness, whichever is greater.

Figure 1: Torispherical Head |

When calculating the volume of fluid in a torispherical head, the head is divided into three parts.Â The first is the lower portion occupied by the knuckle (liquid height 0 to h1 as seen in the diagram).Â The second part is the portion covered by the spherical dish section, which also includes the knuckle at the perimeter.Â The third portion is that which is entirely within the knuckle but above the spherical section.

Each of the three parts are calculated separately and they each require solving an integral.Â Iâ€™ll use the integral for the first part for an example in this article.

Eq. (1) |

where

Â Â Â Â Â Â Â Â Â Â Â k = knuckle radius

Â Â Â Â Â Â Â Â Â Â Â D = vessel diameter

Â Â Â Â Â Â Â Â Â Â Â h = fill height of the vessel

Â Â Â Â Â Â Â Â Â Â Â n = R â€“ kD + (k^{2}D^{2} â€“ x^{2})^0.5, where R = vessel radius (= D/2)

Â Â Â Â Â Â Â Â Â Â Â w = R â€“ h

It looks complicated, but solving in Excel is straightforward if it is done sequentially as described below.

Integrals that are too complex to express implicitly are solved using â€œnumerical methodsâ€.Â A number of numerical methods have been used for integrals with Simpsonâ€™s Rule being one of the best.Â I chose to use Simpsonâ€™s Rule for this problem.Â (Other methods include the trapezoidal rule, Riemann sums, Romberg integration, Gaussian quadratures and the Monte Carlo method).Â There are many webpages that give the derivation of Simpsonâ€™s Rule and explain it in detail; I wonâ€™t repeat that work here.

Simpsonâ€™s rule requires that the integral be broken into intervals.Â Since the integral is evaluating the area under a curve, from x=0 to x= (2kDh-h^{2})^.5, the method first calculates the maximum value for x, divides that by the number of intervals, and then evaluates the function for each value of x.Â In other words, if the maximum value of x was 10 and there were 10 intervals, then the function would be evaluated with x = 0, 1, 2, 3, 4, etc. Â Notice that the term called â€œnâ€ above includes x in its formula.

An even number of intervals is required.Â The more the better. Â I found through trial and error that a good number to use for this particular problem is 1000 intervals.

Implementation could be done in a tabular form on an Excel worksheet.Â However, it is much more elegant to solve Simpsonâ€™s Rule in a Visual Basic for Applications function subroutine. Â Listing 1 gives the function.Â This is located in a VBA Module within TankVolume.Â Each place the calculation is required, the function is called using a cell formula of the form:

= Toris_V1(k, D, h)

where the variables k, D, and h are as defined above.Â Since recalculation takes time, which can be noticable, I put the function call in a conditional statement so it is only used when the tank is horizontal with torispherical heads.Â Assume the variable â€œhead_typeâ€ refers to the type of head and a value of head_type=4 refers to torispherical, the conditional function call becomes:

= if(head_type=4,Toris_V1(k, D, h), 0)

In this case, the function is called only if the heads are torispherical.Â Otherwise, a value of â€œ0â€ is returned.Â This is perfectly fine since the result of this cell in the spreadsheet is only used for torispherical heads.

Function Toris_V1(k, D, H) As Double ' ' Integral solution using Simpson's Rule ' Dim interval As Double Dim i As Integer Dim n As Double, n1 As Double, n2 As Double Dim X As Double, xmax As Double Dim steps As Integer Dim func As Double, sumfunc As Double ' On Error GoTo err_TorisV1 steps = 1000 ' '----------------------------------------------------------- ' Procedure: ' evaluate the maximum value for x (minimum value=0) ' for each value of x from min to max, at each interval, 'Â Â calculate n 'Â Â calculate the function 'Â Â sum results according to Simpson's Rule ' calculate result and return ' ' Derived function (arcsin is not an intrinsic function) 'Â Â Â Arcsin(x) = Atn (x / Sqr(-x * x + 1)) '----------------------------------------------------------- ' ' maximum value of x xmax = (2 * k * D * H - H ^ 2) ^ 0.5 R = D / 2 w = R - H interval = xmax / steps sumfunc = 0 X = 0 n1 = R - k * D n2 = k ^ 2 * D ^ 2 For i = 0 To steps Â Â Â n = n1 + Sqr(n2 - X ^ 2) Â Â Â func = Sqr(n ^ 2 - w ^ 2) / n Â Â Â func = n ^ 2 * Atn(func / Sqr(-func * func + 1)) - w * Sqr(n ^ 2 - w ^ 2) Â Â Â sumfunc = sumfunc + func Â Â Â If i > 0 And i < steps Then Â Â Â Â Â Â Â If (i / 2) = Int(i / 2) Then sumfunc = sumfunc + func Else sumfunc = sumfunc + 3 * func Â Â Â End If Â Â Â X = X + interval Â Â Â If X > xmax Then X = xmax Next i ' Toris_V1 = sumfunc * interval / 3 Exit Function ' err_TorisV1: Toris_V1 = 0 ' End Function |

## 2 Comments