Año 2 · Módulo 15: Dinámica espacial y sistemas continuos

141.1. Cuando la posición importa

La temperatura de una lámina metálica cambia con el tiempo, pero no tiene por qué ser igual en todos sus puntos. Una población de microorganismos puede crecer en una zona del suelo mientras desaparece en otra; una concentración química disminuye cerca de una fuente y aumenta varios centímetros más lejos. En todos estos ejemplos, describir el sistema mediante una sola variable temporal es insuficiente. Necesitamos saber qué ocurre, en cada instante, en lugares diferentes y cómo se comunican esos lugares.

Hasta aquí, buena parte de la dinámica se había formulado mediante variables como x(t)x(t) o vectores x(t)\mathbf{x}(t). Esas representaciones condensan el estado en un número finito de componentes. Si consideramos un fluido dividido en miles de celdas, podríamos asignar una variable a cada una; sin embargo, una partición más fina exigiría todavía más variables. El paso a una descripción continua consiste en representar el estado como una función de la posición (u(x,t)u(x,t) en una dimensión, u(x,y,t)u(x,y,t) en una superficie o u(x,y,z,t)u(x,y,z,t) en un volumen). La posición es ahora una variable independiente adicional.

A esta función la llamamos campo. Un campo escalar asigna un número a cada punto, como la temperatura, la presión o la densidad de una especie. Un campo vectorial asigna un vector, como el campo de velocidad v(r,t)\mathbf{v}\left( \mathbf{r},t \right), que indica la dirección y la rapidez del movimiento local. También existen campos tensoriales, necesarios cuando una magnitud incluye relaciones direccionales más complejas, como el esfuerzo mecánico en un sólido.

Un modelo de campos tiene infinitos grados de libertad en el sentido matemático, porque conocer exactamente u(x,t)u(x,t) en todos los puntos de un intervalo requiere especificar una función completa. En cálculo numérico siempre usamos una representación finita, aunque el objeto que intentamos aproximar siga siendo continuo. Habrá que tener presente esta diferencia entre modelo y discretización cuando aparezcan patrones espaciales que desaparecen al cambiar la resolución.

En una descripción espacial, la observación también tiene una escala, porque un sensor promedia lo que ocurre en un volumen pequeño mientras la variable matemática describe un punto ideal. Si ese volumen es grande respecto del patrón estudiado, las derivadas estimadas reflejarán tanto el fenómeno como la resolución del instrumento.

El primer problema es geométrico. Tenemos que precisar cuál es el dominio Ω\Omega (un intervalo, una superficie, un volumen o una región con obstáculos) y qué ocurre en su frontera ∂Ω\partial\Omega. El segundo, dinámico, consiste en encontrar una ley que relacione la evolución temporal del campo con su distribución espacial. Las ecuaciones en derivadas parciales surgen precisamente de esa relación. Su contenido físico procede de las hipótesis que la justifican.

141.2. Derivadas parciales y gradientes

Consideremos la concentración c(x,t)c(x,t) de una sustancia en un tubo. Si observamos un punto fijo del tubo, su evolución queda descrita por la derivada parcial ∂c/∂t\partial c/\partial t. Para calcularla dejamos constante la posición y variamos ligeramente el tiempo,

∂c∂t(x,t)=lim⁡h→0c(x,t+h)−c(x,t)h.\frac{\partial c}{\partial t}(x,t) = \lim_{h \rightarrow 0}\frac{c(x,t + h) - c(x,t)}{h}.

Si congelamos el tiempo y comparamos posiciones cercanas, obtenemos ∂c/∂x\partial c/\partial x. Ambas derivadas pueden tener signos distintos sin contradicción, porque la concentración puede estar aumentando en un sitio y, en ese mismo instante, ser menor a medida que avanzamos hacia la derecha. Una derivada temporal describe el cambio local a lo largo del tiempo; la espacial caracteriza la forma del perfil en un instante dado.

Por ejemplo, tomemos el campo, deliberadamente sencillo,

u(x,t)=e−t(1+x2),x∈[−1,1],t≥0.u(x,t) = e^{- t}\left( 1 + x^{2} \right),\quad\quad x \in \lbrack - 1,1\rbrack,\quad t \geq 0.

Sus derivadas son

∂u∂t=−e−t(1+x2),∂u∂x=2xe−t,∂2u∂x2=2e−t.\frac{\partial u}{\partial t} = - e^{- t}\left( 1 + x^{2} \right),\quad\quad\frac{\partial u}{\partial x} = 2xe^{- t},\quad\quad\frac{\partial^{2}u}{\partial x^{2}} = 2e^{- t}.

En el origen la pendiente espacial es nula, aunque la curvatura es positiva. La amplitud del campo disminuye con el tiempo en todos los puntos. Esta distinción entre pendiente y curvatura anticipa dos clases de mecanismo, el transporte dirigido por un gradiente y la redistribución que depende de cómo cambia ese gradiente en el espacio.

Las unidades ayudan a interpretar las expresiones. Si cc se mide en mol/m³ y xx en metros, entonces ∂c/∂x\partial c/\partial x tiene unidades mol/m⁴ y la derivada temporal, mol/(m³·s). Si una ecuación iguala términos con dimensiones incompatibles, ningún ajuste estadístico posterior la convertirá en una ley de transporte correcta.

Las derivadas parciales requieren regularidad suficiente. Un perfil con un salto abrupto no tiene derivada clásica en el punto del salto, pero ese hecho no impide modelarlo con una ecuación en derivadas parciales. En aplicaciones se usan soluciones débiles, capaces de describir distribuciones discontinuas y de conservar balances integrales. Conviene aprender primero la interpretación de las derivadas clásicas y después reconocer sus límites, porque confundir una discontinuidad física con un fallo del modelo puede llevar a imponer suavidad artificial donde no corresponde.

El gradiente y la dirección de cambio más rápido

En dos o tres dimensiones, una pendiente aislada ya no basta. Un campo escalar u(r,t)u\left( \mathbf{r},t \right) puede variar de forma diferente según la dirección del desplazamiento. El gradiente reúne sus derivadas espaciales en un vector,

∇u=(∂u∂x,∂u∂y,∂u∂z).\nabla u = \left( \frac{\partial u}{\partial x},\frac{\partial u}{\partial y},\frac{\partial u}{\partial z} \right).

Para un desplazamiento infinitesimal drd\mathbf{r}, el cambio de primer orden del campo, a tiempo fijo, es du≃∇u⋅drdu \simeq \nabla u \cdot d\mathbf{r}. Si recorremos una dirección unitaria e\mathbf{e}, no necesariamente paralela a los ejes de coordenadas, la derivada direccional es Deu=∇u⋅eD_{\mathbf{e}}u = \nabla u \cdot \mathbf{e}.

Por la desigualdad de Cauchy–Schwarz, la derivada direccional máxima vale ∣∇u∣|\nabla u| y se alcanza al avanzar en la dirección del gradiente, siempre que este sea distinto de cero. En un mapa de concentración, el gradiente apunta hacia el incremento más rápido de concentración, pero no indica automáticamente la dirección del flujo material; en una difusión ordinaria, el flujo suele orientarse en sentido contrario. El gradiente tampoco es una velocidad, pues sus unidades son las del campo por unidad de longitud.

Supongamos que u(x,y)=9−x2−4y2u(x,y) = 9 - x^{2} - 4y^{2}. El gradiente es ∇u=(−2x,−8y)\nabla u = ( - 2x, - 8y). En el punto (1,1)(1,1), vale (−2,−8)( - 2, - 8), por lo que el aumento más rápido se dirige hacia el interior de las curvas de nivel elípticas. En el origen el gradiente se anula y hay un máximo, aunque una pendiente nula también podría corresponder a un mínimo o a un punto de silla. La curvatura, descrita localmente por derivadas segundas, discrimina esas posibilidades.

Las curvas de nivel u(x,y)=ku(x,y) = k muestran los lugares donde el campo toma el mismo valor. Allí donde el gradiente no se anula, es perpendicular a la curva de nivel. Por esta propiedad, los mapas de isoconcentración, isotermas o isopresión permiten inferir una dirección local de cambio.

En medios anisótropos, el flujo no tiene por qué ser paralelo al gradiente. Una matriz de transporte D\mathbf{D} puede transformar una dirección espacial en otra. Por eso la forma J=−D∇c\mathbf{J} = - D\nabla c, con DD escalar, expresa una hipótesis de isotropía que no se cumple en todos los medios. La geometría del medio (fibras, estratos, paredes, canales) puede modificar decisivamente el sentido y la intensidad de la transferencia.

141.3. Divergencia, entrada y salida en una región diminuta

Un campo vectorial J(r,t)\mathbf{J}\left( \mathbf{r},t \right) puede representar un flujo de masa por unidad de área y tiempo. La divergencia mide la salida neta local de ese flujo y se escribe

∇⋅J=∂Jx∂x+∂Jy∂y+∂Jz∂z.\nabla \cdot \mathbf{J} = \frac{\partial J_{x}}{\partial x} + \frac{\partial J_{y}}{\partial y} + \frac{\partial J_{z}}{\partial z}.

Una divergencia positiva significa que, en una región muy pequeña alrededor del punto, el flujo saliente supera al entrante. Una divergencia negativa corresponde a convergencia neta hacia esa región. La divergencia se calcula sumando derivadas de las componentes del campo y no mide la magnitud del vector; puede haber una velocidad elevada con divergencia nula.

Imaginemos un segmento entre xx y x+Δxx + \Delta x, de sección transversal constante AA. El flujo de una sustancia que entra por la cara izquierda es AJ(x,t)AJ(x,t), mientras que el que sale por la derecha es AJ(x+Δx,t)AJ(x + \Delta x,t). La diferencia entre ambas cantidades, dividida por el volumen AΔxA\Delta x, se aproxima a −∂J/∂x- \partial J/\partial x. Si el material no aparece ni desaparece por otras vías, esa diferencia debe corresponder al cambio temporal de su concentración. En una dimensión,

∂c∂t=−∂J∂x.\frac{\partial c}{\partial t} = - \frac{\partial J}{\partial x}.

La generalización espacial es ∂tc+∇⋅J=0\partial_{t}c + \nabla \cdot \mathbf{J} = 0. Es una ley de balance; establece que una variación local de concentración debe corresponder a la entrada o salida neta por la frontera, pero todavía no prescribe cómo fluye el material.

La relación entre divergencia local y flujo global se expresa mediante el teorema de Gauss,

∫Ω∇⋅J dV=∫∂ΩJ⋅n dS.\int_{\Omega}^{}\nabla \cdot \mathbf{J}\, dV = \int_{\partial\Omega}^{}\mathbf{J} \cdot \mathbf{n}\, dS.

Aquí n\mathbf{n} es la normal unitaria exterior. El lado derecho mide todo lo que atraviesa la frontera de Ω\Omega, con signo positivo para el flujo saliente. El lado izquierdo suma las divergencias interiores. El teorema permite pasar del balance sobre una región observable a una ecuación local, siempre que el campo y el dominio tengan la regularidad adecuada. En contextos menos regulares hay versiones integrales o débiles que conservan el mismo significado.

Rotación local y divergencia no son lo mismo

La divergencia no agota las propiedades espaciales de un campo vectorial. Para describir la circulación local se introduce el rotacional, el vector ∇×v\nabla \times \mathbf{v}. Una corriente puede girar alrededor de un obstáculo sin acumular fluido en ninguna región; en un modelo bidimensional incompresible podemos encontrar divergencia nula y rotacional distinto de cero. También es posible construir campos con divergencia positiva y rotacional nulo. Son propiedades independientes.

En un campo de velocidad, el rotacional es la vorticidad; la velocidad angular de un elemento material requiere una interpretación específica y, en el régimen usual de continuo, es la mitad de la vorticidad local. Para estudiar la transferencia de sustancias, en cambio, el balance escalar depende directamente de la divergencia del flujo de sustancia. Una circulación intensa puede redistribuir la concentración sin crear ni destruir materia. Esta distinción impide deducir acumulación a partir de la presencia de remolinos o ausencia de movimiento a partir de una divergencia nula.

141.4. Laplaciano y comparación de un punto con su entorno

La combinación ∇⋅(∇u)\nabla \cdot (\nabla u) recibe el nombre de laplaciano. En coordenadas cartesianas,

Δu=∇2u=∂2u∂x2+∂2u∂y2+∂2u∂z2.\Delta u = \nabla^{2}u = \frac{\partial^{2}u}{\partial x^{2}} + \frac{\partial^{2}u}{\partial y^{2}} + \frac{\partial^{2}u}{\partial z^{2}}.

En una dimensión el laplaciano coincide con la segunda derivada espacial. Si el perfil tiene un máximo local estricto y es dos veces diferenciable, la segunda derivada suele ser negativa en la dirección de mayor concavidad; en un mínimo ocurre lo contrario. Para una función suficientemente regular en varias dimensiones, el laplaciano es la suma de las curvaturas alineadas con los ejes cartesianos. Un máximo interior implica Δu≤0\Delta u \leq 0, aunque Δu=0\Delta u = 0 es compatible con curvaturas no nulas en algunas direcciones.

La interpretación del laplaciano como comparación con el vecindario es más intuitiva en una malla. En una dimensión, con separación hh,

u(x+h)−2u(x)+u(x−h)h2→∂2u∂x2\frac{u(x + h) - 2u(x) + u(x - h)}{h^{2}} \rightarrow \frac{\partial^{2}u}{\partial x^{2}}

cuando hh tiende a cero y la función es suficientemente suave. Si el valor del centro supera al promedio de los vecinos, la expresión es negativa. Una ecuación ∂tu=DΔu\partial_{t}u = D\Delta u, con D>0D > 0, hace disminuir ese máximo local. Si el centro está por debajo de sus vecinos, el término es positivo y tiende a elevarlo. Así aparece la tendencia a la suavización que caracteriza la difusión.

En dos dimensiones, la aproximación habitual incluye cuatro vecinos. Al pasar a una malla tridimensional se añaden los pares correspondientes al eje vertical. La forma concreta depende de la geometría y en coordenadas cilíndricas o esféricas incluye factores métricos adicionales. Copiar sin modificaciones ∂xx+∂yy\partial_{xx} + \partial_{yy} a un problema de simetría radial equivale a cambiar la ecuación.

La ecuación Δu=0\Delta u = 0, llamada ecuación de Laplace, describe, bajo hipótesis particulares, estados estacionarios de muchos modelos sin fuentes internas. Sus soluciones son funciones armónicas. Dentro de una región, el valor de una función armónica tiene una propiedad de promedio y, salvo que la función sea constante, no tiene máximos interiores. Estas propiedades reaparecen en problemas de potencial, electrostática y equilibrio térmico y anticipan criterios de consistencia para las soluciones numéricas.

Curvas de nivel elípticas con su gradiente normal y una parábola con pendiente nula y curvatura positiva en el centro.
Figura 141.1. A la izquierda, el gradiente de u = 9 − x² − 4y² en P = (1, 1) es (−2, −8) y resulta perpendicular a la tangente de la curva de nivel. A la derecha, u = 1 + x² tiene pendiente cero en x = 0 y segunda derivada positiva. La segunda diferencia con h = 0,5 reproduce exactamente esa curvatura. Son campos matemáticos estáticos y adimensionales. Ampliar figura ↗

141.5. Del balance integral a una ecuación local

Para comprender por qué aparecen ecuaciones diferenciales con términos tan distintos, conviene empezar por un volumen fijo Ω\Omega y una cantidad por unidad de volumen u(r,t)u\left( \mathbf{r},t \right). Su contenido total es M(t)=∫Ωu dVM(t) = \int_{\Omega}^{}u\, dV. Si el material entra y sale mediante un flujo J\mathbf{J} y hay fuentes o sumideros volumétricos ss, el balance toma la forma

ddt∫Ωu dV=−∫∂ΩJ⋅n dS+∫Ωs dV.\frac{d}{dt}\int_{\Omega}^{}u\, dV = - \int_{\partial\Omega}^{}\mathbf{J} \cdot \mathbf{n}\, dS + \int_{\Omega}^{}s\, dV.

El signo negativo delante del flujo superficial indica que una salida neta reduce la cantidad almacenada. En un volumen de control fijo y bajo condiciones que permiten intercambiar derivada e integral, el teorema de Gauss da

∫Ω(∂tu+∇⋅J−s)dV=0.\int_{\Omega}^{}\left( \partial_{t}u + \nabla \cdot \mathbf{J} - s \right)dV = 0.

Si la igualdad vale para toda subregión adecuada y las funciones tienen regularidad suficiente, concluimos localmente que

∂tu+∇⋅J=s.\boxed{\partial_{t}u + \nabla \cdot \mathbf{J} = s.}

La ecuación es una plantilla para muchas teorías de transporte. La cantidad uu podría ser una densidad de partículas, la masa de un contaminante o la energía por unidad de volumen; de un caso a otro cambian las unidades, la expresión del flujo y el significado de las fuentes. Conviene distinguir una fuente volumétrica de una entrada desde la frontera, porque corresponden a términos diferentes y sumarlas indiscriminadamente duplicaría la contribución física.

Supongamos ahora que el material es transportado por una velocidad v\mathbf{v} y, además, se mueve respecto del fluido mediante otro mecanismo. Podemos escribir J=uv+Jrel\mathbf{J} = u\mathbf{v} + \mathbf{J}_{rel}. Si el transporte relativo es difusivo, una aproximación habitual es Jrel=−D∇u\mathbf{J}_{rel} = - D\nabla u. El balance conduce a

∂tu+∇⋅(uv)=∇⋅(D∇u)+s.\partial_{t}u + \nabla \cdot \left( u\mathbf{v} \right) = \nabla \cdot (D\nabla u) + s.

El término de advección responde al movimiento colectivo; el de difusión redistribuye la cantidad según una relación constitutiva entre flujo y gradiente; las fuentes describen producción o desaparición. La ecuación no implica que esos tres procesos estén siempre presentes ni que los coeficientes sean constantes. En particular, un DD dependiente de la posición debe permanecer dentro de la divergencia, salvo que conozcamos condiciones que permitan simplificarlo.

141.6. Seguir un punto o seguir una partícula

La notación ∂tu\partial_{t}u presupone que observamos el campo en una posición fija. En un fluido, una partícula se desplaza y la variación que experimenta puede ser distinta de la que mide un sensor inmóvil. Si su trayectoria es r(t)\mathbf{r}(t) y satisface dr/dt=v(r(t),t)d\mathbf{r}/dt = \mathbf{v}\left( \mathbf{r}(t),t \right), la regla de la cadena da

ddtu(r(t),t)=∂u∂t+v⋅∇u.\frac{d}{dt}u\left( \mathbf{r}(t),t \right) = \frac{\partial u}{\partial t} + \mathbf{v} \cdot \nabla u.

Esta combinación se llama derivada material o sustancial, Du/DtDu/Dt. El primer término mide el cambio observado en un punto fijo; el segundo aparece porque la partícula visita posiciones con valores diferentes del campo. Un río puede mantener una temperatura estacionaria en cada kilómetro y, pese a ello, calentar un objeto que viaja aguas abajo hacia una región más cálida. En ese caso ∂tT=0\partial_{t}T = 0, pero DT/DtDT/Dt puede ser positivo.

La distinción también interviene al reescribir balances. Al expandir ∇⋅(uv)\nabla \cdot \left( u\mathbf{v} \right) obtenemos v⋅∇u+u∇⋅v\mathbf{v} \cdot \nabla u + u\nabla \cdot \mathbf{v}. Por tanto, el balance advectivo tiene la forma

DuDt+u∇⋅v=−∇⋅Jrel+s.\frac{Du}{Dt} + u\nabla \cdot \mathbf{v} = - \nabla \cdot \mathbf{J}_{rel} + s.

El término u∇⋅vu\nabla \cdot \mathbf{v} describe el efecto de la expansión o compresión del medio y desaparece solo si ∇⋅v=0\nabla \cdot \mathbf{v} = 0, como sucede en una descripción de flujo incompresible. Confundir la derivada material con la derivada parcial en un fluido compresible equivale a perder un mecanismo de acumulación de densidad.

Este ejemplo muestra por qué una ecuación debe leerse antes de resolverse. Dos expresiones parecidas pueden corresponder a observadores diferentes, a cantidades distintas o a una hipótesis adicional sobre el medio. La equivalencia algebraica depende de especificar qué se conserva y de mantener todas las contribuciones de transporte.

141.7. Familias de ecuaciones y naturaleza de sus datos

Las ecuaciones en derivadas parciales se clasifican con frecuencia, al menos para problemas lineales de segundo orden, como elípticas, parabólicas o hiperbólicas. La clasificación general se obtiene de los coeficientes del término diferencial de mayor orden, llamado parte principal. Es una guía para reconocer el tipo de comportamiento esperado y las condiciones que suelen ser apropiadas; en ecuaciones no lineales, degeneradas o de tipo mixto, la clasificación requiere precauciones adicionales.

La ecuación de Laplace, Δu=0\Delta u = 0, es el ejemplo elíptico. Un problema típico prescribe valores o flujos en la frontera de una región y busca la distribución espacial de equilibrio. La información de la frontera afecta globalmente a la solución, sin que intervenga necesariamente una evolución temporal. Cuando agregamos una fuente, se obtiene la ecuación de Poisson, −Δu=f- \Delta u = f, con una convención de signos que puede cambiar según la disciplina.

La ecuación del calor, ∂tu=DΔu\partial_{t}u = D\Delta u, es el paradigma parabólico cuando D>0D > 0 es constante. Evoluciona desde una condición inicial y tiende a suavizar las perturbaciones espaciales. En el modelo matemático ideal, una perturbación localizada tiene influencia no nula a distancias arbitrarias para cualquier tiempo positivo, aunque esa influencia pueda ser extremadamente pequeña. Esta propagación instantánea del soporte es una propiedad matemática del modelo difusivo clásico y no implica una transferencia física con velocidad infinita.

La ecuación de ondas, ∂ttu=c2Δu\partial_{tt}u = c^{2}\Delta u, es hiperbólica. En su régimen ideal describe señales que se propagan con velocidad finita cc. Suele necesitar posición inicial y velocidad inicial, porque contiene una segunda derivada temporal. Sus soluciones pueden conservar frentes pronunciados y mostrar reflexión o interferencia. La ecuación de transporte ∂tu+v∂xu=0\partial_{t}u + v\partial_{x}u = 0 es un ejemplo hiperbólico de primer orden; aunque queda fuera del criterio algebraico usado para las ecuaciones lineales de segundo orden, comparte con la ecuación de ondas la propagación a lo largo de características.

La distinción evita elegir procedimientos inadecuados. Resolver una ecuación del calor con dos condiciones iniciales independientes, como si fuera una ecuación de ondas, suele sobredeterminar el problema. Usar un método que difumina intensamente una onda puede crear disipación numérica espuria. A su vez, interpretar un problema elíptico como una secuencia temporal física puede atribuir significado a un tiempo artificial introducido únicamente para hallar una solución estacionaria.

Condiciones iniciales, fronteras y unicidad

Una EDP escrita sin dominio y sin datos complementarios rara vez define un problema completo. Para la ecuación ∂tu=D∂xxu\partial_{t}u = D\partial_{xx}u en 0<x<L0 < x < L, debemos indicar el perfil inicial u(x,0)=u0(x)u(x,0) = u_{0}(x) y especificar lo que ocurre en x=0x = 0 y x=Lx = L. Las condiciones de frontera describen mecanismos físicos que pueden modificar la masa total, la estabilidad y los estados estacionarios.

En una condición de Dirichlet fijamos el valor del campo, por ejemplo u(0,t)=uLu(0,t) = u_{L}. Puede corresponder a una superficie mantenida a temperatura constante por un reservorio suficientemente grande. La condición de Neumann prescribe una derivada normal o, en problemas de transporte, el flujo normal J⋅n\mathbf{J} \cdot \mathbf{n}. Una frontera impermeable impone flujo nulo; si el medio es isótropo y D>0D > 0, equivale a ∂nu=0\partial_{n}u = 0. Una condición de Robin relaciona valor y flujo, como sucede en una transferencia finita hacia un entorno exterior.

También hay fronteras periódicas, apropiadas para dominios que representan una celda repetida o una geometría cerrada sobre sí misma. La periodicidad iguala valores y flujos correspondientes en los bordes identificados y es distinta de una pared impermeable. En dos y tres dimensiones pueden coexistir segmentos con condiciones distintas, de modo que incluso la intersección entre fronteras requiere atención.

La unicidad depende de la ecuación y de los datos. Para Δu=0\Delta u = 0 con datos de Dirichlet bien planteados, la solución es única bajo hipótesis estándar. Si solo prescribimos la derivada normal, debe satisfacerse una condición de compatibilidad sobre el flujo neto y la solución de Laplace, cuando existe, queda determinada salvo una constante aditiva. La interpretación física de esta libertad es que conocer únicamente los flujos no fija un nivel absoluto de potencial. Conviene reconocerla antes de atribuir a un error del programa la singularidad de una matriz numérica.

141.8. Ejemplo completo de la temperatura a lo largo de una barra

Tomemos una barra homogénea y delgada de longitud LL, con aislamiento lateral ideal, conductividad térmica kk, densidad ρ\rho y calor específico cpc_{p}, todos constantes. Si la sección transversal es AA, la energía térmica de una porción [x,x+dx]\lbrack x,x + dx\rbrack cambia cuando entra o sale calor por sus extremos. La ley de Fourier da el flujo de calor por unidad de área, q=−k∂xTq = - k\partial_{x}T. El balance energético conduce a

ρcp∂T∂t=k∂2T∂x2,α=kρcp.\rho c_{p}\frac{\partial T}{\partial t} = k\frac{\partial^{2}T}{\partial x^{2}},\quad\quad\alpha = \frac{k}{\rho c_{p}}.

La difusividad térmica α\alpha tiene unidades m²/s. Supongamos que los extremos se mantienen a una misma temperatura TbT_{b} y que inicialmente la barra tiene una perturbación sinusoidal,

T(x,0)=Tb+A0sin⁡(πxL),T(0,t)=T(L,t)=Tb.T(x,0) = T_{b} + A_{0}\sin\left( \frac{\pi x}{L} \right),\quad T(0,t) = T(L,t) = T_{b}.

Una solución exacta es

T(x,t)=Tb+A0exp⁡(−απ2tL2)sin⁡(πxL).T(x,t) = T_{b} + A_{0}\exp\left( - \frac{\alpha\pi^{2}t}{L^{2}} \right)\sin\left( \frac{\pi x}{L} \right).

Podemos comprobarla directamente. La derivada temporal multiplica la perturbación por −απ2/L2- \alpha\pi^{2}/L^{2}; la segunda derivada espacial produce −π2/L2- \pi^{2}/L^{2} y la ecuación queda satisfecha. El perfil inicial coincide con el prescrito y la función seno se anula en ambos extremos. Sin observar experimentalmente la barra, hemos verificado que la expresión es una solución matemática del problema idealizado.

La amplitud decae exponencialmente con un tiempo característico τ=L2/(απ2)\tau = L^{2}/\left( \alpha\pi^{2} \right). Duplicar la longitud, conservando el material y el tipo de frontera, multiplica por cuatro este tiempo. Este resultado tiene una explicación geométrica, ya que la redistribución local es eficiente a escalas cortas, mientras que las distancias grandes exigen períodos mayores. No debemos extrapolar el cálculo a cualquier barra real, porque las pérdidas laterales, una conductividad dependiente de la temperatura o una resistencia de contacto introducirían términos y condiciones diferentes.

El ejemplo muestra además cómo influye la forma inicial. Una distribución arbitraria puede descomponerse, bajo condiciones apropiadas, en una suma de modos sinusoidales. Cada modo decae con una tasa proporcional al cuadrado de su número de onda, de modo que las oscilaciones espaciales finas se atenúan antes que las amplias. La suavización que intuimos al observar un campo se convierte así en una afirmación cuantitativa acerca de escalas espaciales.

Perfiles sinusoidales de temperatura que decaen y amplitudes de barras de uno y dos metros; sus tiempos de relajación difieren por un factor cuatro.
Figura 141.2. La solución exacta T = 20 + 10 exp(−απ²t/L²) sen(πx/L) cumple las temperaturas fijas de los extremos. Los perfiles de la izquierda corresponden a L = 1 m. A la derecha, el primer modo de la barra de 2 m necesita cuatro veces más tiempo para reducir su amplitud al factor 1/e. α = 10⁻⁴ m²/s es un parámetro didáctico, no una medición de un material. Ampliar figura ↗

141.9. Adimensionalizar para reconocer escalas

Una ecuación puede contener pocos términos y aun así esconder escalas temporales muy diferentes. Consideremos advección, difusión y reacción lineal en una dimensión,

∂tu+v∂xu=D∂xxu−ku,\partial_{t}u + v\partial_{x}u = D\partial_{xx}u - ku,

con velocidad vv, difusividad D>0D > 0 y tasa de desaparición k≥0k \geq 0. Elegimos una longitud característica LL, una amplitud UU y un tiempo L2/DL^{2}/D. Definimos ξ=x/L\xi = x/L, τ=tD/L2\tau = tD/L^{2} y w=u/Uw = u/U. Al sustituir, obtenemos

∂τw+Pe ∂ξw=∂ξξw−Da w,Pe=vLD,Da=kL2D.\partial_{\tau}w + Pe\,\partial_{\xi}w = \partial_{\xi\xi}w - Da\, w,\quad Pe = \frac{vL}{D},\quad Da = \frac{kL^{2}}{D}.

El número de Péclet PePe compara el transporte advectivo con el difusivo usando la escala escogida. El número DaDa, aquí una forma de número de Damköhler, compara la reacción con la difusión. La misma ecuación dimensional puede comportarse de manera muy distinta según estos números. Con PePe grande, un frente puede avanzar principalmente con el flujo; con DaDa grande, la sustancia puede desaparecer antes de atravesar el dominio por difusión.

Estas razones no son constantes universales del material. Cambiar la longitud característica cambia ambos números, aunque vv, DD y kk permanezcan iguales. Un proceso dominado por difusión a escala microscópica puede estar dominado por advección cuando se observa a escala macroscópica. La adimensionalización permite preguntar qué términos importan para la escala y el fenómeno estudiados.

141.10. Soluciones débiles, mallas y límites del modelo continuo

Las expresiones diferenciales anteriores suponen campos suficientemente regulares. Una sustancia puede estar inicialmente concentrada en una región pequeña, una onda puede presentar un salto y un medio puede tener interfaces entre materiales. Para trabajar con estas situaciones, la teoría de EDP usa formulaciones débiles, en las que se multiplica la ecuación por funciones de prueba, se integra sobre el dominio y, mediante integración por partes, se transfieren derivadas entre funciones. Así pueden definirse soluciones cuya regularidad no permite interpretar todas las derivadas punto por punto.

En una dimensión, la ecuación de balance ∂tu+∂xJ=s\partial_{t}u + \partial_{x}J = s integrada contra una función de prueba suave φ\varphi da, formalmente,

∫0L∂tu φ dx−∫0LJ ∂xφ dx+[Jφ]0L=∫0Lsφ dx.\int_{0}^{L}\partial_{t}u\,\varphi\, dx - \int_{0}^{L}J\,\partial_{x}\varphi\, dx + \lbrack J\varphi\rbrack_{0}^{L} = \int_{0}^{L}s\varphi\, dx.

La expresión conserva explícitamente la influencia de la frontera y en ella se inspiran varios métodos computacionales. En un método de volúmenes finitos, por ejemplo, se sigue la cantidad contenida en cada celda y se calculan los flujos que cruzan sus caras. Si el mismo flujo aparece con signos opuestos en dos celdas vecinas, la suma sobre el sistema respeta la conservación global, salvo aportes externos y errores de redondeo.

La discretización obliga a decidir el tamaño de celda, el paso temporal, la representación de las geometrías curvas y la aproximación de los gradientes. Una malla gruesa puede suavizar una estructura real o crear dependencias artificiales entre regiones. Una malla más fina no siempre soluciona el problema, porque el método temporal puede hacerse inestable y los parámetros físicos pueden haber sido estimados para una escala diferente. Es necesario comprobar la convergencia por refinamiento, los balances conservados y la sensibilidad a las condiciones de frontera. También conviene comparar con una solución exacta sencilla cuando exista.

El continuo es, finalmente, una idealización. Una densidad de población deja de ser una variable suave cuando hay muy pocos individuos por celda; una concentración molecular fluctúa cuando el número de partículas es pequeño. Las ecuaciones de campos pueden surgir como límites aproximados de modelos discretos o estocásticos, pero el tamaño del dominio, las escalas de correlación y las fluctuaciones intervienen al decidir cuándo ese límite es válido. El paso desde las trayectorias aleatorias hacia una ecuación de campos puede revelar regularidades colectivas, a costa de perder información sobre las realizaciones individuales.

Las interacciones que no caben en una derivada local

Una EDP formulada con gradientes y derivadas en un punto presupone alguna forma de interacción local, en la que el cambio en una región depende de lo que ocurre allí y en posiciones arbitrariamente próximas. La hipótesis vale para muchos procesos de transporte, pero una población espacial también puede depender de conexiones a larga distancia. Una semilla desplazada por un animal, una persona que viaja entre dos ciudades o una molécula transportada por una corriente rápida pueden acoplar lugares que no son vecinos euclidianos.

Cuando estos saltos son parte esencial del fenómeno, la ecuación local puede ampliarse con términos integrales, operadores no locales o una estructura de red. La decisión requiere definir qué magnitud se conserva durante el salto y cómo se asigna la salida de un lugar a la entrada de otro. Si un modelo describe el transporte entre regiones con un núcleo integral, ese núcleo debe expresar las distancias, probabilidades o intensidades de conexión pertinentes.

Esta observación ayuda a interpretar las fronteras. El borde geométrico de un mapa puede ser solo el límite de la región que decidimos observar, mientras que una membrana impermeable es una condición física sobre el flujo. Un modelo puede usar el mismo operador diferencial en el interior y dar resultados radicalmente distintos según las transferencias que permita su frontera. La geometría forma parte del problema, aunque las conexiones del sistema no se reducen a ella.

La continuidad conceptual entre ambos niveles es especialmente clara en la difusión. A partir de una regla local de transporte y de un balance de materia se obtiene una ecuación para la concentración espacial. Por otra vía, el límite de muchos desplazamientos aleatorios independientes produce, bajo hipótesis apropiadas, la misma estructura matemática. El capítulo siguiente construye esa conexión y muestra qué propiedades de la ecuación proceden de la conservación, cuáles dependen de la constitución del medio y cuáles pertenecen a la geometría del dominio.

Referencias