richards.f90 File Reference

This routine solves the Richards equation, then compute the new velocities deducted from the gradients of hydraulic head and from the permeability. These velocities are used for post-processing, calculation of dispersion coefficients, convergence criterion of Newton scheme... but not for transport. In order to ensure the exact conservation of mass, the mass fluxes are computed following the procedure of the standard subroutine resopv (See the resopv section of the theory guide for more informations). More...

Functions/Subroutines
subroutine	richards (icvrge, dt)

Detailed Description

This routine solves the Richards equation, then compute the new velocities deducted from the gradients of hydraulic head and from the permeability. These velocities are used for post-processing, calculation of dispersion coefficients, convergence criterion of Newton scheme... but not for transport. In order to ensure the exact conservation of mass, the mass fluxes are computed following the procedure of the standard subroutine resopv (See the resopv section of the theory guide for more informations).

The hydraulic head is contained in the pressure array, while the velocities are contained in the velocity arrays.

The Richards equation for underground flows is : d theta(h) / dt - div( K(h) grad(h) ) = 0 (or source(h) in very particular cases). This equation can also be written, denoting C(h) = d theta(h)/dh : C(h) dh/dt - div( K(h) grad(h) ) = C(h) dh/dt - d theta(h) / dt (1) The right hand term is close to zero and is not considered in ESTEL. We consider it here for exact mass conservation.

(1) is solved with a 'Newton scheme' handled by tridim. The structure used for this Newton scheme is the same as the one used in the loop over nterup for standard flows. For going from time step n to time step n+1, we use sub-iterations indexed by m : C(h^n) (h^(n+1,m+1)-h^n)/detla t - div( K(h^(n+1,m)) grad(h^(n+1,m+1)) ) = C(h^n) (h^(n+1,m)-h^n)/detla t - (theta(h^(n+1,m))-theta(h(n)))/delta t. These sub-iterations, if they converge, converge to the solution i of the problem.

The Darcy velocity q is then computed thanks to the relation : q = -K(h) grad(h).

This routine is essentially inspired from navstv, resopv and codits.

Please refer to the groundwater flows section of the theory guide for more theoretical informations.

Function/Subroutine Documentation

richards()

subroutine richards	(	integer	icvrge,
		double precision, dimension(:), pointer	dt
	)

Parameters

[in]	icvrge	indicator of convergence
[in]	dt	time step (per cell)