Water flow through unsaturated soils: variational formulation of the Richards - Buckingham equation

Constitutive relationships

Water retention curve, van Genuchten (1980):

\[\begin{aligned} \theta(\psi) = \theta_{r} + (\theta_{s} - \theta_{r}) (1+(a_{VG} \psi)^{n_{VG}})^{-m_{VG}} \\\end{aligned}\]

Hydraulic conductivity function, van Genuchten et al. (1991):

\[\begin{aligned} k(\psi) = k_{sat} \frac {(1-((a_{VG} \psi)^{n_{VG}m_{VG}}) (1+(a_{VG} \psi)^{n_{VG}})^{-m_{VG}}))^2} { (1+(a_{VG} \psi)^{n_{VG}})^{m_{VG}l_{VG}}}\end{aligned}\]

Partial differential equation

Mass conservation law.

\[\begin{aligned} \frac{\partial \theta}{\partial t} = - \nabla \cdot \vec{q} \\\end{aligned}\]

Darcy law.

\[\begin{aligned} \vec{q} = k \left( h, P \right) \nabla H \\\end{aligned}\]

Where h is the pore pressure, P is the position in 3D space and H is the total head.

\[\begin{aligned} P = (x, y, z)\\ h= - \psi \\ \nabla H = \nabla h + \nabla z \\\end{aligned}\]

The partial differential equation can be written as:

\[\begin{aligned} \frac{\partial \theta}{\partial t} = -\nabla \cdot k \left( h, P \right) \nabla H + A, \\\end{aligned}\]

where \(A\) is a source/sink term.

Weak formulation

Integrate on both sides and multiply by the test function \(v\).

\[\begin{aligned} \int_{\Omega}\frac{\partial \theta}{\partial t} v\,d\Omega = -\int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega + \int_{\Omega}A v\,d\Omega \\\end{aligned}\]

Integrate by parts.

\[\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \\ \int_{\Omega}\nabla \cdot \left( k \left( h, P \right) \nabla H v\right) \,d\Omega \\\end{aligned}\]

Apply Gauss divergence theorem.

\[\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N} \cup \Gamma_{D}} k \left( h, P \right) \nabla H v n\,d\Gamma \\ \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N}} k \left( h, P \right) \nabla H v n\,d\Gamma + \int_{\Gamma_{D}} k \left( h, P \right) \nabla H v n\,d\Gamma\end{aligned}\]

Because \(v = 0\) on \(\Gamma_{D}\).

\[\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N}} k \left( h, P \right) \nabla H v n\,d\Gamma \\\end{aligned}\]

Put in the initial integral.

\[\begin{aligned} \int_{\Omega} \frac{\partial \theta}{\partial t} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\\end{aligned}\]

\(\theta \left( h \right)\) is a constitutive relationship. The weak formulation can thus be written on a \(h\) basis.

\[\begin{aligned} \int_{\Omega} \frac{\partial \theta}{\partial t} \times \frac{\partial h}{\partial h} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\ \\ \frac{\partial \theta}{\partial h} \int_{\Omega} \frac{\partial h}{\partial t} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\\end{aligned}\]

Notes

Integration by parts.

\[\begin{aligned} \int_{\Omega} \left(\nabla \cdot u \right) v\,d\Omega = -\int_{\Omega} u \nabla v\,d\Omega + \int_{\Omega} \nabla \cdot \left( u v \right) n\,d\Omega \\\end{aligned}\]

Gauss divergence theorem.

\[\begin{aligned} \int_{\Omega} \nabla \cdot F\, d\Omega = \int_{\partial \Omega = \Gamma_{N} \cup \Gamma_{D}} F \cdot n\, d\Gamma\end{aligned}\]

The problem in steady-state with no Neumann boundary conditions can be defined in Sfepy.

[Someone else is editing this]

You are editing this file