Ley de Darcy para flujo multifásico

Morris Muskat et al. ^{[1] [2]} desarrollaron las ecuaciones que rigen el flujo multifásico (una ecuación vectorial para cada fase del fluido ) en medios porosos como una generalización de la ecuación de Darcy (o ley de Darcy ) para el flujo de agua en medios porosos. Los medios porosos suelen ser rocas sedimentarias como rocas clásticas (principalmente arenisca ) o rocas carbonatadas .

${\displaystyle \mathbf {u} _{a}=-\mu _{a}^{-1}K_{ra}\mathbf {K} \cdot \left(\nabla P-\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

Las fases fluidas actuales son agua, petróleo y gas, y se representan por el subíndice a = w,o,g respectivamente. La aceleración gravitacional con dirección se representa como o o . Nótese que en ingeniería petrolera el sistema de coordenadas espaciales está orientado hacia la derecha con el eje z apuntando hacia abajo. La propiedad física que vincula las ecuaciones de flujo de las tres fases fluidas es la permeabilidad relativa de cada fase fluida y la presión. Esta propiedad del sistema fluido-roca (es decir, sistema agua-petróleo-gas-roca) es principalmente una función de las saturaciones del fluido , y está vinculada a la presión capilar y al proceso de flujo, lo que implica que está sujeto al efecto de histéresis . ${\displaystyle \mathbf {g} }$ ${\displaystyle g\nabla z}$ ${\displaystyle g{\mathbf {e} _{z}}{^{3}}}$

En 1940, MC Leverett ^[3] señaló que para incluir los efectos de la presión capilar en la ecuación de flujo, la presión debe depender de la fase. La ecuación de flujo se convierte entonces en

${\displaystyle \mathbf {u} _{a}=-\mu _{a}^{-1}K_{ra}\mathbf {K} \cdot \left(\nabla P_{a}-\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

Leverett también señaló que la presión capilar muestra efectos de histéresis significativos . Esto significa que la presión capilar para un proceso de drenaje es diferente de la presión capilar de un proceso de imbibición con las mismas fases de fluido. La histéresis no cambia la forma de la ecuación de flujo gobernante, pero aumenta (generalmente duplica) el número de ecuaciones constitutivas para las propiedades involucradas en la histéresis.

Entre 1951 y 1970, las computadoras comerciales entraron en escena para los cálculos científicos y de ingeniería y las simulaciones de modelos. La simulación por computadora del comportamiento dinámico de los yacimientos de petróleo pronto se convirtió en un objetivo para la industria petrolera, pero la capacidad de procesamiento era muy limitada en ese momento.

Con un poder computacional débil, los modelos de yacimientos eran correspondientemente burdos, pero la ampliación de escala de los parámetros estáticos era bastante simple y compensaba parcialmente la burdosidad. La cuestión de la ampliación de escala de las curvas de permeabilidad relativa de las curvas de roca derivadas a escala de tapón de núcleo (que a menudo se denota la microescala) a las celdas de cuadrícula gruesas de los modelos de yacimiento (que a menudo se llama la macroescala) es mucho más difícil, y se convirtió en un campo de investigación importante que todavía está en curso. Pero el progreso en la ampliación de escala fue lento, y no fue hasta 1990-2000 que la dependencia direccional de la permeabilidad relativa y la necesidad de representación tensorial se demostró claramente, ^{[4] [5]} a pesar de que al menos un método capaz ^[6] ya se desarrolló en 1975. Un caso de aumento de escala de este tipo es un yacimiento inclinado donde el agua (y el gas) se segregarán verticalmente en relación con el petróleo además del movimiento horizontal. El tamaño vertical de una celda de cuadrícula también suele ser mucho menor que el tamaño horizontal de una celda de cuadrícula, creando áreas de flujo pequeñas y grandes respectivamente. Todo esto requiere diferentes curvas de permeabilidad relativa para las direcciones x y z. Las heterogeneidades geológicas en los yacimientos, como láminas o estructuras de permeabilidad entrecruzadas en la roca, también causan permeabilidades relativas direccionales. Esto nos dice que la permeabilidad relativa debería, en el caso más general, estar representada por un tensor. Las ecuaciones de flujo se convierten entonces en

${\displaystyle \mathbf {u} _{a}=-\mu _{a}^{-1}\mathbf {K} _{ra}\cdot \mathbf {K} \cdot \left(\nabla P_{a}-\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

El caso mencionado anteriormente reflejó la inyección de agua en dirección descendente (o inyección de gas en dirección ascendente) o la producción por agotamiento de presión . Si se inyecta agua en dirección ascendente (o gas en dirección descendente) durante un período de tiempo, se generarán diferentes curvas de permeabilidad relativa en las direcciones x+ y x-. Este no es un proceso de histéresis en el sentido tradicional y no se puede representar con un tensor tradicional. Se puede representar con una declaración IF en el código del software y ocurre en algunos simuladores de yacimientos comerciales. El proceso (o más bien, la secuencia de procesos) puede deberse a un plan de respaldo para la recuperación del campo, o el fluido inyectado puede fluir a otra formación rocosa del yacimiento debido a una parte abierta inesperada de una falla o un cemento que no sella detrás de la carcasa del pozo de inyección. La opción de permeabilidad relativa rara vez se usa y solo observamos que no cambia (la forma analítica de) la ecuación gobernante, sino que aumenta (generalmente duplica) el número de ecuaciones constitutivas para las propiedades involucradas.

La ecuación anterior es una forma vectorial de la ecuación más general para el flujo de fluidos en medios porosos y ofrece al lector una buena descripción general de los términos y las cantidades involucradas. Antes de continuar y transformar la ecuación diferencial en ecuaciones diferenciales para que las utilicen las computadoras, debe escribir la ecuación de flujo en forma de componentes. La ecuación de flujo en forma de componentes (usando la convención de suma ) es

${\displaystyle u{_{a}}^{\sigma }=-\mu _{a}^{-1}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  donde a = w, o, g  donde = 1,2,3 ${\displaystyle \sigma }$

La velocidad de Darcy no es la velocidad de una partícula de fluido, sino el flujo volumétrico (representado frecuentemente por el símbolo ) de la corriente de fluido. La velocidad del fluido en los poros (o velocidad de poro, abreviada pero incorrectamente llamada) está relacionada con la velocidad de Darcy por la relación ${\displaystyle \mathbf {u} _{a}}$ ${\displaystyle \mathbf {q} _{a}}$ ${\displaystyle \mathbf {v} _{a}}$

${\displaystyle \mathbf {v} _{a}=\phi ^{-1}\mathbf {q} _{a}=\phi ^{-1}\mathbf {u} _{a}}$  donde a = w, o, g

El flujo volumétrico es una cantidad intensiva, por lo que no es buena para describir cuánto fluido ingresa por tiempo. La variable preferida para comprender esto es la cantidad extensiva llamada caudal volumétrico, que nos dice cuánto fluido sale (o ingresa) de un área determinada por tiempo, y está relacionada con la velocidad de Darcy por la relación

${\displaystyle Q_{a}=\mathbf {A} \cdot \mathbf {u} _{a}}$  donde a = w, o, g

Observamos que el caudal volumétrico es una cantidad escalar y que la dirección está determinada por el vector normal de la superficie (área) y el flujo volumétrico (velocidad de Darcy). ${\displaystyle Q_{a}}$

En un modelo de yacimiento, el volumen geométrico se divide en celdas de cuadrícula y el área de interés ahora es el área de intersección entre dos celdas adyacentes. Si se trata de celdas vecinas reales, el área es la superficie lateral común y, si una falla divide las dos celdas, el área de intersección suele ser menor que la superficie lateral completa de ambas celdas adyacentes. Por lo tanto, una versión de la ecuación de flujo multifásico, antes de que se discretice y se utilice en simuladores de yacimientos, es

${\displaystyle Q_{a}=-\mu _{a}^{-1}\mathbf {A} \cdot \mathbf {K} _{ra}\cdot \mathbf {K} \cdot \left(\nabla P_{a}-\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

En forma expandida (componente) se convierte en

${\displaystyle Q{_{a}}=-\mu _{a}^{-1}A{_{\sigma }}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  donde a = w, o, g

La presión hidrostática (inicial) a una profundidad (o nivel) z por encima (o por debajo) de una profundidad de referencia z ₀ se calcula mediante

${\displaystyle P{_{a}}=P_{a0}+\int \limits _{z_{0}}^{z}\rho _{a}\left(z\right)g\left(z\right)dz}$  donde a = w, o, g

Cuando se realizan cálculos de presión hidrostática, normalmente no se aplica un subíndice de fase, sino que se cambia la fórmula/cantidad según la fase que se observa a la profundidad real, pero hemos incluido el subíndice de fase aquí para mayor claridad y coherencia. Sin embargo, cuando se ejecutan cálculos de presión hidrostática, se puede utilizar una aceleración de la gravedad que varía con la profundidad para aumentar la precisión. Si no se necesita una precisión tan alta, la aceleración de la gravedad se mantiene constante y la presión calculada se denomina presión de sobrecarga . No se necesita una precisión tan alta en las simulaciones de yacimientos, por lo que la aceleración de la gravedad se trata como una constante en esta discusión. La presión inicial en el modelo de yacimiento se calcula utilizando la fórmula para la presión de sobrecarga (inicial), que es

${\displaystyle P{_{a}}=P_{a0}+g\int \limits _{z_{0}}^{z}\rho _{a}\left(z\right)dz}$  donde a = w, o, g

Para simplificar los términos dentro del paréntesis de la ecuación de flujo, podemos introducir un potencial de flujo llamado potencial -, pronunciado psi-potencial, que se define por ${\displaystyle \psi }$

${\displaystyle \psi _{a}=P{_{a}}-g\int \limits _{z_{0}}^{z}\rho _{a}\left(z\right)dz}$  donde a = w, o, g

Consta de dos términos que son la presión absoluta y la carga gravitatoria. Para ahorrar tiempo de cálculo, la integral se puede calcular inicialmente y almacenar como una tabla para utilizarla en la consulta de tabla, que es más económica desde el punto de vista computacional. La introducción del potencial implica que ${\displaystyle \psi }$

${\displaystyle \nabla \psi _{a}=\nabla P_{a}-\rho _{a}g\nabla z}$  donde a = w, o, g

El potencial psi también se denomina con frecuencia "presión de referencia", ya que la función representa la presión en cualquier punto del yacimiento después de ser transferida al plano de referencia/profundidad z ₀ . En el trabajo práctico de ingeniería es muy útil referir las presiones medidas en los pozos a un nivel de referencia o mapear la distribución de las presiones de referencia en todo el yacimiento. De esta manera, la dirección del movimiento del fluido en el yacimiento se puede ver de un vistazo, ya que la distribución de la presión de referencia es equivalente a la distribución del potencial. Dos ejemplos simples aclararán esto. Un yacimiento puede constar de varias unidades de flujo que están separadas por capas de esquisto apretadas. El fluido de un yacimiento o unidad de flujo puede ingresar a una falla a una profundidad y salir de la falla en otro yacimiento o unidad de flujo a otra profundidad. Del mismo modo, el fluido puede ingresar a un pozo de producción en una unidad de flujo y salir del pozo de producción en otra unidad de flujo o yacimiento.

La ecuación de flujo multifásico para medios porosos ahora se convierte en

${\displaystyle Q_{a}=-\mu _{a}^{-1}\mathbf {A} \cdot \mathbf {K} _{ra}\cdot \mathbf {K} \cdot \nabla \psi _{a}}$  donde a = w, o, g

Esta ecuación de flujo multifásico ha sido tradicionalmente el punto de partida para el programador de software cuando comienza a transformar la ecuación de ecuación diferencial a ecuación de diferencias con el fin de escribir un código de programa para un simulador de yacimientos que se utilizará en la industria petrolera. Las variables dependientes desconocidas han sido tradicionalmente la presión del petróleo (para los yacimientos petrolíferos) y las cantidades volumétricas de los fluidos involucrados, pero se puede reescribir el conjunto total de ecuaciones del modelo que se deben resolver para la presión del petróleo y las cantidades de masa o moles de los componentes del fluido involucrados. ^[7]

Las ecuaciones anteriores están escritas en unidades del SI y asumimos que todas las propiedades de los materiales también están definidas en unidades del SI. Como resultado de esto, las versiones anteriores de las ecuaciones no necesitan ninguna constante de conversión de unidades. La industria petrolera aplica una variedad de unidades, de las cuales al menos dos tienen cierta prevalencia. Si desea aplicar unidades distintas a las unidades del SI, debe establecer constantes de conversión de unidades correctas para las ecuaciones de flujo multifásico.

Conversión de unidades

Las ecuaciones anteriores están escritas en unidades del SI (SI abreviado) suprimiendo que la unidad D (darcy) para la permeabilidad absoluta está definida en unidades no pertenecientes al SI. Es por eso que no hay constantes relacionadas con las unidades. La industria petrolera no utiliza las unidades del SI. En su lugar, utilizan una versión especial de las unidades del SI que llamaremos unidades del SI aplicadas, o utilizan otro conjunto de unidades llamadas unidades de campo que tienen su origen en los EE. UU. y el Reino Unido. La temperatura no está incluida en las ecuaciones, por lo que podemos utilizar el método de etiqueta de factor (también llamado método de factor de unidad) que dice que si tenemos una variable/parámetro con unidad H, multiplicamos esta variable/parámetro por una constante de conversión C y luego la variable obtiene la unidad G que queremos. Esto significa que aplicamos la transformación H*C = G, y el efecto no perteneciente al SI de la definición de permeabilidad se incluye en el factor de conversión C para la permeabilidad. La transformación H*C = G se aplica a cada dimensión espacial, por lo que nos concentramos en los términos principales, ignorando los signos, y luego completamos el paréntesis con el término de gravedad. Antes de comenzar con la conversión, observamos que tanto la ecuación de flujo original (monofásica) de Darcy como las ecuaciones de flujo multifásico generalizadas (o extendidas) de Muskat et al. utilizan la velocidad del yacimiento (flujo volumétrico), la tasa volumétrica y las densidades. Las unidades de estas cantidades reciben el prefijo r (o R) para distinguirlas de sus contrapartes en condiciones de superficie estándar, que reciben el prefijo s (o S). Esto es especialmente importante cuando convertimos las ecuaciones a unidades de campo. La razón por la que entramos en detalles en el tema aparentemente simple de la conversión de unidades es que muchas personas cometen errores al realizar conversiones de unidades.

Ahora estamos listos para comenzar el trabajo de conversión. Primero, tomamos la versión de flujo de la ecuación y la reescribimos como

${\displaystyle 1={\frac {K\nabla _{\gamma }P_{a}}{\mu _{a}u{_{a}}}}}$  donde a = w, o, g

Queremos colocar el factor de conversión compuesto junto con el parámetro de permeabilidad. Aquí notamos que nuestra ecuación está escrita en unidades del SI, y que el grupo de variables/parámetros (en adelante llamados parámetros para abreviar) en el lado derecho constituye un grupo adimensional. Ahora convertimos cada parámetro y reunimos estas conversiones en una única constante de conversión. Ahora notamos que nuestra lista con constantes de conversión (las C) va de la unidad aplicada a las unidades del SI, y esto es muy común para tales listas de conversión. Por lo tanto, asumimos que nuestros parámetros se ingresan en unidades aplicadas y los convertimos (de nuevo) a unidades del SI.

${\displaystyle {\frac {C_{k}KC_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}{C_{\mu }\mu _{a}C_{u}u{_{a}}}}=C_{f}{\frac {K\nabla ^{\gamma }P_{a}}{\mu _{a}u{_{a}}}}=\left({\frac {K\nabla ^{\gamma }P_{a}}{\mu _{a}u{_{a}}}}\right)_{SI}=1}$  donde a = w, o, g

Observe que hemos eliminado la permeabilidad relativa, que es un parámetro adimensional. Este factor de conversión compuesto se denomina constante de Darcy para la ecuación formulada del flujo y es

${\displaystyle C_{f}={\frac {C_{k}C_{\nabla }C_{p}}{C_{\mu }C_{u}}}={\frac {C_{k}C_{\nabla }}{C_{\mu }C_{u}/C_{p}}}}$

Dado que nuestro grupo de parámetros es adimensional en unidades básicas del SI, no necesitamos incluir las unidades del SI en las unidades para nuestro factor de conversión compuesto, como puede ver en la segunda tabla. A continuación, tomamos la versión de tasa de la ecuación y la reescribimos como

${\displaystyle 1={\frac {AK\nabla _{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}}$  donde a = w, o, g

Ahora convertimos cada parámetro y recopilamos estas conversiones en una única constante de conversión.

${\displaystyle {\frac {C_{A}AC_{k}KC_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}{C_{\mu }\mu _{a}C_{Q}Q{_{a}}}}=C_{r}{\frac {AK\nabla ^{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}=\left({\frac {AK\nabla ^{\gamma }P_{a}}{\mu _{a}Q{_{a}}}}\right)_{SI}=1}$  donde a = w, o, g

Observe que hemos eliminado la permeabilidad relativa, que es un parámetro adimensional. Este factor de conversión compuesto se denomina constante de Darcy para la ecuación formulada del flujo y es

${\displaystyle C_{r}={\frac {C_{A}C_{k}C_{\nabla }C_{p}}{C_{\mu }C_{Q}}}={\frac {C_{k}C_{\nabla }C_{A}}{C_{\mu }C_{Q}/C_{p}}}}$

El gradiente de presión y el término de gravedad son idénticos para las ecuaciones de flujo y velocidad y, por lo tanto, se analizarán solo una vez. La tarea aquí es tener un término de gravedad que sea consistente con las unidades aplicadas ("unidades H") para el gradiente de presión. Por lo tanto, debemos colocar nuestro factor de conversión junto con los parámetros de gravedad. Escribimos "el paréntesis" en unidades del SI como

${\displaystyle \nabla ^{\gamma }P_{a}=\rho _{a}g}$  donde a = w, o, g

y reescribirlo como

${\displaystyle 1={\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}}$  donde a = w, o, g

Ahora convertimos cada parámetro y reunimos estas conversiones en una única constante de conversión. En primer lugar, observamos que nuestra ecuación está escrita en unidades del SI y que el grupo de parámetros del lado derecho constituye un grupo adimensional. Por lo tanto, suponemos que nuestros parámetros se ingresan en unidades aplicadas y los convertimos (de nuevo) a unidades del SI.

${\displaystyle {\frac {C_{\rho }\rho _{a}C_{g}g}{C_{\nabla }\nabla ^{\gamma }C_{p}P_{a}}}=C_{pdg}{\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}=\left({\frac {\rho _{a}g}{\nabla ^{\gamma }P_{a}}}\right)_{SI}=1}$  donde a = w, o, g

Esto proporciona el factor de conversión compuesto para la conversión de consistencia como

${\displaystyle C_{pdg}={\frac {C_{\rho }C_{g}}{C_{\nabla }C_{p}}}}$

Dado que nuestro grupo de parámetros es adimensional en unidades SI, no necesitamos incluir las unidades SI en las unidades para nuestro factor de conversión compuesto, como puede ver en la segunda tabla.

Esto es todo para las ecuaciones analíticas, pero cuando el programador transforma la ecuación de flujo en una ecuación de diferencias finitas y luego en un algoritmo numérico, está ansioso por minimizar el número de operaciones computacionales. Aquí hay un ejemplo con dos constantes que se pueden reducir a una mediante la fusión

${\displaystyle C_{cg}=C_{pdg}\cdot g}$

Usando unidades industriales, la versión de flujo de la ecuación de flujo en forma vectorial se convierte en

${\displaystyle \mathbf {u} _{a}=-C_{f}\mu _{a}^{-1}\mathbf {K} _{ra}\mathbf {K} \cdot \left(\nabla P_{a}-C_{pdg}\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

y en forma de componente se convierte en

${\displaystyle u{_{a}}^{\sigma }=-C_{f}\mu _{a}^{-1}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-C_{pdg}\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  donde a = w, o, g  donde = 1,2,3 ${\displaystyle \sigma }$

Utilizando unidades industriales, la versión de tasa de la ecuación de flujo en forma vectorial se convierte en

${\displaystyle Q_{a}=-C_{r}\mu _{a}^{-1}\mathbf {A} \cdot \mathbf {K} _{ra}\cdot \mathbf {K} \cdot \left(\nabla P_{a}-C_{pdg}\rho _{a}\mathbf {g} \right)}$  donde a = w, o, g

y en forma de componente se convierte en

${\displaystyle Q{_{a}}=-C_{r}\mu _{a}^{-1}A{_{\sigma }}K_{ra}{^{\sigma }}_{\beta }K{^{\beta }}_{\gamma }\left(\nabla ^{\gamma }P_{a}-C_{pdg}\rho _{a}g{e_{z}}{^{\gamma }}\right)}$  donde a = w, o, g

La conversión de unidades es una actividad bastante rara, incluso para los profesionales técnicos, pero esa es también la razón por la que la gente olvida cómo hacerlo correctamente.

Véase también

• Yacimiento de petróleo
• Ingeniería de yacimientos
• Modelado de yacimientos
• Simulación de yacimientos
• Ecuaciones del petróleo negro
• Conversión de unidades : incluye tablas de factores de conversión
• Números adimensionales en mecánica de fluidos
• Problema de Fermi : se utiliza para enseñar análisis dimensional
• Método de análisis dimensional de Rayleigh
• Similitud (modelo) : una aplicación del análisis dimensional
• Sistema de medición
• Unidades de medida
• Lista de magnitudes adimensionales
• Órdenes de magnitud (números)
• Análisis dimensional
• Normalización (estadística) y momento estandarizado , conceptos análogos en estadística
• Teorema π de Buckingham

Referencias

1. Muskat M. y Meres MW 1936. El flujo de fluidos heterogéneos a través de medios porosos. Artículo publicado en J. Appl. Phys. 1936, 7, págs. 346-363. https://dx.doi.org/10.1063/1.1745403
2. Muskat M. y Wyckoff RD y Botset HG y Meres MW 1937. Flujo de mezclas de gas y líquido a través de arenas. Publicado en Transactions of the AIME 1937, 123, págs. 69-96. El identificador del documento SPE es SPE-937069-G. https://dx.doi.org/10.2118/937069-G
3. Leverett MC 1941. Comportamiento capilar en sólidos porosos. Documento presentado en la reunión de AIME en Tulsa en octubre de 1940. Publicado en Transactions of the AIME 1941, 142, pp 159–172. El identificador del documento SPE es SPE-941152-G. https://dx.doi.org/10.2118/941152-G
4. Pickup GE y Sorbie KS 1996. La ampliación del flujo bifásico en medios porosos mediante tensores de permeabilidad de fases. Artículo SPE-28586 presentado por primera vez en la Conferencia y exposición técnica anual de la SPE celebrada en Nueva Orleans, EE. UU., del 25 al 28 de septiembre de 1994. Artículo SPE-28586-PA publicado en SPEJ en diciembre de 1996. https://dx.doi.org/10.2118/28586-PA
5. Kumar ATA y Jerauld GR 1996. Impacts of Scaleup on Fluid Flow from Plug to Grid Block Scale in Reservoir Rock. Documento SPE/DOE 35452 presentado en el 10.º Simposio SPE/DOE de 1996 sobre recuperación mejorada de petróleo celebrado en Tulsa, Oklahoma, EE. UU., del 21 al 24 de abril de 1996. Documento SPE-35452-MS publicado por SPE 1996. https://dx.doi.org/10.2118/35452-MS
6. Kyte JR y Berry DW 1975. Nuevas pseudofunciones para controlar la dispersión numérica. Artículo SPE 5105 presentado por primera vez en la 49.ª reunión anual de otoño de la SPE-AIME celebrada en Houston del 6 al 9 de octubre de 1974. Artículo SPE-5105-PA publicado por la SPEJ en agosto de 1975, págs. 269-275. https://dx.doi.org/10.2118/5105-PA
7. Young LC y Stephenson RE 1983. Un enfoque composicional generalizado para la simulación de yacimientos. Documento SPE 10516 presentado por primera vez en el Simposio de Yacimientos de la SPE celebrado en Nueva Orleans, del 31 de enero al 3 de febrero de 1982. Documento SPE-10516-PA publicado por primera vez por SPEJ en octubre de 1983, vol. 23, n.º 05, pp. 727-742; Transacción AIME 275; ahora por ResearchGate; https://www.researchgate.net/publication/244956766 y por https://dx.doi.org/10.2118/10516-PA