Applications¶

[hpl 1: Remaining: comment on boundary conditions for diffusion of substance and heat conduction (include Robin/cooling).]

Diffusion of a substance¶

The first process to be considered is a substance that gets transported through a fluid at rest by pure diffusion. We consider an arbitrary volume $V$ of this fluid, containing the substance with concentration function $c(\x,t)$. Physically, we can think of a very small volume with centroid $\x$ at time $t$ and assign the ratio of the volume of the substance and the total volume to $c(\x, t)$. This means that the mass of the substance in a small volume $\Delta V$ is approximately $\varrho c\Delta V$, where $\varrho$ is the density of the substance. Consequently, the total mass of the substance inside the volume $V$ is the sum of all $\varrho c\Delta V$, which becomes the volume integral $\int_V\varrho cdV$.

Let us reason how the mass of the substance changes and thereby derive a PDE governing the concentration $c$. Suppose the substance flows out of $V$ with a flux $\q$. If $\Delta S$ is a small part of the boundary $\partial V$ of $V$, the volume of the substance flowing out through $dS$ in a small time interval $\Delta t$ is $\varrho \q\cdot\boldsymbol{n} \Delta t \Delta S$, where $\boldsymbol{n}$ is an outward unit normal to the boundary $\partial V$, see Figure. We realize that only the normal component of $\q$ is able to transport mass in and out of $V$. The total outflow of the mass of the substance in a small time interval $\Delta t$ becomes the surface integral

$$ \int\limits_{\partial V} \varrho\q\cdot\boldsymbol{n} \Delta t\, dS\thinspace . $$

Assuming conservation of mass, this outflow of mass must be balanced by a loss of mass inside the volume. The increase of mass inside the volume, during a small time interval $\Delta t$, is

$$ \int\limits_V \varrho (c(\x,t+\Delta t) - c(\x,t)) dV, $$

assuming $\varrho$ is constant, which is reasonable. The outflow of mass balances the loss of mass in $V$, which is the increase with a minus sign. Setting the two contributions equal to each other ensures balance of mass inside $V$. Dividing by $\Delta t$ gives

$$ \int\limits_V \varrho \frac{c(\x,t+\Delta t) - c(\x,t)}{\Delta t} dV = - \int\limits_{\partial V} \varrho\q\cdot\boldsymbol{n}\, dS\thinspace . $$

Note the minus sign on the right-hand side: the left-hand side expresses loss of mass, while the integral on the right-hand side is the gain of mass.

Now, letting $\Delta t\rightarrow 0$, we have

$$ \frac{c(\x,t+\Delta t) - c(\x,t)}{\Delta t} \rightarrow \frac{\partial c}{\partial t}, $$

so

$$ \begin{equation} \int\limits_V \varrho \frac{\partial c}{\partial t}dV + \int\limits_{\partial V} \varrho\q\cdot\boldsymbol{n}\, dS = 0\thinspace . \label{diffu:app:substance:integral} \tag{1} \end{equation} $$

To arrive at a PDE, we express the surface integral as a volume integral using Gauss' divergence theorem:

$$ \int\limits_V (\varrho \frac{\partial c}{\partial t} + \nabla\cdot(\varrho \q))dV = 0\thinspace . $$

Since $\varrho$ is constant, we can divide by this quantity. If the integral is to vanish for an arbitrary volume $V$, the integrand must vanish too, and we get the mass conservation PDE for the substance:

$$ \begin{equation} \frac{\partial c}{\partial t} + \nabla\cdot \q = 0\thinspace . \label{diffu:app:substance:massconv} \tag{2} \end{equation} $$

An arbitrary volume of a fluid.

A fundamental problem is that this is a scalar PDE for four unknowns: $c$ and the three components of $\q$. We therefore need additional equations. Here, Fick's law comes at rescue: it models how the flux $\q$ of the substance is related to the concentration $c$. Diffusion is recognized by mass flowing from regions with high concentration to regions of low concentration. This principle suggests that $\q$ is proportional to the negative gradient of $c$:

$$ \begin{equation} \q = -\dfc\nabla c, \label{diffu:app:substance:Fick} \tag{3} \end{equation} $$

where $\dfc$ is an empirically determined constant. The relation (3) is known as Fick's law. Inserting (3) in (2) gives a scalar PDE for the concentration $c$:

$$ \begin{equation} \frac{\partial c}{\partial t} = \dfc\nabla^2 c\thinspace . \label{diffu:app:substance:PDE} \tag{4} \end{equation} $$

Heat conduction¶

Heat conduction is a well-known diffusion process. The governing PDE is in this case based on the first law of thermodynamics: the increase in energy of a system is equal to the work done on the system, plus the supplied heat. Here, we shall consider media at rest and neglect work done on the system. The principle then reduces to a balance between increase in internal energy and supplied heat flow by conduction.

Let $e(x,t)$ be the internal energy per unit mass. The increase of the internal energy in a small volume $\Delta V$ in a small time interval $\Delta t$ is then

$$ \varrho (e(\x,t+\Delta t) - e(\x,t))\Delta V, $$

where $\varrho$ is the density of the material subject to heat conduction. In an arbitrary volume $V$, as depicted in Figure, the corresponding increase in internal energy becomes the volume integral

$$ \int\limits_V \varrho (e(\x,t+\Delta t) - e(\x,t))dV\thinspace . $$

This increase in internal energy is balanced by heat supplied by conduction. Let $\q$ be the heat flow per time unit. Through the surface $\partial V$ of $V$ the following amount of heat flows out of $V$ during a time interval $\Delta t$:

$$ \int\limits_{\partial V} \q\cdot\boldsymbol{n}\Delta t\, dS\thinspace . $$

The simplified version of the first law of thermodynamics then states that

$$ \int\limits_V \varrho (e(\x,t+\Delta t) - e(\x,t))dV = - \int\limits_{\partial V} \q\cdot\boldsymbol{n}\Delta t\, dS\thinspace . $$

The minus sign on the right-hand side ensures that the integral there models net inflow of heat (since $\boldsymbol{n}$ is an outward unit normal, $\q\cdot\boldsymbol{n}$ models outflow). Dividing by $\Delta t$ and notifying that

$$ \lim_{\Delta t\rightarrow 0} \frac{e(\x,t+\Delta t) - e(\x,t)}{\Delta t} = \frac{\partial e}{\partial t}, $$

we get (in the limit $\Delta t\rightarrow 0$)

$$ \int\limits_V \varrho \frac{\partial e}{\partial t} dV + \int\limits_{\partial V} \q\cdot\boldsymbol{n}\Delta t\, dS = 0\thinspace . $$

This is the integral equation for heat conduction, but we aim at a PDE. The next step is therefore to transform the surface integral to a volume integral via Gauss' divergence theorem. The result is

$$ \int\limits_V\left( \varrho \frac{\partial e}{\partial t} + \nabla\cdot\q\right) dV = 0\thinspace . $$

If this equality is to hold for all volumes $V$, the integrand must vanish, and we have the PDE

$$ \begin{equation} \varrho \frac{\partial e}{\partial t} = -\nabla\cdot\q\thinspace . \label{diffu:app:heat:PDE1} \tag{5} \end{equation} $$

Sometimes the supplied heat can come from the medium itself. This is the case, for instance, when radioactive rock generates heat. Let us add this effect. If $f(\x,t)$ is the supplied heat per unit volume per unit time, the heat supplied in a small volume is $f\Delta t\Delta V$, and inside an arbitrary volume $V$ the supplied generated heat becomes

$$ \int\limits_V f\Delta t dV\thinspace . $$

Adding this to the integral statement of the (simplified) first law of thermodynamics, and continuing the derivation, leads to the PDE

$$ \begin{equation} \varrho \frac{\partial e}{\partial t} = -\nabla\cdot\q + f\thinspace . \label{diffu:app:heat:PDE2} \tag{6} \end{equation} $$

There are four unknown scalar fields: $e$ and $\q$. Moreover, the temperature $T$, which is our primary quantity to compute, does not enter the model yet. We need an additional equation, called the equation of state, relating $e$, $V=1/\varrho=$, and $T$: $e=e(V,T)$. By the chain rule we have

$$ \frac{\partial e}{\partial t} = \left.\frac{\partial e}{\partial T}\right\vert_{V} \frac{\partial T}{\partial t} + \left.\frac{\partial e}{\partial V}\right\vert_{T} \frac{\partial V}{\partial t}\thinspace . $$

The first coefficient $\partial e/\partial T$ is called specific heat capacity at constant volume, denoted by $c_v$:

$$ c_v = \left.\frac{\partial e}{\partial T}\right\vert_{V}\thinspace . $$

The specific heat capacity will in general vary with $T$, but taking it as a constant is a good approximation in many applications.

The term $\partial e/\partial V$ models effects due to compressibility and volume expansion. These effects are often small and can be neglected. We shall do so here. Using $\partial e/\partial t = c_v\partial T/\partial t$ in the PDE gives

$$ \varrho c_v\frac{\partial T}{\partial t} = -\nabla\cdot\q + f\thinspace . $$

We still have four unknown scalar fields ($T$ and $\q$). To close the system, we need a relation between the heat flux $\q$ and the temperature $T$ called Fourier's law:

$$ \q = -k\nabla T, $$

which simply states that heat flows from hot to cold areas, along the path of greatest variation. In a solid medium, $k$ depends on the material of the medium, and in multi-material media one must regard $k$ as spatially dependent. In a fluid, it is common to assume that $k$ is constant. The value of $k$ reflects how easy heat is conducted through the medium, and $k$ is named the coefficient of heat conduction.

We now have one scalar PDE for the unknown temperature field $T(\x,t)$:

$$ \begin{equation} \varrho c_v\frac{\partial T}{\partial t} = \nabla\cdot(k\nabla T) + f\thinspace . \label{diffu:app:heat:PDE} \tag{7} \end{equation} $$

Porous media flow¶

The requirement of mass balance for flow of a single, incompressible fluid through a deformable (elastic) porous medium leads to the equation

$$ S\frac{\partial p}{\partial t} + \nabla\cdot(\q - \alpha\frac{\partial\u}{\partial t}) = 0, $$

where $p$ is the fluid pressure, $\q$ is the fluid velocity, $\u$ is the displacement (deformation) of the medium, $S$ is the storage coefficient of the medium (related to the compressibility of the fluid and the material in the medium), and $\alpha$ is another coefficient. In many circumstances, the last term with $\u$ can be neglected, an assumption that decouples the equation above from a model for the deformation of the medium. The famous Darcy's law relates $\q$ to $p$:

$$ \q = -\frac{K}{\mu}(\nabla p - \varrho\g), $$

where $K$ is the permeability of the medium, $\mu$ is the dynamic viscosity of the fluid, $\varrho$ is the density of the fluid, and $\g$ is the acceleration of gravity, here taken as $\g = -g\kk$. Combining the two equations results in the diffusion model

$$ \begin{equation} S\frac{\partial p}{\partial t} = \mu^{-1}\nabla(K\nabla p) + \frac{\varrho g}{\mu}\frac{\partial K}{\partial z}\thinspace . \label{diffu:app:porous:PDE} \tag{8} \end{equation} $$

Boundary conditions consist of specifying $p$ or $\q\cdot\boldsymbol{n}$ (i.e., normal velocity) at each point of the boundary.

Potential fluid flow¶

Let $\v$ be the velocity of a fluid. The condition $\nabla\times\v =0$ is relevant for many flows, especially in geophysics when viscous effects are negligible. From vector calculus it is known that $\nabla\times\v =0$ implies that $v$ can be derived from a scalar potential field $\phi$: $\v = \nabla\phi$. If the fluid is incompressible, $\nabla\cdot\v = 0$, it follows that $\nabla\cdot\nabla\phi = 0$, or

$$ \begin{equation} \nabla^2\phi = 0\thinspace . \label{_auto1} \tag{9} \end{equation} $$

This Laplace equation is sufficient for determining $\phi$ and thereby describe the fluid motion. This type of flow is known as potential flow. One very important application where potential flow is a good model is water waves. As boundary condition we must prescribe $\v\cdot\boldsymbol{n} =\partial\phi/\partial n$. This gives rise to what is known as a pure Neumann problem and will cause numerical difficulties because $\phi$ and $\phi$ plus any constant are two solutions of the problem. The simplest remedy is to fix the value of $\phi$ at a point.

Streamlines for 2D fluid flow¶

The streamlines in a two-dimensional stationary fluid flow are lines tangential to the flow. The stream function $\psi$ is often introduced in two-dimensional flow such that its contour lines, $\psi = \hbox{const}$, gives the streamlines. The relation between $\psi$ and the velocity field $\v=(u,v)$ is

$$ u = \frac{\partial\psi}{\partial y},\quad v = - \frac{\partial\psi}{\partial x}\thinspace . $$

mathcal{I}t follows that $\nabla\v = \psi{yx} - \psi_{xy}=0$, so the stream function can only be used for incompressible flows. Since

$$ \nabla\times\v = \left(\frac{\partial v}{\partial y} - \frac{\partial u}{\partial x}\right)\kk \equiv \omega\kk, $$

we can derive the relation

$$ \begin{equation} \nabla^2\psi = -\omega, \label{_auto2} \tag{10} \end{equation} $$

which is a governing equation for the stream function $\psi(x,y)$ if the vorticity $\omega$ is known.

The potential of an electric field¶

Under the assumption of time independence, Maxwell's equations for the electric field $\boldsymbol{E}$ become

$$ \begin{align*} \nabla\cdot\boldsymbol{E} &= \frac{\rho}{\epsilon_0},\\ \nabla\times\boldsymbol{E} &= 0, \end{align*} $$

where $\rho$ is the electric charge density and $\epsilon_0$ is the electric permittivity of free space (i.e., vacuum). Since $\nabla\times\boldsymbol{E}=0$, $\boldsymbol{E}$ can be derived from a potential $\varphi$, $\boldsymbol{E} = -\nabla\varphi$. The electric field potential is therefore governed by the Poisson equation

$$ \begin{equation} \nabla^2\varphi = -\frac{\rho}{\epsilon_0}\thinspace . \label{_auto3} \tag{11} \end{equation} $$

If the medium is heterogeneous, $\rho$ will depend on the spatial location $\boldsymbol{r}$. Also, $\epsilon_0$ must be exchanged with an electric permittivity function $\epsilon(\boldsymbol{r})$.

Each point of the boundary must be accompanied by, either a Dirichlet condition $\varphi(\boldsymbol{r}) = \varphi_D(\boldsymbol{r})$, or a Neumann condition $\frac{\partial\varphi(\boldsymbol{r})}{\partial n} = \varphi_N(\boldsymbol{r})$.

[hpl 2: What are appropriate boundary conditions?] [hpl 3: Svein, can you get this generalized to electric field inside some heterogeneous medium? Think we get some variable coefficient equation in that case (anisotropic conduction?). Can you figure it out?] [sl 4: is this what you were thinking of?]

Development of flow between two flat plates¶

Diffusion equations may also arise as simplified versions of other mathematical models, especially in fluid flow. Consider a fluid flowing between two flat, parallel plates. The velocity is uni-directional, say along the $z$ axis, and depends only on the distance $x$ from the plates; $\u = u(x,t)\kk$. The flow is governed by the Navier-Stokes equations,

$$ \begin{align*} \varrho\frac{\partial\u}{\partial t} + \varrho\u\cdot\nabla\u &= -\nabla p + \mu\nabla^2\u + \varrho\f,\\ \nabla\cdot\u &= 0, \end{align*} $$

where $p$ is the pressure field, unknown along with the velocity $\u$, $\varrho$ is the fluid density, $\mu$ the dynamic viscosity, and $\f$ is some external body force. The geometric restrictions of flow between two flat plates puts restrictions on the velocity, $\u = u(x,t)\ii$, and the $z$ component of the Navier-Stokes equations collapses to a diffusion equation:

$$ \varrho\frac{\partial u}{\partial t} = - \frac{\partial p}{\partial z} + \mu\frac{\partial^2 u}{\partial z^2} + \varrho f_z, $$

if $f_z$ is the component of $\f$ in the $z$ direction.

The boundary conditions are derived from the fact that the fluid sticks to the plates, which means $\u=0$ at the plates. Say the location of the plates are $z=0$ and $z=L$. We then have

$$ u(0,t)=u(L,t)=0\thinspace . $$

One can easily show that $\partial p/\partial z$ must be a constant or just a function of time $t$. We set $\partial p/\partial z = -\beta(t)$. The body force could be a component of gravity, if desired, set as $f_z = \gamma g$. Switching from $z$ to $x$ as independent variable gives a very standard one-dimensional diffusion equation:

$$ \varrho\frac{\partial u}{\partial t} = \mu\frac{\partial^2 u}{\partial x^2} + \beta(t) + \varrho\gamma g,\quad x\in [0,L],\ t\in (0,T]\thinspace . $$

The boundary conditions are

$$ u(0,t)=u(L,t)=0, $$

while some initial condition

$$ u(x,0) = I(x) $$

must also be prescribed.

The flow is driven by either the pressure gradient $\beta$ or gravity, or a combination of both. One may also consider one moving plate that drives the fluid. If the plate at $x=L$ moves with velocity $U_L(t)$, we have the adjusted boundary condition

$$ u(L,t) = U_L(t)\thinspace . $$

[hpl 5: Exercises based on this diffusion model.]

Flow in a straight tube¶

Now we consider viscous fluid flow in a straight tube with radius $R$ and rigid walls. The governing equations are the Navier-Stokes equations, but as in the section Development of flow between two flat plates, it is natural to assume that the velocity is directed along the tube, and that it is axi-symmetric. These assumptions reduced the velocity field to $\u = u(r,x,t)\ii$, if the $x$ axis is directed along the tube. From the equation of continuity, $\nabla\cdot\u = 0$, we see that $u$ must be independent of $x$. Inserting $\u = u(r,t)\ii$ in the Navier-Stokes equations, expressed in axi-symmetric cylindrical coordinates, results in

$$ \begin{equation} \varrho\frac{\partial u}{\partial t} = \mu\frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial u}{\partial r}\right) + \beta(t) + \varrho\gamma g,\quad r\in [0,R],\ t\in (0,T]\thinspace . \label{diffu:app:pipeflow:pde} \tag{12} \end{equation} $$

Here, $\beta(t) = -\partial p/\partial x$ is the pressure gradient along the tube. The associated boundary condition is $u(R,t)=0$.

Tribology: thin film fluid flow¶

Thin fluid films are extremely important inside machinery to reduce friction between gliding surfaces. The mathematical model for the fluid motion takes the form of a diffusion problem and is quickly derived here. We consider two solid surfaces whose distance is described by a gap function $h(x,y)$. The space between these surfaces is filled with a fluid with dynamic viscosity $\mu$. The fluid may move partially because of pressure gradients and partially because the surfaces move. Let $U\ii + V\jj$ be the relative velocity of the two surfaces and $p$ the pressure in the fluid. The mathematical model builds on two principles: 1) conservation of mass, 2) assumption of locally quasi-static flow between flat plates.

The conservation of mass equation reads $\nabla\cdot\u$, where $\u$ is the local fluid velocity. For thin films the detailed variation between the surfaces is not of interest, so $\nabla\cdot\u = 0$ is integrated (average) in the direction perpendicular to the surfaces. This gives rise to the alternative mass conservation equation

$$ \nabla\cdot\q = 0,\quad \q = \int\limits_0^{h(x,y)}\u dz, $$

where $z$ is the coordinate perpendicular to the surfaces, and $\q$ is then the volume flux in the fluid gap.

Locally, we may assume that we have steady flow between two flat surfaces, with a pressure gradient and where the lower surface is at rest and the upper moves with velocity $U\ii + V\jj$. The corresponding mathematical problem is actually the limit problem in the section Development of flow between two flat plates as $t\rightarrow\infty$. The limit problem can be solved analytically, and the local volume flux becomes

$$ \q(x,y,z) = \int\limits_0^{h}\u(x,y,z) dz = -\frac{h^3}{12\mu}\nabla p + \frac{1}{2} Uh\ii + \frac{1}{2} Vh\jj\thinspace . $$

The idea is to use this expression locally also when the surfaces are not flat, but slowly varying, and if $U$, $V$, or $p$ varies in time, provided the time variation is sufficiently slow. This is a common quasi-static approximation, much used in mathematical modeling.

Inserting the expression for $\q$ via $p$, $U$, and $V$ in the equation $\nabla\q = 0$ gives a diffusion PDE for $p$:

$$ \begin{equation} \nabla\cdot\left(\frac{h^3}{12\mu}\nabla p\right) = \frac{1}{2}\frac{\partial}{\partial x}(hU) + \frac{1}{2}\frac{\partial}{\partial x}(hV)\thinspace . \label{_auto4} \tag{13} \end{equation} $$

The boundary conditions must involve $p$ or $\q$ at the boundary.

[hpl 6: Exercise!]

Propagation of electrical signals in the brain¶

One can make a model of how electrical signals are propagated along the neuronal fibers that receive synaptic inputs in the brain. The signal propagation is one-dimensional and can, in the simplest cases, be governed by the Cable equation:

$$ \begin{equation} c_m \frac{\partial V}{\partial t} = \frac{1}{r_l}\frac{\partial^2 V}{\partial x^2} - \frac{1}{r_m}V \label{} \end{equation} $$

where $V(x,t)$ is the voltage to be determined, $c_m$ is capacitance of the neuronal fiber, while $r_l$ and $r_m$ are measures of the resistance. The boundary conditions are often taken as $V=0$ at a short circuit or open end, $\partial V/\partial x=0$ at a sealed end, or $\partial V/\partial x \propto V$ where there is an injection of current.