Let us define p as the capillary pressure needed to force water out of a certain sample of porous rock until it occupies only the fraction s (saturation) of the pore space. How p depends on s is an important question in soil science and in the characterization of aquifers and oil reservoirs.

Question

Let us define p as the capillary pressure needed to force water out of a certain sample of porous rock until it occupies only the fraction s (saturation) of the pore space. How p depends on s is an important question in soil science and in the characterization of aquifers and oil reservoirs. To measure p versus s, one spins in a centrifuge a water-filled cylindrical sample of the porous rock and measures the volume of water removed (spun out) as a function of the spinning rate. The average saturation $\bar{s}$ of water left in the sample at the spin rate w is given by

$\bar{s} = \frac{r_{1} + r_{2} }{2r_{2}p } \int_{0}^{p}{\frac{s(p^{'} )}{(1 – Bp^{'} /p)^{1/2} } dp^{'}} ,$ (9.4.72)

[latex]\bar{s} = \frac{r_{1} + r_{2} }{2r_{2}p } \int_{0}^{p}{\frac{s(p^{'} )}{(1 - Bp^{'} /p)^{1/2} } dp^{'}} ,[/latex] (9.4.72)

Accepted Answer

where [latex]r_{1}[/latex] and [latex]r_{2}[/latex] are the distance of the two ends of the cylinder from the axis of rotation in the centrifuge and

[latex]p = \frac{1}{2} \Delta ho w^{2} (r^{2}_{2} - r^{2}_{1} ) and B = 1 - \frac{r^{2}_{1}}{r^{2}_{2}} , [/latex] (9.4.73)

where [latex]\Delta ho[/latex] is the difference between water and air densities.
Equation (9.4.72) is a Volterra equation of the first kind. One varies p to vary [latex]\bar{s} [/latex] and solves Eq. (9.4.72) for [latex]s(p^{'})[/latex] versus [latex]p^{'}[/latex]. The problem is that experiments are accompanied by error and a Volterra equation of the first kind is ill conditioned (see Linz, 1982), which means that many solution techniques for solving integral equations are overly sensitive to error in the data ([latex]\bar{s}[/latex] and p). For example, if the data are available on a uniform mesh [latex]h = p_{i+1} - p_{i}, [/latex] i being the i th measurement or the i th setting of rotation rate w, then the integral in Eq. (9.4.72) can be approximated numerically by

[latex]\bar{s}_{i} = \frac{h(r_{1} + r_{2} )}{2r_{2} p_{i} } [\frac{s(p_{0} )}{2} + \sum\limits_{j = 1}^{i-1}{\frac{s(p_{j} )}{\sqrt{1 - B_{p_{j} }/p_{i} } }+ \frac{s(p_{i}) }{2\sqrt{1 - B} } } ][/latex] (9.4.74)

for i = 1 , . . . , M, where M is the total number of measurements made. This is a linear algebraic system that can be solved by forward substitution. However, since the left-hand side of Eq. (9.4.74) is multiplied by h, it follows that, upon grid refinement (smaller h), error in [latex]\overline{s_{i} }[/latex] , is amplified as l/h in the solution. The problem can be alleviated by using a least squares technique to smooth the data (Linz, 1982). Other techniques can also be used (Ayappa, 1989).

To illustrate the numerical problem in solving Eq. (9.4.74) when [latex]\overline{s_{i} }[/latex] contains error, consider the theoretical capillary pressure curve

[latex]s (p) = \begin{cases} 1, & 0 \lt p \leq 2,\ \frac{1.5}{p} + 0.25, & p\gt 2,\end{cases} [/latex] (9.4.75)

where p is in an appropriate set of units. Suppose [latex]r_{1} / r_{2} = 0.5.[/latex] Then [latex](r_{1} + r_{2})/ 2 r_{2} = 0.75[/latex] and [latex]B = 0.75. [/latex] From Eq. (9.4.72), it follows that

[latex]\overline{s} (p) =\frac{0.75}{p} \int_{0}^{p} \frac{dp^{'} }{\sqrt{1 - 0.75 p^{'} / p } } , 0 \lt p \leq 2,[/latex] (9.4.76)

and

[latex]\overline{s} (p) = \frac{0.75}{p} \int_{0}^{2} \frac{dp^{'}}{\sqrt{1 - 0.75 p^{'} / p }} + \frac{0.75}{p} \int_{2}^{p} \frac{1.5/p^{'} + 0.25}{\sqrt{1 - 0.75 p^{'} / p }} dp^{'} = \frac{7}{4} - \frac{3}{2} \sqrt{1 - \frac{1.5}{p} } + \frac{2.25}{p} \left( anh^{-1} \sqrt{1 - \frac{1.5}{p} } - anh^{-1} \frac{1}{2} ight) , p\gt 2,[/latex] (9.4.77)

and Eq. (9.4.74) becomes

[latex]\overline{s}_{i} = \frac{0.75}{i} \left[\frac{1}{2} + \sum\limits_{j = 1}^{i-1}{\frac{s_{j} }{\sqrt{1 - 0.75 j / i} } + s_{i} } ight] ,[/latex] (9.4.78)

Where [latex]\overline{s}_{i} = \overline{s} (ih)[/latex] and [latex]{s}_{i} = {s} (ih)[/latex]. We can write Eq. (9.4.78) in matrix form as

[latex]\overline{s} = As + b,[/latex] (9.4.79)

Where [latex]b_{i} = 0.75/2i[/latex] and

[latex]a_{ij} = \begin{cases} \frac{0.75}{i \sqrt{1 - 0.75i / j} } , & j \lt i,\ \frac{0.75}{2i \sqrt{1 - 0.75} } , & j = i, \ 0, & j \gt i,\end{cases}[/latex] (9.4.80)

for i = 1 , . . . , M. Suppose now that we generate "experimental" values for [latex]\overline{s} [/latex] by adding a fixed percentage of random fluctuation to the theoretical values from Eqs. (9.4.76) and (9.4,77). We can generate such data on a uniform mesh of values for p, resulting in the vector [latex]\overline{s} [/latex], which can subsequently be used in Eq. (9.4.79) to solve for s. A Mathematica program is provided in the Appendix for this example, where we show that even a fluctuation of 0.01% results in large deviations from the theoretical function in Eq. (9.4.75).

Question 9.4.2: Let us define p as the capillary pressure needed to force wa...

Related Answered Questions

Consider in L2 (-∞, ∞) the operator Kv = ∫-∞^∞ exp (-(x²+y²)/2) v (y) dy. (9.5.79) (a) Prove that K is self-adjoint. (b) Prove that K is completely continuous. (c) Find the eigenvectors of K. ...

Verified Answer:

Consider the equation du/dt = iKu, u(t = 0) = u0, (9.5.126) Where i = √-1 and K is a completely continuous, normal operator in H. Assume that the imaginary part of the eigenvalue λn goes as 1/n^1/2 and examine the asymptotic behavior of u(r). ...

Verified Answer:

Find the Volterra equation corresponding to the initial value problem d² y/d x² + exp (-x) d y / d x + x² y = x, (9.4.64) Where d y/ d x = 1 and y = 2 at x = 0. (9.4.65) ...

Verified Answer:

Solve the equation u + Ku = f, (9.3.49) where u, f ∈ L2 (-∞, ∞; exp(-t²)) and K = aa^† (9.3.50) with a(t) = t² and f(t) = t^4. (9.3.51) ...

Verified Answer:

Solve the equation of the form u + a1 b1^† u + a2 b2^† u = f, (9.3.43) given by u (t) + 1/4 sin (Πt) ∫-1^1 cos (Πs) ds + 1/4 sin (2Πt) ∫-1^1 (cos 2Πs) u (s) ds = t^2.(9.3.44) ...

Verified Answer:

Solve the equation Ku = f for u, f ∈ L2 (-l, 1), where K is the two-term dyadic K = P0 P1^† + β P1 P1^†, (9.2.50) β is a number, and Pl is a Legendre polynomial defined by Pl (x) = 1 / 2^l l! d^l / d x^l (x^2 – 1)^l. (9.2.51) ...

Verified Answer:

Solve the equation Ku = f, u, f ∈ L2{-∞, ∞), where ∫-∞^∞ exp (-(x-t)²)u(t)dt = exp (-x²/2). (9.2.36) ...

Verified Answer:

Solve the equation u – ab^† u = f, where u, f ∈ L2 (0, ∞). ...

Verified Answer:

Solve the equation Ku = f, where u, f ∈ L2 (R3), and R3 is an unbounded three-dimensional Euclidean vector space, i.e., r = x i + y j + z i, -∞ < x, y, z < ∞. In polar coordinates, r = r sin θ cos Φ i + r sin θ sin Φ j + r cos θ k, Where 0 < r < ∞, 0 < θ < Π, and 0 < Φ < 2Π. ...

Verified Answer:

Solve the equation u + Ku = f, where K is a dyad operator, i.e., K = ab†, and u, f ∈L2(0, 1). The equation is given by U(X) – 1/2 ∫0^1 x t u (t) d t = 5x/6. (9.2.31) ...

Verified Answer: