# Solution to 1D consolidation problem python implementation

A solution to the 1D consolidation problem is given by

$$\frac{\partial}{\partial t} p = c_{v} \frac{\partial^{2}}{\partial y^{2}} p$$

where $p$ is the pore water pressure, $c_v$ is the consolidation coefficient, $t$ is the elapsed time and $y$ is the height of the column. Further background is provided by this document.

The solution is reproduced below:

$$p(y,t) = \sum_{n=1}^{n = \infty}\frac{2p_0}{n\pi}(1-\cos(n\pi))\sin(\frac{n \pi y}{2h}\exp(\frac{-n^2\pi^2c_v t}{4h^2}))$$

with initial conditions

$$p(y,0) = p_0$$

and boundary conditions

$$p(0,t) = 0$$ $$p(2h,t) = 0$$

To my knowledge that initial condition is a straight line equivalent to the initial loading applied on the soil layer. At $t = 0$ the graph of $y$ versus $p(y,t)$ should look like this for $p_0 = 5000000$

Here is my python implementation of the above solution.

from numpy import cos, exp,pi, linspace, array,sin
import matplotlib.pyplot as plt

#PRESSURE VERSUS HEIGHT

G = 5000000 # Pa = N/m2
E = 10000000# Pa = N/m2
K = E*G/3/(3*G-E)#Bulk modulus of soil skeleton Pa = N/m2
Kl = 2000000000 #Bulk modulus of liquid Pa = N/m2
n = 0.33 #porosity
gammal = 9810 #unit weight of liquid N/m3
k = 1.157E-17 #intrinsic permeability liquid m2
h = 0.5 #m

alpha = 1/(K+4*G/3) #m*s2/kg
beta = 1/Kl #m*s2/kg

cv = k/(gammal*(alpha+n*beta))#m^2/s

n = 200#number of terms to include in series
timevalues = array([0])
yvalues = linspace(0,1,100)
h = 0.5
p0 = 5000000

for t in timevalues:
pvalues = []
for y in yvalues:
SUM = 0
SUM2 = 0
for j in range(1,n,1):
SUM = SUM+(2*p0/j/pi)*(1-cos(j*pi))*sin(j*pi*y/2/h)*exp(-j**2*pi**2**(cv*t/h**2)/4)
pvalues.append(SUM)
plt.plot(pvalues, yvalues)
plt.grid()
plt.show()


I am using arbitrary values for all the parameters and a special relation between liquid bulk modulus, soil bulk modulus and porosity to determine the consolidation coefficient. These things are not relevant however. What is relevant it the initial loading p0 = 5000000 and the timevalues array assignment of 0. I am creating this array and looping over it because later I wish to use other time values. For now however, using the initial time of $t = 0 s$ I am getting the nonsensical result as shown below:

I am plotting the values for $t = 0$. Why am I getting this weird sine-like function rather than a constant function? I am also quite sure the analytical result is correct because I have tested it with sympy. Please help thanks.

• Do you realize that your initial condition does not satisfy the boundary conditions? Since your equation is a diffusion equation, I don't think that the solution is a straight line, at least not in the stationary case. The solution should smooth out from your maximum value to 0 in the borders. Commented Sep 12, 2017 at 17:09
• I agree with you that strictly speaking the initial condition doesn't satisfy the boundary conditions and that that is a problem. However please check page 7 of the link I provided in the original post. You'll see that they are using the same initial and boundary conditions. I was thinking of it more as a step function that goes to zero at the top and bottom but maintains a constant value of p0 throughout the height Commented Sep 12, 2017 at 20:15
• When $t\rightarrow \infty$ the solution should tend to zero. Commented Sep 12, 2017 at 23:25
• According to the document you linked to the solution is: $$p(y,t) = \sum_{n=1}^{n = \infty}\frac{2p_0}{n\pi}(1-\cos(n\pi))\sin(\frac{n \pi y}{2h})\exp(\frac{-n^2\pi^2c_v t}{4h^2})$$, i.e. the exponent should not be in the argument of the sine function. Commented Sep 13, 2017 at 14:32
• Its not in my script! Commented Sep 16, 2017 at 15:19