Hello,
I’m trying to solve a different problem: I have a poromechanic Material and I’m pulling on the right side.
I’m not sure, if my implementation of the initial Condition is correct, because I define, that the pressure at
p0 = 10000
but the pressure of my next time step is in a range from 0 (because of my boundary Condition) to -100. I dont understand the big difference between t0 and t1. Is there a mistake in my implementation? Or do I have to implement an initial pressure differently?
Here is the section with the initial condition of my code (and at the bottom is a running example)
P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2) # displacement
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1) # pressure
# Make a mixed space
TH = P2 * P1
V = FunctionSpace(Netz, TH)
# Class representing the intial conditions
class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 0.0 # ux
values[1] = 0.0 # uy
values[2] = 10000.0 # p
def value_shape(self):
return (3,)
# Define functions
(u, p) = TrialFunctions(V) # new Timestep
(v, r) = TestFunctions(V)
du = Function(V) # old Timestep
u0, p0 = split(du)
# Create intial conditions and interpolate
u_init = InitialConditions()
du.interpolate(u_init)
'''2D Poromechanik'''
from fenics import *
from numpy import cos , pi
from ufl import nabla_div
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
# Elasticity parameters
E= 220*1000 # E in Pa
nu = 0.45 # Querkontraktionszahl
b = 1 # Biot Koeffizient b = 1 für inkompr. Fluid + Solid (Castelletto)
k = 2.88*10**(-14) # Permeabilität in m^4/Ns
K_f = 3.3*10**9 #von Wasser in Pa Detournay S- 24
phi = 0.54 # Porosität
M = Constant(K_f/phi)
p_ini = 10000
# Timevariables
T = 10 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
# Define Constants
b = Constant(b) # Biot Koeff
k = Constant(k) # Permeabilität
nu = Constant(nu) # Querkontraktionszahl
G = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))
K = Constant((1+nu)/(3*(1-nu))*2*G*(1-nu)/(1-2*nu))
B_big = Constant(b*M/(K+pow(b,2)*M))
nu_u = Constant((3*nu+b*B_big*(1-2*nu))/(3-b*B_big*(1-2*nu)))
##############################################################################
# Create mesh Beginning #
##############################################################################
# Create mesh and define function space
nx0 = 0
ny0 = 0
nx1 = 0.090 # in m
ny1 = 0.030 # in m
#nx1 = 9*10^-3
#ny1 = 3*10^-3
def mesh(lx,ly, Nx,Ny):
m = UnitSquareMesh(Nx, Ny)
x = m.coordinates()
#Refine on the left and right sides
x[:,0] = (x[:,0] - 0.5) * 2.
x[:,0] = 0.5 * (cos(pi * (x[:,0] - 1.) / 2.) + 1.)
#Scale
x[:,0] = x[:,0]*lx
x[:,1] = x[:,1]*ly
return m
Netz = mesh(nx1,ny1, 40,40)
nCells = Netz.num_cells()
print('Das Netz besteht aus', nCells,'Elementen')
##############################################################################
# Create mesh Ende #
##############################################################################
P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2) # displacement
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1) # pressure
# Make a mixed space
TH = P2 * P1
V = FunctionSpace(Netz, TH)
##############################################################################
# Boundary Stuff: Anfang #
##############################################################################
# Mark boundary subdomians
def boundary(x, on_boundary): # fuer Druck-RB
return on_boundary
# Mark boundary subdomians
class rechts(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],nx1,tol)
class links(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],nx0,tol)
class oben(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],ny1,tol)
class unten(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],ny0,tol)
right = rechts()
left = links()
top = oben()
bottom = unten()
boundaries = MeshFunction('size_t',Netz,Netz.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries,1)
right.mark(boundaries,2)
top.mark(boundaries,3)
bottom.mark(boundaries,4)
ds_test = ds(domain=Netz, subdomain_data=boundaries)
n = FacetNormal(Netz)
# Boundary Conditions
bcp_O = DirichletBC(V.sub(1), Constant(0), top)
bcp_U = DirichletBC(V.sub(1), Constant(0), bottom)
bcu_L = DirichletBC(V.sub(0), Constant((0, 0)), left)
# pulling on the right side
Dehnrate = 1*10**(-3)
u_RB_lin = Expression(("Dehnrate*t"), degree=2, Dehnrate=Dehnrate, t=0)
RB = DirichletBC(V.sub(0).sub(0), u_RB_lin, right)
bc_ges = [bcu_L, bcp_U, bcp_O, RB]
# Class representing the intial conditions
class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 0.0 # ux
values[1] = 0.0 # uy
values[2] = 10000.0 # p
def value_shape(self):
return (3,)
##############################################################################
# Boundary Stuff: Ende #
##############################################################################
# Define functions
Tensor_grad = 2
(u, p) = TrialFunctions(V) # new Timestep
(v, r) = TestFunctions(V)
du = Function(V) # old Timestep
u0, p0 = split(du)
# Create intial conditions and interpolate
u_init = InitialConditions()
du.interpolate(u_init)
theta = 0.5 # 1 -> Euler Rückwärts, 0.5 -> Crank Nicolson
p_mid = (1-theta)*p0+ theta*p
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
# linearer strain tensor
def epsilon(u):
# return 0.5*(nabla_grad(u) + nabla_grad(u).T)
return sym(nabla_grad(u))
def tr_epsilon(u):
return tr(epsilon(u))
# Stress tensor
def sigma(u,p):
return 2*G* epsilon(u) + 2*G * nu/(1-2*nu) *(nabla_div(u))*I - b*p*I
# Momentum balance + Darcy und Fluidmassbalance
Imp =inner(sigma(u,p),nabla_grad(v))*dx
Fl = b*(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx + 1/M * (p-p0)/dt *r*dx +\
k*inner(nabla_grad((p_mid)),nabla_grad(r))*dx
F = Imp + Fl
l = lhs(F)
r = rhs(F)
t = 0
##############################################################################
# Lösungsschleife #
##############################################################################
for n in range(num_steps):
# Update current time
t += dt
u_RB_lin.t=t
w = Function(V)
if n ==0:
w.interpolate(u_init)
solve(l == r, w, bc_ges)
(u,p) = w.split(True)
du.assign(w)