Hello everyone,
I am trying to simulate the Richards equations for unsaturated flow in soil. The equation can be written as
where theta describes soil water retention.
I can’t add more than one embedded image, but we also define the hydraulic conductivity as a function of pressure head (h) :
Se^\tau(1-(1-Se^{1/m})^m)^2
where Se is given by:
Se=\frac{\theta-\theta_r}{\theta_s-\theta_r}
The soil water retention theta and the pressure head h are related in the following way:
theta(h)=\theta_r+\frac{\theta_s-\theta_r}{(1 + |\alpha h|^n)^m}
when h < 0. All other parameters theta_s, theta_r, alpha, n, m, tau, and K_s$ are constants.
I tried to implement a Picard iteration, per Unstable solutions when solving Richards’ equation., but the solver converges a lot faster than expected, and I’m pretty sure the solution it finds is incorrect. Additionally, it will only run for sufficiently large dt (around 50-60s minimum), which seems very strange. For any smaller dt, the Newton solver doesn’t converge.
I was feeling pretty confident in my variational form, but I am fairly new to fenics, so I am wondering if I the set-up is incorrect, or if Picard is not implemented correctly in the time loop. My current code is below:
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import math
# Define Functions
def h(theta):
Se = (theta - theta_r)/(theta_s - theta_r)
return -(((Se**(-1/m))-1)**(1/n))/alpha
def K(theta):
Se = (theta - theta_r)/(theta_s - theta_r)
return K_s*(Se**tau)*(1-(1-(Se**(1/m)))**m)**2
# Define Parameters
theta_r = 0.102
theta_s = 0.5
n = 2
alpha = 0.0335
m = 1 - 1/n
K_s = 0.00922
tau = 0.4
# Time and Spatial Variables
nx = 100 # number of intervals
min_x = -1 # min x value
max_x = 0 # max x value
T = 720
dt = 60
num_steps = math.floor(T/dt)
# Define Boundary Conditions
theta_top = Expression('0.2', degree=0) # Must be between theta_r and theta_s
theta_bottom = Expression('0.4', degree=0)
def boundary_bottom(x):
return bool(x[0]-min_x < DOLFIN_EPS)
def boundary_top(x):
return bool(-x[0] < DOLFIN_EPS)
# Create a Mesh and Function Space
mesh = IntervalMesh(nx, min_x, max_x)
V = FunctionSpace(mesh, 'P', 1)
bc1 = DirichletBC(V, theta_top, boundary_top)
bc2 = DirichletBC(V, theta_bottom, boundary_bottom)
bc = [bc1, bc2]
# Define Initial Value
theta_0 = Expression('-0.2*x[0] + 0.2', degree=1) #Linear profile
theta_n = interpolate(theta_0, V)
theta_k = interpolate(theta_0, V)
# Define Variational Problem
theta = Function(V)
v = TestFunction(V)
z_ex=Expression('x[0]',degree=1)
z=project(z_ex,V)
G = ((theta-theta_n)/dt)*v*dx + (K(theta_k))*dot(grad(h(theta)+z), grad(v))*dx
# Picard Set-Up
tol = 1.0E-5
max_iter = 25
# Time Loop
temp_theta_90 = [0]*101
time = [0]*101
t = 0
count = 0
for n in range(num_steps):
eps = 10000.0
iter = 0
temp_theta = theta_n.vector().get_local()
if count < 101:
temp_theta_90[count] = temp_theta[90]
time[count] = t
while eps > tol and iter < max_iter:
iter += 1
solve(G == 0, theta, bc)
diff = theta.vector().get_local() - theta_k.vector().get_local()
eps = np.linalg.norm(diff, ord=np.Inf)
theta_k.assign(theta)
theta_n.assign(theta)
t += dt
count += 1
plt.plot(time, temp_theta_90)
I would greatly appreciate any insights you might have. Thank you for your time and assistance!