Dear everyone, I am trying to solve the Richards’ equation which describes the variably saturated groundwater flow. The equation can be written as (sorry for the coarse equaion expression but I can only upload one picture)
d(theta) / dt=div [ K ( theta ) *grad (h( theta ) + y) ]
where theta is the unknown value and refers to the soil water content, t is time, K and h are soil hydraulic parameters and can be calculated by theta, y refers to the elevation.
I have tried to solve the equation in a UnitSquareMesh with quadrilateral cell type, together with a Dirichlet boundary condition of constant theta (theta_D=0.4) at y=1 and an initial theta of 0.2 among the space to simulate 1D soil water infiltration. Indeed, I have successfully run my code with reference to the non-linear Poisson equation demo, but the result looks strange as the theta at a certain elevation always fluctuates before it begins to raise, while in theroy the theta should stay constant before raising. Below is a picture presenting the results of my code and the results from Hydrus-1D (a mature software for simulating 1D soil water movement by method of finite difference) with totally same boundary condition and soil parameters.
Following are part of my code:
#define time step and dt
dt=60
T=dt
Tmax=7200
#define functionspace
mesh=UnitSquareMesh.create(20,20,CellType.Type.quadrilateral)
V=FunctionSpace(mesh,'CG',2)
#define function and test function
theta=Function(V)
v=TestFunction(V)
#define boundary conditions
theta_D1=Expression('0.4',degree=0)
theta_D2=Expression('0.2',degree=0)
def boundary1(x):
return 1-x[1]<DOLFIN_EPS
def boundary2(x):
return x[1]<DOLFIN_EPS
bc1=DirichletBC(V,theta_D1,boundary1)
bc2=DirichletBC(V,theta_D2,boundary2)
bc=[bc1,bc2]
#define initial condition
theta_init=Expression('0.2',degree=0)
theta0=project(theta_init,V)
#define variational problem
y_ex=Expression('x[1]',degree=1)
y=project(y_ex,V)
F=(theta-theta0)/dt*v*dx+K*dot(grad(h+y),grad(v))*dx
#K,h are soil parameters and have been defined as functions of theta
#solving and time stepping
output=File('2D_Richards_test/solution.pvd')
while T<=Tmax:
solve(F==0,theta,bc)
output << theta
theta0.assign(theta)
T+=dt
It seems my code cannot get stable solutions but I have no idea where the problem is and how to improve my code. Looking foward to replies.Thanks a lot!