Method for estimating physical characteristics of two materials
10788347 ยท 2020-09-29
Assignee
- United States Of America As Represented By The Secretary Of The Air Force (Wright-Patterson AFB, OH)
Inventors
Cpc classification
G06T11/008
PHYSICS
G01F1/74
PHYSICS
International classification
G01F1/74
PHYSICS
G01F1/64
PHYSICS
Abstract
An estimate of interfacial areas between the liquid and gas, liquid and wall, and gas and wall in two phase flow is determined using a standard 2D sensor in a fashion to infer 3D information about the liquid/vapor profile when the sensor length is much longer than the diameter. Cross-sectional flow areas for the gas and liquid are also estimated as a function of the axial dimension of the sensor, and the centroid of the mass in the sensor element can also be estimated. An electric capacitance tomography (ECT) system creates tomograms of the flow inside a sensor, and estimates of 3D physical area information are produced from the tomograms.
Claims
1. A method of estimating a physical characteristic of a material, comprising: disposing a flowing mixture of a first material and a second material in a three-dimensional sensor volume; defining within the three-dimensional sensor volume a matrix of parallel voxels, each voxel having x, y and z dimensions; measuring at least one parameter of the mixture of the first and second material within the three-dimensional sensor volume and producing a tomogram in the form of a two-dimensional matrix that is defined within a perimeter and contains multiple values with one value being associated with each voxel, each value representing an amount of the first material in the associated voxel; calculating multiple hypothetical points having x, y and z coordinates in the three-dimensional sensor volume based on the multiple values of the tomogram, the z coordinate of each point being calculated from at least one value of the tomogram and the x and y coordinates of each point corresponding to at least the associated voxel of the at least one value of the tomogram; calculating at least one estimated physical characteristic of the first material based on the multiple points; and changing the flowing mixture of the first and second materials in response to at least one estimated physical characteristic of the first material to change a physical characteristic of the flowing mixture.
2. The method of claim 1 wherein calculating at least one estimated physical characteristic comprises calculating a surface area based on the multiple hypothetical points that is an estimate of the surface area of an interface between the first and second materials within the volume.
3. The method of claim 1 wherein the three-dimensional sensor volume is defined by a circumferential sensor wall and by length dimensions of a sensor and wherein calculating at least one estimated physical characteristic comprises calculating a hypothetical surface area of the interface between the circumferential wall and the first material based on the values in the tomogram along the perimeter of the tomogram.
4. The method of claim 1 wherein the three-dimensional sensor volume is defined by a circumferential sensor wall and dimensions of a sensor and wherein calculating at least one estimated physical characteristic comprises calculating a hypothetical surface area of an interface between the circumferential wall and the second material based on the multiple hypothetical points that are disposed adjacent to the circumferential wall.
5. The method of claim 1 wherein calculating at least one estimated physical characteristic comprises calculating a volumetric void fraction based on the multiple hypothetical points that is an estimate of a volume fraction of at least one of the materials within the three-dimensional sensor volume.
6. The method of claim 1 wherein calculating at least one estimated physical characteristic comprises providing mass densities of the first and second materials within the three-dimensional sensor volume and calculating an estimated centroid of mass of the mixture of the first and second material within the three-dimensional sensor volume based upon the multiple hypothetical points and the mass densities of the first and second materials.
7. The method of claim 1 wherein calculating at least one estimated physical characteristic comprises calculating an estimated cross sectional area of one of the materials within the three-dimensional sensor volume based upon the multiple hypothetical points.
8. The method of claim 1 wherein calculating multiple hypothetical points in the volume comprises: calculating the z coordinate of each point based on a value of the tomogram; and determining the x and y coordinates of each point based on the position of a voxel.
9. The method of claim 1 wherein calculating at least one estimated physical characteristic comprises: mapping a surface to the multiple points to produce a mapped surface; and calculating the at least one estimated physical characteristic of the first material based on the mapped surface.
10. The method of claim 9 wherein calculating the at least one estimated physical characteristic comprises calculating a surface area of the mapped surface to produce a calculated surface area that is an estimate of the surface area of an interface between the first and second materials within the volume.
11. The method of claim 9 wherein mapping comprises: defining a plurality of triangular planar surfaces using a plurality of combinations of three adjacent points within the multiple points whereby the mapped surface is the plurality of triangular planar surfaces; calculating the surface area of each triangular surface to produce multiple area values; and adding the multiple area values to produce the estimate of the surface area of the interface between the first and second materials within the volume.
12. The method of claim 9 wherein mapping comprises: defining a plurality of elements based upon the points, each element including four points; defining a three-dimensionally curved surface within each element; calculating the surface area of each three-dimensionally curved surface to produce multiple area values; and adding the multiple area values to produce the estimate of the surface area of the interface between the first and second materials within the volume.
13. The method of claim 9 wherein calculating the at least one estimated physical characteristic comprises calculating an estimated volume of the first material based on the mapped surface.
14. A method of estimating a physical characteristic of a material, comprising: providing a flowing mixture of a first phase of the material and a second phase of the material in a conduit wherein the flowing mixture is controlled by at least one pump and at least one valve; disposing the flowing mixture moving in a flow direction within a three-dimensional sensor volume defined by a circumferential sensor wall and length dimensions of a sensor; disposing a plurality of capacitive sensors in a side-by-side relationship around and adjacent to the circumferential sensor wall, each sensor having a width and a length, and the length of each sensor being parallel to the flow direction of the flowing mixture; defining within the three-dimensional sensor volume a matrix of parallel voxels, each voxel having x, y and z dimensions with the z dimension being parallel to the flow direction; measuring capacitance within the three-dimensional sensor volume using the plurality of capacitive sensors and producing a tomogram in the form of a two-dimensional matrix containing multiple values with one value being associated with each voxel and corresponding to the electrical permittivity of the material within the voxel that is associated with the value, each value also corresponding to an amount of the first phase of the material in the associated voxel; calculating multiple hypothetical points having x, y and z coordinates in the three-dimensional sensor volume based on the multiple values of the tomogram, the z coordinate of each point being calculated from at least one value of the tomogram and the x and y coordinates of each point corresponding to at least the associated voxel of the at least one value; calculating at least one estimated physical characteristic of the material based on the multiple points and based on the assumption that the multiple points lie on a hypothetical interface between the first and second phases within the three-dimensional sensor volume; and controlling at least one of the pump or the valve to modify the flowing mixture based upon the estimated physical characteristic.
15. The method of claim 14 wherein the estimated physical characteristic is at least one of the following: a surface area of the interface between the first and second phases, a surface area of the sensor wall in contact with the first phase, a surface area of the sensor wall in contact with the second phase, a mass centroid of the material within the three-dimensional sensor volume, the fraction of the three-dimensional sensor volume occupied by the first phase, and the cross sectional area of the first phase within the three-dimensional sensor volume.
16. The method of claim 14 wherein a surface area of the sensor wall in contact with the first phase is estimated using only the outermost values in the tomogram.
17. The method of claim 14 wherein a surface area of the sensor wall in contact with the first phase is estimated using only the outermost points of the multiple points.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
(1) The invention may best be understood by reference to various embodiments and variations of the invention, examples of which are described below in conjunction with the drawings in which:
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)
(10)
(11)
(12)
(13)
(14)
(15)
(16)
(17)
(18)
(19)
(20)
(21)
(22)
(23)
(24)
(25)
DETAILED DESCRIPTION
(26) Referring now to the drawings in which like reference characters refer to like or corresponding parts or elements throughout the several views, a technique is disclosed to estimate interfacial areas between two materials, such as a liquid and a gas. A technique is also disclosed for estimating interfacial areas between a material and a wall, such as between a liquid and wall and between a gas and a wall in two phase flow. The technique uses a standard 2D sensor in a fashion to infer 3D information about the liquid/vapor profile when the sensor length is much longer than the diameter. It will also allow the cross-sectional flow areas for the gas and liquid to be estimated as a function of the axial dimension of the sensor. The centroid of the mass in the sensor element can also be determined. Tomograms of the flow inside a sensor may, in some embodiments, be created by commercially available electric capacitance tomography (ECT) systems. The methods disclosed herein provide a quantitative interpretation of the tomogram providing estimates of 3D physical area information.
(27) Two phase flows consisting of a liquid and gas phase are common in many applications such as air conditioning, chemical and petroleum industries. Predicting pressure drop and heat transfer rates in these flows typically depends on knowledge of the void fraction which can be defined on an area or volumetric basis. In many cases the flow may be pulsating and chaotic which leads to difficulty in characterizing the interfacial areas between the tube wall and each phase and the interfacial area between each phase. The impact of the interfacial areas, A.sub.lw, A.sub.gw, A.sub.lg can be seen in the transient one dimensional momentum equations for separated flow given as: [1] (The bracketed numbers refer to references listed at the end of this Detailed Description.)
(28)
(The numbers enclosed within parentheses at the end of equations are equation numbers.) As may be observed in the equations, the transient liquid one dimensional momentum decreases in response to an increase in the area of the liquid/wall interface and increases in response to an increase in the area of the liquid/gas interface. Also, the transient gas one dimensional momentum decreases in response to increases in both the area of the gas/wall interface and the area of the liquid/gas interface. Thus, to understand and predict these momentums, it is important to know or estimate the aforementioned interfacial areas, A.sub.lw, A.sub.gw, A.sub.lg, and such information can be used separately or in conjunction with void fraction information to better predict physical characteristics of the flowing material, such as pressure drop and heat transfer rates.
(29) In this discussion, the following nomenclature is used:
Nomenclature
(30) a=x.sub.2x.sub.1 the width of a pixel (m)
(31) A.sub.i area of the ith pixel (m.sup.2)
(32) A.sub.gw interfacial area between gas and wall (m.sup.2)
(33) A.sub.l cross sectional flow area of the liquid (m.sup.2)
(34) A.sub.lg interfacial area between liquid and gas (m.sup.2)
(35) A.sub.lw interfacial area between liquid and wall (m.sup.2)
(36) A.sub.g cross sectional flow area of the gas (m.sup.2)
(37) ,
(38) b=y.sub.2y.sub.1 the height of a pixel (m)
(39) C.sub.i i.sup.th row in the connectivity matrix gives pixel number and the corner node numbers
(40) D diameter (m)
(41) ECT Electric capacitance tomography
(42) H(x) Heaviside step function given in eq. 27
(43) J.sub.i i.sup.th Jacobian given by eq. 7
(44) L sensor length (m)
(45) L.sub.liq length of a voxel that is occupied by liquid (m)
(46) m.sub.i mass of liquid and gas in the i.sup.th voxel (kg)
(47) n1.sub.i, n2.sub.i, n3.sub.i, n4.sub.i corner node numbers for the i.sup.th element
(48) N.sub.frame number of frames used in a temporal average
(49) NP number of active pixels used in a tomogram (e.g., 812 in the illustrated embodiment)
(50) NBP number of boundary elements (e.g., 88 in the illustrated embodiment)
(51) P pressure (Pa)
(52) R radius of tube (3.5 mm)
(53) S.sub.i area of the i.sup.th element on the liquid/gas interface
(54) t time (s)
(55)
(56) x.sub.i x coordinate of the ith node (m)
(57) x.sub.ci, y.sub.ci, z.sub.ci coordinates of the centroid of the mass in the i.sup.th voxel
(58)
(59) z axial coordinate (m)
(60)
(61)
(62) Subscripts
(63) g gas
(64) l liquid
(65) w wall
(66) Greek Variables
(67) specific weight (N/m.sup.3)
(68) x, y lengths of the sides of the square elements (m)
(69) .sub.i angle subtended by the ith area element of wall wetted by liquid (radians)
(70) relative permittivity
(71) .sub.* normalized relative permittivity ratio
(72)
(73) <
(74) angle (radians)
(75) vertical direction in master element
(76) horizontal direction in master element
(77) density (kg/m.sup.3)
(78) shear stress (Pa)
(79) , .sub.max angle and maximum value of the angle used in
(80) .sub.ij .sub.j.sub.i;=x, y, z notation used in eq. 14a.
(81) .sub.i i.sup.th shape function used in the interpolation functions given by eq. 8
(82) Many techniques have been used to estimate the void fraction including optical, gamma ray attenuation, and techniques based on either electric resistance or capacitance. The present approach uses an electric capacitance tomographic (ECT) technique to estimate the liquid/vapor interface in flows that may have different physical characteristics as shown in
(83) Some embodiments of the invention may use an ECT sensor 32 from Industrial Tomographic Systems [2].
(84)
(85) .sub.g is the relative permittivity of the gas phase
(86) .sub.l is the relative permittivity of the liquid phase
(87) is the relative permittivity measured at a given pixel in the tomogram.
(88) Thus a value of .sub.*=0 corresponds to gas and .sub.*=1 corresponds to a liquid.
(89) A representative tomogram 42 of a two phase mixture of the refrigerant R134a is shown in
(90)
NP=812 is the number of active pixels used in the tomogram
The spatial value is a function of time and the temporal average given as
(91)
(92) Here N.sub.frame is the number of frames or tomograms that are used in the time average. The void fraction is then estimated as 1<
(93) However, this information can be used in a new method to estimate the liquid profile and thus the interfacial areas as well as the centroids of the gas and liquid regions. It also can be used to calculate the estimated volumetric void fraction of the mixture in the sensor volume.
(94) Liquid/Vapor Interfacial Area
(95) To estimate the liquid profile in the sensor volume, the new method assumes that the pixel value represents the volume fraction of liquid in a rectangular voxel bounded by the pixel area times the length of the sensor as shown in
(96)
(97) To estimate surface areas, the values from a tomogram, such as tomogram 80 shown in
(98) This procedure is discussed in more detail below in connection with
(99) The following discussion for an exemplary embodiment of the invention will use a square mesh, although the described methods may be applied to any non-uniform mesh (such as shown in
(100) In summary, the tomogram represents a first square mesh 89 with each square 90 in the mesh representing a pixel. The first step of the disclosed method is to create a second square mesh where the centroids of the first square mesh form the corners of each square element 94 in the second mesh. The second mesh is not fully shown in
(101) A connectivity matrix is created which lists the elements 94 and the four corner node numbers going clockwise around the element 94, C.sub.i=(i, n1.sub.i, n2.sub.i, n3.sub.i, n4.sub.i), where i is the element number and n1.sub.i, n2.sub.i, n3.sub.i, n4.sub.i are the node numbers for the element.
(102) A matrix corresponding to the x,y,z values of the nodes is also created, where the nodes are the corners of the elements 94, which are also the centroids of the pixels 90. The z value is the measured permittivity ratio from the tomogram times the length 50 of the sensor. The surface elements 94 are thus created and are then mapped to a master element 100 as shown in
(103)
where a=x.sub.2x.sub.1; b=y.sub.2y.sub.1 are the lengths of the sides of the pixels 90. The area of the i.sup.th element is then given as: [4]
S.sub.i=.sub.1.sup.1.sub.1.sup.1{square root over (J.sub.1.sup.2+J.sub.2.sup.2+J.sub.3.sup.2dd)}(6)
J.sub.i are the Jacobians and depend on the mapping function and the element geometry. The Jacobians are given by:
(104)
(105) For the chosen geometry and nodes, first order interpolation functions are used to describe the coordinates in terms of the transformed variables and the nodal coordinates. Higher order functions can be described but additional nodes would be needed for each element. Thus,
(106)
The x.sub.i, y.sub.i, z.sub.i values are the coordinates for the i.sup.th node.
(107) The Jacobians can now be found from eq. 7. This will hold for the uniform mesh shown in
(108) As illustrated in
(109)
(110) This can be shown to be equivalent to using an expression from Larson et al page 1032. [5] The expression for the surface is given as z=f(x,y) over a square region R given by x.sub.1xx.sub.2; y.sub.1yy.sub.2. The area is given as:
(111)
(112) Here the same mapping to a master element 100 and linear interpolation functions are used and the interface written as z(,). The partial derivatives are evaluated using the chain rule and the area becomes:
(113)
(114) Comparing eq. 11 with eq. 6 after eq. 7 and 8 have been substituted it can be seen to be the same for the special case of square elements. The area of a given element 94 is found using Gauss-Legendre quadrature. [6]
(115)
(116) The integration reduces to the sum of 4 integrand evaluations at the locations given in equation 13b. For the uniform mesh, only differences in z at the nodes need to be calculated for each element. The other terms are the same for all elements and only need to be calculated once. The total area is the sum of the areas of the elements 94.
(117) Alternative Method to Finding the Area
(118) The above described method of mapping a surface to the points 82 of
(119)
The differences in the x and y directions are constant for the uniform mesh and would only need to be evaluated once. For a non-uniform mesh they would need to be calculated for each element.
Liquid Wetted Area
(120) Another surface area that may be calculated using the pixel data from the tomogram 80 is the surface area of the interface between the liquid (or the gas) and the wall of the tube 38 within the sensor 32 (
(121)
(122) The liquid/wall interface surface area 118 is found by integration around the tube perimeter and is demonstrated using the trapezoidal rule although other methods could be used.
(123)
NBP is the number of boundary pixels 110. From the uniform
(124)
For the last segment the first boundary pixel is used for z.sub.i+1. This completes the circle around the tube.
(125) With the uniform mesh, the arc lengths corresponding to the segments defined by the boundary pixels 110 aren't the same. In some meshes .sub.i may be a constant but the following approach will find the appropriate value using the boundary pixel 110 locations. The pixels 110 are first sorted in a clockwise or counter clockwise fashion. The location of the centroid (e.g. centroids 122a and 122b) of each pixel represents a vector (e.g., vectors 124a and 124b) extending radially from the center of the tube 38 to the pixel as shown in
(126)
and the area of the wall of tube 38 wetted by liquid is approximated as:
(127)
The area of the tube in contact with vapor (gas) is just the total tube area minus the liquid wetted area or
(128)
In practice, A.sub.wg would calculated as A.sub.wg=DLA.sub.wl.
(129) The error in using this technique has two components. The first is associated with using the simple trapezoidal rule for the integration as opposed to a more complex formula. The second is approximating the values of .sub.*,i at the tube wall as the values given by the centroids of the boundary pixels 110. The error associated with using the trapezoid rule can be estimated by looking at a case where the liquid/vapor interface forms a plane 126 that intersects the sensor volume 128 as shown in
(130) Volumetric Void Fraction
(131) It is recognized that the information taken from a tomogram such as tomogram 42 of
V=.sub.A.sup.1f(x,y)dxdy=.sub.A.sup.1f(x,y)J.sub.3dd(19)
The volume under a surface f(x.y) is illustrated in
Using the same bilinear mapping for the elements 94 as that used to find the liquid/vapor interfacial area the volume under an element 94 is given as:
(132)
As with the area calculation, this can be written for a non-uniform mesh with higher order interpolation functions as well. Evaluating eq. 20 gives
(133)
(134) This is just the average value,
(135) Centroid of the Fluid Mixture in the Sensor Volume
(136) The centroid of the fluid mixture can also be found from the tomogram data. Referring to
.sub.i=.sub.g+.sub.*,i(.sub.l.sub.g)(22)
The coordinates for the centroid of the mass in the sensor volume is given by
(137)
To find
(138)
Substituting into 23c gives
(139)
Cross Sectional Flow Areas
(140) The cross-sectional flow area 150 for the gas, A.sub.g, and the cross sectional flow area 152 for the liquid, A.sub.l, phases are shown in
(141)
Eq. 27 assumes the area of each pixel is the same. If they are different the expression can be expressed as:
(142)
(143) Here A.sub.i is the area of the i.sup.th pixel and would be calculated once for a given mesh using Gauss-Legendre quadrature. These areas can be can be retained in this discrete form or curve fit to provide a smoother approximation to the area. This also allows an approximation of the gradients in the axial direction of the flow areas either in a finite or continuous fashion. The cross sectional area of a particular material (e.g. liquid) in the sensor (such as sensor 32) will vary depending on the z position within the sensor. One method of estimating an average cross sectional area of a material in the flow would be to determine the cross sectional area at multiple z positions within a sensor. For example, the cross sectional area of the liquid could be calculated at every 0.1 inch along the length of the sensor and then those areas could be averaged to determine an average cross sectional area.
(144) The cross sectional area of a material at one z position in the sensor may also provide useful information. For example, the cross sectional area of a liquid in a gas/liquid flow may be calculated for the center of the sensor and such area may be calculated repeatedly over time. Each of the calculated areas may be compared to the others. If the cross sectional area of the liquid is fluctuating or oscillating over time, the frequency and amplitude of the oscillations in the cross sectional area would be a measure of the frequency and magnitude of waves in the flowing mixture within the sensor.
(145) Test Results
(146) The concept was tested using a vertical tube 160 with salt 162 and air 164 as two phase media in a vertical orientation as shown in
(147)
Exemplary Apparatus
(148)
(149) The sensor 180 provides an output to a controller 200 which includes data processors and communication devices for implementing the methods. The controller 200 powers the sensor 180 and receives communications from the sensor 180 through lines 202, which are communication lines and power lines. The controller 200 is also connected to power and control pump 192, valve 198, mixing valve 176 and pump 172 through the lines 202. The pump 192 is connected to a supply conduit 190 and, in this example, is supplied with water. The output of pump 192 flows through conduit 194 to valve 198 as indicated by the flow arrow 196, and the valve 198 controls the flow through conduit 199 to the mixing valve 176.
(150) The sensor 180 measures capacitance and those capacitance measurements are provided to the controller 200 which calculates a tomogram as discussed above with respect to sensor 32. The controller 200 repetitively samples the sensor 180 and repetitively produces tomograms at a rate that is sufficient for a particular application, which will vary widely. In this application, the controller 200 is configured to produce tomograms at a rate of one sample per second. The controller 200 is also configured to calculate one or more of the hypothetical physical characteristics discussed above in less than one half a second. So, for example, the controller 200 may calculate a hypothetical surface area of the interface between the water and steam within the sensor 180. In addition, the sensor 180 may calculate additional hypothetical physical characteristics, such as the hypothetical area of a wall of the sensor 180 in contact with steam. Then, the controller 200 compares the hypothetical physical characteristics against predefined limits and transmits control commands when the hypothetical physical characteristics meet or exceed the predefined limits. So, for example if the hypothetical surface area between the wall of sensor 180 and steam exceeds its predefined limit, the controller 200 issues control commands to the pump 192 and the valve 198 causing a desired amount of flow through the conduit 199 and water is introduced through the mixing valve 176 into the conduit 178. The supply of water through the mixer 176 will decrease the amount of steam in the flowing mixture and will decrease the surface area of the sensor wall that is contacted by steam.
(151) The controller 200 may also be calculating the hypothetical surface area between the water and gas within the sensor 180. Also, it may be saving each such calculation and calculating a rate of change in the hypothetical surface area between the water and gas. When this rate of change exceeds a predefined limit, that circumstance in this particular embodiment can be predicting the formation of oscillations within the flowing mixture in the conduit 178. In this particular embodiment, such waves would constitute a dangerous or catastrophic event. Thus, the controller 200 in response to such condition issues commands to stop the pump 172 and 192. In addition, it will command the valve 198 and the mixing valve 176 to stop all flow through the conduit 178, and the apparatus 170 is shut down. Alternatively, when the rate of change exceeds a predefined limit, the controller 200 may be programmed to take corrective action. For example, the pump 192, valve 198 and mixing valve 176 may receive commands to introduce more water into the flow within conduit 178. By increasing the water, hopefully, the rate of change in the surface area between the water and steam will reverse or stabilize. The controller 200 will continue to monitor such rate of change and will allow the embodiment to continue functioning so long as the rate of change remains below the predefined limit. It will be understood that all of the various hypothetical physical characteristics discussed herein may be calculated by the controller 200 and compared against one or more predefined limits, and in each case corrective actions may be executed when any of the hypothetical physical characteristics exceed their limits, and one of those corrective actions could be a complete shutdown of the apparatus 170.
(152) Having described several embodiments and variations of the invention in the above Detailed Description, it will be understood that the invention is capable of numerous modifications, rearrangements and substitutions of parts without departing from the spirit of the invention as defined in the Claims.
REFERENCES
(153) [1] Wallis, G., One-Dimensional Two Phase Flow, 1969, McGraw-Hill Inc. [2] Industrial Tomography Systems plc, Sunlight House, 85 Quay Street, Manchester, M3 3JZ, UK [3] Kreitzer, P., Hanchak, M. and Byrd, L., Horizontal Two Phase Flow Regime Identification: Comparison of Pressure Signature, ECT and High Speed Visualization, presented at 2012 ASME IMECE, Houston, Tex. [4] Taylor, A. E., Mann, W. R., Advanced Calculus, 2.sup.nd ed., 1972, Xerox College Publishing [5] Larson, R. E., Hostetler, R. P., Edwards, B. H., Calculus with Analytical Geometry 6.sup.th ed., 1998, Houghton Mifflin Co. [6] Carnahan, B., Luther, H. A., Wilkes, J. O., Applied Numerical Methods, 1969, J. Wiley & Sons, Inc.