US 6807493 B2 Abstract A method of estimating liquid flow rates in capillary structures has application in determining the surface tension-driven flow rates of liquid fuel in propellant management devices in zero-gravity conditions. Analytic equations governing an assumed open-channel geometry are simplified by reasonable approximations, which allow a modified set of analytic equations to be derived. This modified set of analytic equations is derived by assuming an artificial taper in the capillary passage. Flow rate can be determined from the modified set of equations, which closely approximates the flow rate in the open-channel geometry.
Claims(32) 1. A method of estimating the flow rate of a liquid through an open-channel geometry, the method comprising steps of:
defining (i) an open-channel geometry having an inlet end and an outlet end, between which there is a capillary passage having two parallel walls, (ii) a set of parameters associated with the open-channel geometry, and (iii) analytical equations involving the set of parameters which govern the flow of liquid through the open-channel geometry; and
deriving a modified set of analytic equations corresponding to the governing analytic equations of the open-channel geometry, and solving said modified set of analytic equations to calculate a representative flow rate defined by the modified set of equations;
wherein said modified set of equations are derived from a modified open-channel geometry based on an artificial assumption that an elongate dimension of the open-channel geometry is tapered, so that the representative flow rate can be calculated as an approximation of the flow rate in the open-channel geometry.
2. The method as claimed in
3. The method as claimed in
_{1 }and outlet compartment volume V_{2 }associated with respective inlet and outlet compartments.4. The method as claimed in
5. The method as claimed in
_{1}, outlet liquid pressure P_{2 }and ambient pressure P_{0}.6. The method as claimed in
_{1}−P_{2 }by neglecting the pressure drops across the inlet and outlet compartments in the modified set of equations.7. A method as claimed in
_{i }and R_{o }of the liquid at the inlet and outlet ends of the open-channel geometry is assumed to be independent of the inlet or outlet end geometry, in the modified set of equations.8. The method as claimed in
_{1 }and R_{2 }at the inlet and outlet compartments of the open-channel geometry is calculated using the contact angle α and the inlet and outlet compartment volumes V_{1 }and V_{2}.9. The method as claimed in
_{1}−P_{2 }is determined from the Laplace-Young equation from the respective inlet and outlet meniscus radii R_{1 }and R_{2}, and the surface tension coefficient σ.10. The method as claimed in
11. The method as claimed in
12. The method as claimed in
13. The method as claimed in
_{i }and H/(L −X_{o})) are assumed to be vanishingly small, so that these terms can be disregarded in the modified set of equations.14. The method as claimed in
15. The method as claimed in
16. The method as claimed in
17. The method as claimed in
18. The method as claimed in
19. The method as claimed in
20. The method as claimed in
_{i }and H_{o }are given by: 21. The method as claimed in
22. The method as claimed in
_{i }and R_{o }of the free surface at the inlet and outlet ends of the modified open-channel geometry can be respectively calculated while maintaining the relationship between the contact angle α and the heights H_{i }and H_{o }at the inlet end and outlet end.23. The method as claimed in
_{i}−P_{o }by neglecting the velocity of the liquid through the open-channel geometry.24. The method as claimed in
25. The method as claimed in
26. The method as claimed in
27. The method as claimed in
28. The method as claimed in
29. The method as claimed in
30. The method as claimed in
31. An apparatus for estimating the flow rate of a liquid through an open-channel geometry, the apparatus comprising:
means for defining (i) an open-channel geometry having an inlet end and an outlet end, between which there is a capillary passage having two parallel walls, (ii) a set of parameters associated with the open-channel geometry, and (iii) analytical equations involving the set of parameters which govern the flow of liquid through the open-channel geometry; and
means for deriving a modified set of analytic equations corresponding to the governing analytic equations of the open-channel geometry, and solving said modified set of analytic equations to calculate a representative flow rate defined by the modified set of equations;
wherein said modified set of equations are derived from a modified open-channel geometry based on an artificial assumption that an elongate dimension of the open-channel geometry is tapered, so that the representative flow rate can be calculated as an approximation of the flow rate in the open-channel geometry.
32. A computer program product having a computer readable medium having a computer program recorded therein for estimating the flow rate of a liquid through an open-channel geometry, said computer program comprising:
code means for recording a definition of (i) an open-channel geometry having an inlet end and an outlet end, between which there is a capillary passage having two parallel walls, (ii) a set of parameters associated with the open-channel geometry, and (iii) analytical equations involving the set of parameters which govern the flow of liquid through the open-channel geometry; and
code means for storing a modified set of analytic equations corresponding to the governing analytic equations of the open-channel geometry, and solving said modified set of analytic equations to calculate a representative flow rate defined by the modified set of equations;
wherein said modified set of equations are derived from a modified open-channel geometry based on an artificial assumption that an elongate dimension of the open-channel geometry is tapered, so that the representative flow rate can be calculated as an approximation of the flow rate in the open-channel geometry.
Description The invention relates to estimating liquid flow rates in open-channel geometries and relates particularly, though not exclusively, to an improved method of effectively estimating liquid fuel flow rates and other parameters associated with the flow of liquid fuel through vane-type propellant management devices under zero-gravity conditions, that is, in fuel delivery systems used by space equipment. In vane-type surface tension propellant management devices (PMD) commonly used in satellite fuel tanks, the propellant is transported along guiding vanes from a reservoir at the inlet of the device to a sump at the outlet from where the fuel is pumped to the satellite engine. The pressure gradient driving this free-surface flow under zero-gravity (zero-g) conditions is generated by surface tension and is related to the differential curvatures of the propellant-gas interface at the inlet and outlet of the PMD. Constructional details of exemplary designs and systems can be found in the existing patent literature. In this respect, by way of example only, U.S. Pat. Nos. 4,733,531, 5,293,895, 4,553,565 and 4,743,278 disclose typical constructions. In a real-life situation, after one set of satellite orbiting manoeuvres are completed, the designer/operator would like to know for how long one has to wait before the sump gets refilled, so that the next set of manoeuvres may be performed. Given the complexity of the governing equations in this free-surface flow problem, the only possible options are to generate the required data for the drainage times either experimentally or via numerical simulation. If one decides to proceed experimentally, the options are to either perform the experiment in space (which would obviously be expensive) or on the ground, where one would be faced with the cost and difficulties of setting up a zero-g environment. Experimental approaches involve obvious difficulties and limitations, a notable one of which is expense. An overview of some of the other issues involved in using experimental techniques are discussed (in the context of space shuttle fuel tanks) in Dale A. Fester, Ashton J. Villars, and Preston E. Uney, Surface tension propellant acquisition system technology for space shuttle reaction control tanks. J. Spacecraft 9, 522 (1976). The second option, of using computer simulation, presents equally difficult challenges. One computes an unsteady, three-dimensional free-surface flow in a complex geometry. Further, the flow rates involved in a zero-g environment are extremely small; typically, the drainage time for 2 litres of propellant could be anywhere from 8 to 24 hours. These flow rates are several orders of magnitude less than those in a 1-g environment. Computing these zero-g flows by a brute-force unsteady 3-D simulation would not only be prohibitively expensive, but the issues of convergence and accuracy would be difficult to settle because of the extremely small flow rates involved. If one attempts to take care of the complexity of the geometry by ignoring the entry and exit regions, then the problem would be to accurately estimate the entry and exit conditions in a straight section of the vane. To the inventor's knowledge there do not currently exist any published numerical simulation results for a zero-g PMD. Estimates of drainage times, presumably obtained experimentally, are available in several references, typically papers of the American Institute of Aeronautics and Astronautics, such as: (a) A. Kerebel and D. Baralle, A low-cost surface tension tank optimised for telecommunication satellites. AIAA-85-1131; (b) A. Kerebel and P. Durgat, Development of a telecommunications spacecraft propellant tank. AIAA-86-1502; or (c) D. Baralle and J. P. Fournier, Propellant tank for telecommunication platforms. AIAA-89-2761. Owing to the substantial cost and the technical hurdles involved in accurately estimating these minuscule flow rates by either direct numerical simulation or by experimental methods which simulate zero-g conditions in the lab, any solution which offers advantages of any kind over existing techniques would provide benefits to those working in the field, particularly in relation to the design and general operation of satellite fuel tanks. In view of the above, a need clearly exists for an improved method of calculating flow rates in PMD devices that at least attempts to address one or more of the limitations of existing techniques. The aspects of the invention involve a recognition that an accurate and computationally tractable solution to the problem of calculating the fuel flow rates in a zero-g PMD can be achieved by using a semi-analytical procedure under certain reasonably idealized conditions. As the inventive technique uses exact analytical solutions (via eigenfunction expansions) of a suitably perturbed version of the Stokes flow equations, the issue of convergence is confined to the accuracy with which boundary data are satisfied. Further, due to relatively conservative computational and memory requirements, the effect of varying parameters (such as, for example, aspect ratio and contact angle) can be more readily investigated. Accordingly, a first aspect of the invention provides a method of estimating the flow rate of a liquid through an open-channel geometry, the method comprising steps of: defining (i) an open-channel geometry having an inlet end and an outlet end, between which there is a capillary passage having two parallel walls, (ii) a set of parameters associated with the open-channel geometry, and (iii) analytical equations involving the set of parameters which govern the flow of liquid through the open-channel geometry; and deriving a modified set of analytic equations corresponding to the governing analytic equations of the open-channel geometry, and solving said modified set of analytic equations to calculate a representative flow rate defined by the modified set of equations; In the modified set of equations, the taper of the channel geometry is assumed to be linear in the height of the open-channel geometry. The representative flow rate is calculated on the assumption that the flow rate is independent of the geometries of the inlet and outlet ends. Further, time derivatives are assumed to be negligible so that the flow rate is directly related to the instantaneous pressure gradient between the inlet end and the outlet end. The artificial assumption of a tapered channel height obviates the need to take into account the calculation of wall layers in the open-channel geometry for the analytic governing equations, so that the free surface of the liquid has approximately a constant radius of curvature which is able to satisfy requirements of a defined contact angle between the liquid and the open-channel geometry. The calculated meniscus radii of curvature of the liquid at the inlet and outlet ends of the open-channel geometry is assumed to be independent of the inlet or outlet end geometry, in the modified set of equations. The instantaneous pressure difference along the open-channel geometry is determined from the Laplace-Young equation and from the respective inlet and outlet meniscus radii and the surface tension coefficient of the liquid. A parameter δ defines a linear taper in average channel height H, such that the inlet and outlet heights H In the modified equations, the pressure difference for the open-channel geometry is equated with the pressure difference for the modified open-channel geometry to determine the parameter δ. The modified set of equations are preferably solved with a finite difference scheme using a weighted least-squares algorithm. The parallel walls of the open-channel geometry represent the fuel tank wall and a vane in a propellant management device, with the other two sides being open. The invention also provides a computer program product having a computer readable medium having a computer program recorded therein for estimating the flow rate of a liquid through an open-channel geometry, said computer program comprising: means for recording a definition of (i) an open-channel geometry having an inlet end and an outlet end, between which there is a capillary passage having two parallel walls, (ii) a set of parameters associated with the open-channel geometry, and (iii) analytical equations involving the set of parameters which govern the flow of liquid through the open-channel geometry; and means for storing a modified set of analytic equations corresponding to the governing analytic equations of the open-channel geometry, and solving said modified set of analytic equations to calculate a representative flow rate defined by the modified set of equations; The computer program is suitable for execution on existing general-purpose computing hardware. FIGS. 1A and 1B are side and end views of a schematic representation of a flow geometry and coordinate system used in conjunction with an embodiment of the invention. FIGS. 2A and 2B are side and end views of a revised schematic representation of a portion of the channel represented in FIG. FIG. 3 is an end view of a schematic representation of a flow geometry of FIG. 1 in an assumed static case of no fuel flow through the geometry. FIGS. 4A and 4B are side and end views of a schematic representation of a modified flow geometry in a linear gradient is assumed in the channel height. FIGS. 5A and 5B are side and end views of a non-dimensionalized schematic representation of the modified flow geometry of FIG. FIG. 6 is a flowchart of a procedure for calculating flow rate and drainage time in accordance with an embodiment of the invention. FIG. 7 is a schematic representation of computer hardware and associated equipment able to perform the calculation procedure of FIG. An embodiment of the invention is described in relation to the calculation of fuel flow rates in a PMD device as described below in further detail. The salient features of a PMD device are illustrated in the idealized geometry shown in FIGS. 1A and 1B, where all quantities in upper case are dimensional (length L, height H and width B). Here the reservoir Under zero gravity (zero-g) conditions, the pressure gradient driving the flow in the positive X direction is generated by surface tension and is related to the difference in the curvatures of the gas-propellant interfaces at the reservoir In FIGS. 1A and 1B, the pressure jumps discontinuously from the ambient P In the situation shown in FIGS. 1A and 1B, the contact angle α satisfy 0≦α<π/2, corresponding to a wetting fluid (which is the case for most propellants) and P The geometry of FIGS. 1A and 1B are an idealized version in the sense that (a) the vane is assumed to be straight in the flow direction (X) and (b) the inlet and exit compartments are simplified to be conical sections. The latter assumption is not essential for the described analysis and is made only for convenience. In conical sections, the meniscus curvature for a given volume of contained fuel and the rate of change of this curvature for a given drainage/accumulation rate of the contained fuel are both easily calculated analytically, while in a more complicated geometry these quantities would have to be obtained by numerical quadrature/differentiation. Neglect of the vane curvature is justified as follows. In practice, the fuel tank wall The crucial steps in idealization of the geometry are now described, which enable for the analysis. In actual practice, the length L of the channel is large, typically of the order of 100H. This suggests that the flow at locations far away from the entry and exit regions would be “fully developed”, that is, independent (to leading order) of the actual conditions prevalent at the inlet and outlet of the device. Thus the fully developed flow can be analyzed in a straight segment of length L* which is still assumed to be large compared to H, but whose entry and exit points are “far away” from the entry and exit regions of the actual device shown in FIGS. 1A and 1B in conclusion, the flow in the geometry of FIGS. 2A and 2B is to be analyzed, where
The flow rates are extremely small because of the small pressure differentials driving the flow. Therefore, as described hereinafter, the “quasi-steady” approximation is made: the time derivatives drop out to leading order in a suitably non-dimensionalized version of the governing equations. FIGS. 2A and 2B represent these assumptions. The flow rate in the device of FIG. 2A depends only on the (instantaneous) pressure gradient, which itself changes slowly with time. This pressure gradient may be obtained as follows. Let the radii of curvature of the meniscus at the inlet and outlet compartments be R Given the volumes of fluid at the inlet/outlet compartments and the contact angle α (which is a thermodynamic parameter dependent only on the propellant and the material of the PMD, usually titanium or stainless steel), R where R Because the velocity of the fluid is small, upon using Bernoulli's equation, the pressure drops across the inlet/exit compartments can be neglected and thus the pressure drop across the channel of length L in FIG. 1A is P The entry/exit regions are assumed to be small compared to the length of the device and may therefore be neglected; that is, the flow is “fully developed” over most of the length L. This is an excellent approximation because the length L is large compared to H and the flow rates/pressure gradients are small. The approximation L*≈L can be made without significant loss of accuracy and take the (instantaneous) leading-order pressures at the inlet and outlet of the device of FIG. 2 to be P The problem now is to (i) figure out the flow rate at this instant; (ii) to repeat this calculation for various other pressure drops till the inlet/outlet pressures in FIG. 1 equalize (that is, the meniscus radii satisfy R To gain some insight, first consider the static case when the inlet and outlet pressures in FIG. 2A are equal and there is no flow, as represented in FIG. Now consider the situation in which there is a small steady pressure gradient and therefore a small flow rate in FIG. Intuitively, the free surface shape still tends to be “nearly circular”; in other words, the pressure does not vary much across most of the cross section in FIG. 3, at any given location X. Further, assuming that the solid surfaces are highly polished, the contact angle must still be maintained. The only way to reconcile these two requirements is to assume that the free surface is circular (to leading order) everywhere except for small layers near the walls, where the contact angle requirement would force a deviation. The radius of curvature of this circular part at any instant of time is be determined by the corresponding leading-order pressure at any given location X and the Laplace-Young equation. In general, the analysis of the wall-layers is cumbersome analytically. The question that naturally arises is whether the (instantaneous) flow rates can be estimated without having to deal with these layers. If the contact angle condition is ignored and it is assumed that the wall layers are non-existent (and therefore that the free surface has uniform curvature throughout the cross-section), the flow rate becomes indeterminate, as demonstrated below. Instead, to deal with these conflicting requirements, the “fully developed” instantaneous flow rates are calculated for a modified geometry in which a small linear gradient is introduced in the channel height H, as shown in FIG. This assumption obviates the wall layer issues and imposes a free surface that has everywhere, to leading order, a constant radius of curvature that does satisfy the contact angle requirement of (3). At the inlet and outlet of this modified device, the channel heights are given by Here δ is a parameter satisfying 0<δ≦1 and is chosen as follows. Note that δ is fixed so that the average channel height is H (which occurs at the middle section of the device because of the linear gradient). The radii of curvature of the free surface of the device of FIG. 4B at the inlet and outlet is obtained from (3): Assume that the curvature of the free surface in the flow direction (X) in the geometry of FIGS. 4A and 4B is negligibly small (i.e., the radii of curvature in the flow direction are much larger as compared to R Equating this result with the right hand side of (2) and upon using (4) and (5), a quadratic equation for δ is obtained, which is solved as follows: Note that δ does satisfy 0<δ≦1 as required. Typically it is expected that T→0 and δ→1, because H is very small as compared to the compartment meniscus radii R Consider the rationale for modifying the geometry of FIG. 2A to that in FIG. Since the free surface in the bulk of the cross-section is unaffected and the contact angle is also maintained, the computed flow in the modified geometry of FIG. 4A is expected to be close to that in FIG. In summary, the problem now is to determine the flow rate at each instant in the idealized geometry of FIGS. 4A and 4B given the volumes of fluid at the inlet and outlet compartments, the meniscus radii R Above the flow is assumed to be quasi-steady. As a result, the instantaneous time dependence of the flow equations are ignored. A “slowly varying” (with respect to time) pressure drop is being imposed on the device of FIGS. 4A and 4B; δ starts with a value close to 1 and very slowly reaches the value of 1, when the pressures at the inlet and outlet of the PMD in FIG. 1A get equalized. This assumption holds well in practice because the inlet and outlet compartments have much larger cross-section areas as compared to that of the channel. As the flow rates are extremely small, the temporal rates of change of the compartment levels and therefore of the meniscus radii and the corresponding fluid pressures are also extremely small. Secondly the length L of the device is assumed to be large as compared to the channel height H and all quantities in the “fully developed” region of FIG. 2A (and hence FIG. 4A) are also “slowly varying” in the flow direction, X. All lengths are non-dimensionalized with respect to H
All variables in the flow equations are assumed to be functions of εx in order to translate the above assumption into practice. The equations for the channel height and the (leading-order) free-surface shape in non-dimensional form are given by: FIGS. 5A and 5B illustrate the coordinate system and the geometry of FIGS. 4A and 4B in non-dimensional form. Equations (8) and (9) are in tune with the above assumptions, namely that a linear gradient in the channel height and a uniformly circular shape of the leading-order free surface (which maintains the contact angle at all x) are imposed. Now let the velocities, pressure and surface tension coefficient be non-dimensionalized by V*, P* and σ*, where Here μ is the fluid viscosity and (2) and (4) define the other quantities. Now let (u, v, w) be the non-dimensional velocity components in the (x, y, z) coordinate directions and let p and λ be the non-dimensional pressure and surface tension coefficient (λ=σ/σ*). Note that λ may be obtained from the following: The second equality follows from follows from (4)-(6). The above non-dimensionalization yields a Reynolds number (R The non-dimensional fluid pressure p may now be obtained from (1) as follows, upon neglecting the curvature of the free surface in the flow direction (which, from (9), may easily seen to be of higher order): where {circumflex over (p)} is a higher-order term satisfying {circumflex over (p)}(ξ, y, z)=O(ε). Here p
The inertial terms in the perturbed flow equation can be neglected. Assuming that the component velocities (u, v, w) are each functions of (ξ, y, z) and upon neglecting all higher-order terms (in particular, the time derivatives, the inertial terms and all second derivatives with respect to x) in the Navier-Stokes equations, the following equations are obtained to leading order: where and the continuity equation is given by
Here the subscripts denote partial derivatives. Solutions to (14) and (15) in which u is symmetric about the planes y=0 and z=0 and v and w satisfy the corresponding obvious symmetry conditions in y and z as indicated, for example, by (15). Only the domain in which y and z are non-negative is considered. These symmetry conditions are specified by:
At the solid wall y=ω, the no-slip and no-normal velocity boundary conditions must be satisfied:
At the free surface, the normal component of the vector velocity u [=(u, v, w)] must vanish. This may be stated to leading order as follows:
Here ∇ is the gradient vector (∂/∂x, ∂/∂y, ∂/∂z). In the ensuing section, the terms εu, v and w in the above equation are shown to be all O(ε
where τ Finally, note that the perturbation pressure {circumflex over (p)} (where {circumflex over (p)} represents the O(ε) component of the pressure indicated in (12)) induces a higher-order perturbation to the free surface that, in general, is of non-uniform curvature. The condition that the perturbed free surface maintain the contact angle along the contact line (y, z)=(ω, b/2) is also required to hold. This condition is somewhat complicated because of the fact that the leading-order free surface also exhibits an O(ε(1−δ)) deviation from the actual physical contact angle, which is a three-dimensional entity (defined as the angle between the normal to the free surface and the normal to the solid boundary at the contact line). The solution procedure is now outlined, skipping the obvious algebraic details. An attempt is made to provide qualitative justification for some of the steps in this procedure. Define In what follows the functional dependencies of various quantities (for example on ξ) will be assumed to be understood when not indicated. A solution u to (14) and (15) that satisfies (16) is assumed as follows:
where the two velocity vectors on the right hand side each satisfy (14) and (15) (but not necessarily (16)) and are three-dimensional (3-D) and planar respectively. The planar velocity (0, v First the three-dimensional solution for the flow rate calculation is described. The most general solution for u that satisfies (14), (16) and (17) is given by: where A where C where B
where D Having satisfied the continuity and momentum equations, the boundary conditions are examined. The no-normal velocity condition (17) is imposed on v and for j≧1, where and H The terms w Note that it cannot be taken for granted that Equation (28) is automatically satisfied even if the condition is not imposed, because the 3-D solution fails to satisfy the symmetry constraint v Here the expressions for P where and Equations (29) and (30) are critically important in the described solution scheme in the sense that their satisfaction guarantees the existence of the planar solution; this issue is discussed further below. The problem now is to solve (29) and (30) for the unknown functions A The 3-D component of the term {circumflex over (p)} in (12) may now be stated as: where C
where: Here F and E As described hereinafter, a weighted least-squares algorithm is used to implement (29) and (30). In the least-squares procedure the requirement R(ξ)=0 in (34) is included, with an extremely small weightage. The idea is that the deviations from R(ξ)=0 in the 3-D solution are kept at manageable levels but without compromising on (29) and (30); if this condition is not imposed, the flow rate is not significantly affected, but the stresses at y=ω especially near the contact line, become large. If the weightage given to R(ξ)=0 is sufficiently small, the performance of the algorithm with respect to (29) and (30) also actually improves. In general, the 3-D solution exhibits deviations from (19) that are O(ε) everywhere (except possibly at the contact line). Neglect of the boundary condition (19) to leading order is justified by the fact that normal stress at the free surface, (namely, the pressure to leading order), is O(1/(1−δ)), as is evident from (10) and (12); this dominates the tangential component, which is O(ε). Equation (19) may be taken care of at the next level of approximation, when the free surface is perturbed. This issue is not examined in detail. The basic idea is as follows: consider a class of perturbations to the free surface, of O(ε), such that the non-planar contribution to the tangential stress of the 3-D solution on the free surface at every cross-section in ξ gets cancelled out for every member of this class except possibly for a small layer near the contact line. The planar solution, described in the following section, determines the appropriate member of the above class such that (19) holds, except in the above-mentioned layer. In the layer near the contact line, it is, in general, not be possible to perturb the free surface as described above, unless slip is allowed (the so-called Navier slip condition) and the contact line becomes dynamic rather than static. This is so because of the constraint that the class of perturbed free surfaces must maintain the contact angle. Therefore Equation (19) does not hold at points on the free surface that are sufficiently close to the contact line, here assumed static. This is not surprising because the shear stress, in general, exhibit singular behaviour at the contact line. One interesting possibility is that the possible deviations from (19) near the contact line may be handled by permitting a non-linear gradient in the channel height ω(ξ), as follows:
Equation (8) may be thought of as a linearization of the above as δ→1which does not significantly influence the leading-order flow rate and this justifies the linear gradient imposed on the channel height in (8). For the planar solution to exist, the net (two-dimensional) flow into every cross-section must vanish. If (29) and (30) are satisfied by the 3-D solution (and this is assumed to be the case), the object is to show that the boundary conditions for the planar solution is compatible with this requirement. Denote by ψ(y, z; ξ) the stream function for the planar solution, with
Then ψ satisfies the biharmonic equation:
Here {overscore (ω)} is the vorticity of the planar solution. Let the free surface be perturbed as follows:
For a suitable initial guess ƒ
z=f _{p}(ξ, y): ψ=0, (τ_{2})_{t}=−(τ_{3})_{t}. (36)
Here (τ Clearly if the planar solution exists, the full solution indeed satisfies all the boundary conditions (except for (19) near the contact line). A crucial requirement for the existence of ψ is apparent from (36). Denote
To avoid a discontinuity in ψ at the corner point (y, z)=(0, ƒ It is easy to show that (37) is true provided the 3-D solution satisfies (29) and (30); this follows from an elementary mass balance on a segment of the channel of length dx in FIG. The numerical algorithm for the planar solutions is not described, as the implementation of such algorithms and techniques are well understood and widely practiced in the art. The perturbation solution may be further pursued by considering the O(ε) pressure term, which induces an O(ε A couple of important assumption are involved in the solution procedure outlined above. The first of these is that u Secondly, a family of perturbations of O(ε) are assumed to exist to the leading-order free surface, which maintains the contact angle to this order and for which the non-planar component of the 3-D tangential stress at the free surface vanishes everywhere except possibly near the contact line. The planar velocity component v The assumption made is that there does exist a member of the above family of perturbations which satisfies the above requirements. These assumptions are plausible, but obviously difficult to justify. The conclusion that there is significant viscous dissipation due to rapid variations in y and z in the planar solution near the free surface may be drawn. The dissipation significantly reduces the flow rate obtainable for the given average pressure gradient. Numerical solutions later described support this conclusion. It can be considered whether requiring the 3-D solution to satisfy (29) and (30) is the only way to ensure that the planar solution exists. This would appear to be extremely likely. An alternative method, if it exists, must work at every location in ξ. However, there are other possible solutions than the one proposed here which is discussed below. Consider the following modification of (23): where Ξ The justification for setting Ξ Let us now consider an alternative solution procedure that may at first sight look easier and more logical. Let w The question is whether Ξ This solution procedure fails for the special case in which Ξ Energy considerations are used to derive an upper bound for the flow rate. Consider the geometry of FIG. 2 in which the channel width is B, the height is H (with no gradient imposed) and with the instantaneous pressure difference over the length L being given by (2). Assume further that the “free surface” is in fact a solid boundary with slip allowed (but no normal velocity) and has the shape of a circular arc that maintains the contact angle. If δ is calculated as in (7) and using the same non-dimensionalization as before, the “fully developed” velocity profile in this device is:
The solution for u is obtained from (22) by setting A The flow rate in this hypothetical device, given by εΓ Therefore a significant portion of the energy supplied to maintain the pressure gradient is dissipated in sending the fluid along a spiralling path from inlet to outlet, as compared to the straight path in the hypothetical device. This dissipation due to cross-flow is not negligible, as explained above. Further, for a given contact angle, the ratio of the actual flow rate to this upper bound of εΓ The computed solution satisfies these two important checks (in addition to others), as later shown. Note also that the computed flow rate must be the maximum possible for the given pressure gradient, given that the cross flow effects have been minimized. Finally, an important question that arises is how small should the parameters ε and (1−δ) be in order for valid results to be obtained. One way to estimate this a posteriori is to check whether the perturbation pressure {circumflex over (p)} is sufficiently smaller than the leading term in (12). Since λ as defined in (10) is O(1/[1−δ]), the criterion for validity of these results may be stated as:
This equation must hold throughout the domain, but is particularly important at z=ƒ, otherwise the assumed circular shape of the leading-order free surface would be in error. Note that {circumflex over (p)} in (39) is given by
where {circumflex over (p)} Equation (39) can be checked a posteriori, after the coefficients A Before presenting the numerical algorithm below, the significant advantages gained with the described semi-analytical approach are considered. First consider am important invariance property of the equations (29) and (30). Multiply (30) by (1−δ) and consider the δ-dependent terms in (29) and (30). Observe that on the left hand side of these equations, the terms A
Every other function of ξ in (29) and (30) now becomes a “hatted” function of {circumflex over (ξ)} with no explicit δ-dependence. An explicit δ-dependent terms are easily seen to drop out of the transformed equations; the only δ-dependence of these “hatted” equations is via the fact that the range of the independent variable {circumflex over (ξ)} is now on {circumflex over (ξ)}ε[0, (1−δ)]. Once these transformed equations (and hence (29) and (30)) have been solved for one particular value of δ, say the initial (in time) value δ=δ The advantages gained with the described approach can now be more clearly appreciated. As seen directly above, if successfully solved for the flow rate at the initial instant, the flow rate for all instants is solved. Secondly it is required only to discretize the first order derivatives of A Further, a higher-order accurate finite difference scheme can be used without running into numerical instability. Note that the “fast variations” in y and z (especially near the free surface) have been taken care of analytically. Because the shape of the free surface is known to leading order, no iterative procedures are needed. Also because exact solutions of the equations are used, the only concern involves the accuracy to which the boundary conditions are satisfied, in order to be reasonably certain that the small flow rates are accurately computed. Contrast this approach with the standard “brute force” numerical simulation approach. A full-fledged three-dimensional time-dependent simulation in a complex geometry and to a very demanding accuracy (because of the small flow rates, of the order of 10 Therefore the three-dimensional grid must be fine near these regions; in view of the accuracy demand, the fine spatial discretizations needed would make the computation extremely costly in terms of both time and memory. Further, the time discretiztion must also be fine in view of the well-known stability limitations of such algorithms (via the “CFL” criterion); for this reason, a higher-order spatially accurate scheme would make the computation costly in time as well. The cost factor is compounded by the fact that the free surface must be determined iteratively. Even if these formidable obstacles are overcome, the accuracy of the computed flow cannot be assured except via prohibitively expensive grid refinement studies. In short, a full-fledged three-dimensional numerical simulation is fraught with difficulties and requires highly sophisticated numerical techniques implemented on a super-computer. When it is considered that with the present method, the total drainage time can be estimated in a single simulation (or maybe a few simulations, to optimize the numerical parameters) of a few seconds implemented on personal computer equipment as described below, the cost advantages alone are apparent. The numerical algorithm is now described in sufficient detail to instruct an implementation. A change of the dependent variables from A
Upon differentiating with respect to ξ, the following is obtained: The purpose of this transformation is to ensure that the coefficients of the matrix (in the system of linear equations to be solved) are all bounded as n→∞; this makes the matrix relatively better conditioned. Secondly, the finite difference discretization errors in Ã′ Upon substituting (41) in to (29), (30) and R=0 in (34), the following is obtained: where
In what follows, the tildes are dropped and the algorithm is described with respect to equations (29), (30) and (34). It is understood that the transformed equations, to which the algorithm equally applies, are actually being solved in practice. The interval ξε[0, 1] is discretized into I equally spaced points (ξ The details of this scheme are given below. Use the notation ω A weighted least-squares algorithm is used to implement (29), (30) and (34). This algorithm is simple to describe: form the weighted sum {overscore (S)} of the squares of the errors in (29), (30) and (34) with weight factors θ The least-squares approach has been used successfully in recent times by Shankar (P. N. Shankar, The eddy structure in Stokes flow in a cavity. J. Fluid Mech. 250, 371 (1993); P. N. Shankar, Three-dimensional eddy structure in a cylindrical container. J. Fluid Mech. 342, 97 (1997)) to solve biharmonic and Stokes flow problems in two and three dimensions. Define and Here the partial derivatives are interpreted in the usual way, treating the IN+1 unknowns A
The least-squares algorithm for the discretization given below is as follows: for j=1,2, . . . , I and m=1,2, . . . , N. Equations (42) and (30) (evaluated at ξ=0 and multiplied by the weight factor θ
with Γ as the last unknown. The equations are ordered in similar fashion (with j and m taking the place of i and n in the above scheme), with (30) as the last equation. In addition to the change of variables mentioned at the beginning of this section, the matrix was pre-conditioned by normalizing the maximum magnitude in each column to 1. Some comments on this algorithm follow. The weight factors used are:
Setting θ The condition number test applied was that the machine should not give the same values for the condition number C and C+1, when implemented in double precision. Although this test did not fail in the examples presented, it is nevertheless recommended that either more sophisticated pre-conditioning methods be used or that this algorithm be implemented in quadruple precision. For a given value of N (typically 4-10 for the examples considered here) some experimentation is needed to fix the optimal values of I, J, K and M. Due to the high condition numbers, the convergence is extremely sensitive to the values of these parameters, as well as those of δ, α and the aspect ratio defined as B/H. It should be emphasized that the exact solution itself is not necessarily sensitive as indicated above. The point is that if the algorithm yields an acceptable solution for a certain combination of parameters, making very small changes in these parameters may tend to make the algorithm fail. A different combination of parameters that work (in the sense described in the ensuing paragraph) may need to be found. Typically, I ranged from 5-20; J from 30 to 120; K from 15 to 60; and M from 100 to 200 in the examples presented in the ensuing section. The aspect ratios considered ranged from 1 to 8; the contact angle α, from 0° to 60°; and δ was set at 0.7. The problem is solved for one value of δ, which may be chosen according to convenience. In practice, δ is close to 1. The computation times on a personal computer are of the order of a few seconds or minutes; the memory requirements of this algorithm are negligible. For estimating the convergence of the obtained solutions, the flow rates at various stations in ξ for a given combination of parameters are calculated and compared for equality. If the no-normal velocity boundary conditions at the free surface and the solid wall are more accurately satisfied, the discrepancy in the calculated flow rates at various ξ is lowered. An algorithm that computes these boundary conditions as well as the error in (37) was also implemented. Secondly, as N and/or I are increased, the flow rate should exhibit convergence. Finally, the finite difference scheme has a local truncation error of O([h(1−δ)] A special solution is considered that is of relevance here. Suppose the harmonic solution in (22) is to vanish identically; that is, A The leading-order flow rate Θ(ξ) (=εΓ) is as defined in (28). The flow rate must, of course, be independent of ξ by continuity. If the quadratic velocity u The addition of the harmonic correction term in (22) to this quadratic solution should therefore lead to a much better match of the flow rates at various stations than provided by the quadratic solution. This is one of the criteria used to estimate the convergence of the solution. Define the residual (i.e. the left hand side minus the right hand side) of (29), divided by ε For the quadratic and computed solutions, a comparison of this residual as well as that of the normal velocity v
Four values of the contact angle α, namely, 0°, 15°, 30° and 60°, are considered. It must be kept in mind that for a given propellant and the material of the PMD (usually titanium or stainless steel), α is a constant which is close to zero for most propellants. In the tables that follow, β is the aspect ratio B/H and ζ
where ψ The ratio υ=Θ/Θ Θ Θ Obviously the object is to demonstrate that Θ FIG. 6 summarizes the entire procedure for calculation of the drainage time. In step If the comparison of step In the next step In Table 1 directly below shows the results of a simulation for α=0° and β=2. Here the notation aEb means a×10
Note the excellent agreement of the computed flow rates at various stations in ξ against the relatively poor agreement for the quadratic flow rates. Observe that the residuals in the free surface boundary condition have been reduced by two orders of magnitude as compared to the quadratic case. Importantly, ζ
Note that the coefficients A
Note that the convergence is fairly satisfactory and that N=6, 7 produce the best results. In general, the number of Fourier modes (J and K) needs to be increased as N increases. The minimum number of modes required depends on the values of the coefficients A
Finally, some results for other aspect ratios and contact angles are presented in Table 5 below.
Note that υ<1 in all instances (refer also to the third and fourth tables) and for given α, υ decreases as β increases. These facts support the remarks made at the end of Sec. 4. An additional interesting phenomenon is worth pointing out. For β=2, as α is increased through (0°, 15°, 30°, 60°), it is found from the third and fifth tables that the respective values of υ are given by (0.557, 0.465, 0.424, 0.606). This behaviour of υ may be explained as follows. As the contact angle increases, the curvature of the free surface decreases and therefore by continuity, the dissipation due to cross-flow effects may be expected to decrease (in view of the slower variations in y and z). But as α increases, the total cross-section area increases and this factor tends to increase the dissipation due to cross flow because of the longer spiralling paths traced out by the fluid. Thus a minimum in υ is found at some positive valued α because of the influence of these two competing factors. Finally, the criterion for validity of the stated results of this section, namely Eq. (39), was checked at the free surface with {circumflex over (p)} replaced by {circumflex over (p)} _{n}, Q_{n }and RIn the preceding analysis, the expressions for P Initially, define (using appropriate formulae in the main text and omitting functional dependencies; in particular, ƒ is defined in Eq. (9)). The formulae for P In the above, a finite difference scheme for numerical computation can be used as now described. Consider five equally spaced points x where Ω In this case, if (ξ The technique described above is now summarised in overview with reference to FIG. Initially, an open-channel geometry is specified in steps Once the parameters of open-channel geometry and the fluid are appropriately specified, respective meniscus radii R At this stage, a determination is made as to whether ε In step Other terms R If the solution is determined to have converged, the flow rate is calculated in step Other variants of the described algorithm are possible. For example, still higher-order finite difference schemes can be considered and (38) could have a role to play for higher contact angles and/or aspect ratios. A problem of interest is to include the effects of small gravity as well in the present method. This is straightforward if gravity is in the flow direction. But if the gravity force is in any other direction, the symmetry of the calculated flows in y and z would be destroyed. In this case the method proposed is still applicable; the eigenfunctions and Fourier modes corresponding to other symmetries in y and z must also be included and boundary conditions on the full domain must be considered. The algebraic complications would increase considerably. The above described process can be implemented using a computer program product in conjunction with a computer system The computer system The computer The video interface The method steps for are effected by instructions in the software that are carried out by the computer system In particular, the software may be stored in a computer readable medium, including the storage device The computer system In some instances, the program may be supplied to the user encoded on a CD-ROM or a floppy disk (both generally depicted by the storage device Further to the above, the described methods can be realised in a centralised fashion in one computer system Computer program means or computer program in the present context mean any expression, in any language, code or notation, of a set of instructions intended to cause a system having an information processing capability to perform a particular function either directly or after either or both of the following: a) conversion to another language, code or notation or b) reproduction in a different material form. Other variants of the described algorithm are possible. For example, one could consider still higher-order finite difference schemes and (38) could have a role to play for higher contact angles and/or aspect ratios. A problem of interest is to include the effects of small gravity as well in the present method. This is straightforward if gravity is in the flow direction. But if the gravity force is in any other direction, the symmetry of the calculated flows in y and z would be destroyed. In this case the proposed method is still applicable; the eigenfunctions and Fourier modes corresponding to other symmetries in y and z must also be included and boundary conditions on the fill domain must be considered. The algebraic complications would increase considerably. The techniques described are more broadly applicable to other applications besides propellant fuel flow through propellant management devices. In particular, the described embodiment is directly applicable to problems which involve estimating transfer flow rates for any liquid by capillary pumping in zero-gravity conditions (that is, in “space”). This may arise in the context of satellites, as described, or in relation to the delivery of liquid fuels in spacecraft from fuel tank to the spacecraft's thrusters or to the fuel tank of another space vehicle. It is understood that the invention is not limited to the embodiment described, but that various alterations and modifications, as would be apparent to one skilled in the art, are included within the scope of the invention. Patent Citations
Non-Patent Citations
Classifications
Legal Events
Rotate |