Hello cdaversin,
Thank you for your answer.
You were right, I wasn’t paying attention with the names for my trial function and lhs. I decomposed my functions further and now it is running without any errors!!!
Although I still have an error, because my pressure is zero (10^(-12)) for all the time steps, which I don’t really understand.
I used the assign(du.sub(0), u) to updated my u^n to my u^n+1 at the end of the time step.
And I don’t update my p0 because I just wanted it as an inital Condition. I define the function p0 and in the equations with p I replaced p with p-p0, and set p0=0 after the first time step so that it is only an initial condition. (At least that was my Plan)
Here is my updated Code:
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
# Elasticity parameters
E= 100*10**6 # E in Pa
nu = 0.25 # Querkontraktionszahl
k0= 1*10**(-14)
visk = 1 * 10**(-5)
k = k0/visk # Permeabilität in m^4/Ns
b = 1 # Biot Koefficent
M = 10**14 # Biot Modul
F = 2 # Force on the plate in N
a = 1 # length of the plate
G = Constant(E/(2*(1 + nu))) # shear modulus
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)))
T = 1 # final time
num_steps = 1000 # number of time steps
dt = T / num_steps # time step size
# Create mesh and define function space
nx0 = 0
ny0 = 0
nx1 = 1
ny1 = 1
mesh = RectangleMesh(Point(nx0, ny0), Point(nx1, ny1), 20,20)
# Mark boundary subdomians
class rechts(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],1,tol)
class links(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0],0,tol)
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],1,tol)
class unten(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[1],0,tol)
marker = MeshFunction('size_t',mesh,mesh.topology().dim()-1 ,0)
TopBoundary().mark(marker,1)
submesh_top = MeshView.create(marker,1)
right = rechts()
left = links()
bottom = unten()
left.mark(marker,2)
right.mark(marker,3)
bottom.mark(marker,4)
# Function spaces
LM = FiniteElement('Lagrange', submesh_top.ufl_cell(), 1) # Lagrange multiplier
#LM = VectorElement('Lagrange', submesh_top.ufl_cell(), 2) # Lagrange multiplier
P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2) # für die Verschiebung quadratisch
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1) # für den Druck linear
# Create mixed space
#TH = MixedElement([P2, P1, LM])
#V = FunctionSpace(mesh, TH)
V = MixedFunctionSpace(FunctionSpace(mesh, P2),
FunctionSpace(mesh, P1),
FunctionSpace(submesh_top, LM))
# Define functions
(u, p, l) = TrialFunctions(V) # Verschiebung, Druck, Lagrangsch-Multiplyer
(v, r, e) = TestFunctions(V)
du = Function(V) # Hilfsgröße zum splitten
u0, p0, l0 = du.sub(0), du.sub(1), du.sub(2)
# Define Constants
n = FacetNormal(mesh)
k = Constant(k) # Permeabilität
nu = Constant(nu) # Querkontraktionszahl
b = Constant(b)
M = Constant(M)
F = Constant(F)
a = Constant(a)
mu = Constant(E/(2*(1 + nu)))
lam = Constant(E*nu/( (1+nu)*(1-2*nu)))
# Fluidflux
q = 0 # -q_vektor*n= q
# Kinematics
d = u.geometric_dimension()
I = Identity(d) # Identity tensor
# linearer Straintensor
def epsilon(u):
# return 0.5*(nabla_grad(u) + nabla_grad(u).T)
return sym(nabla_grad(u))
# Spur des Spannungstensors
def tr_epsilon(u):
return tr(epsilon(u))
#def n_por(u0): # Porositiy
# return 0.2/(1+tr_epsilon(u0))
n_por = 0.1
# linear solid strain
def sigma_s(u):
return lam*tr_epsilon(u)*Identity(d) + 2*mu*epsilon(u)
# strain from both solid and fluid
def sigma(u,p):
return (1-n_por) * sigma_s(u) -n_por*(p)*I # Anteil: n_por*(-p0)*I wird bei der rhs hinzugefügt
# Measures
ds_test = ds(domain=mesh, subdomain_data=marker)
dL = Measure('dx', domain=submesh_top)
# Stress-Vector
T1 = Expression(("0", "-F"),degree =2, F=F)
#Imp =inner(sigma(u,p),nabla_grad(v))*dx -inner(T1,v)*ds_test(1)
#Fl = - k*inner(nabla_grad((p-p0)),nabla_grad(r))*dx -(tr_epsilon(u) - tr_epsilon(u0))/dt *r *dx
#Lag = inner(nabla_div(u),e)*dL #inner(l,v)*dL + inner(u,e)*dL
#eqn = Imp + Fl + Lag
L_eqn1 = inner(sigma(u,p),nabla_grad(v))*dx # momentum balance
R_eqn1 = inner(T1,v)*ds_test(1) - inner(n_por*p0*I,nabla_grad(v))*dx # momentum balance
L_eqn2 = - k*inner(nabla_grad((p)),nabla_grad(r))*dx - tr_epsilon(u)/dt *r *dx # Fluidmass balance with Darcy
R_eqn2 = - tr_epsilon(u0)/dt *r *dx - k*inner(nabla_grad((p0)),nabla_grad(r))*dx
L_eqn3 = inner(nabla_div(u),e)*dL # Lagrangian multiplier
L_eqn = L_eqn1 + L_eqn2 + L_eqn3
R_eqn = R_eqn1 + R_eqn2
# free-slip BC
bcu_B = DirichletBC(V.sub_space(0).sub(1), Constant((0)), bottom)
bcu_L = DirichletBC(V.sub_space(0).sub(0), Constant((0)), left)
bcp_R = DirichletBC(V.sub_space(1), Constant((0)), right)
bc_ges = [bcu_L, bcu_B, bcp_R]
# Initial Conditions:
e_p0 = Expression('1/(3*a) *B_big*(1+nu_u)*F', degree =2, F=F, a=a, B_big=B_big,nu_u=nu_u)
e_u0 = Expression(('(F*nu_u/(2*G))*x[0]/a' , '-F*(1-nu_u)/(2*G)*(x[1]/a)'), degree =2, F=F, a=a, G=G ,nu_u=nu_u)
p_0 = interpolate(e_p0, V.sub_space(1))
u_0 = interpolate(e_u0, V.sub_space(0))
# Solve variational problem
# Step in time
t = 0
vtkfile_p = File('Ergebnisse/p_Mandel_v10.pvd')
vtkfile_u= File('Ergebnisse/u_Mandel_v10.pvd')
for n in range(num_steps):
if n == 0:
# Assign initial Conditions to the functions
assign(p0, p_0)
assign(u0, u_0)
else:
p0 = 0
# Update current time
t += dt
w = Function(V)
solve(L_eqn == R_eqn, w, bc_ges)
# (u,p,l) = w.split(True)
u,p,l = w.sub(0), w.sub(1), w.sub(2)
u.rename('u','u') # that Paraview doesnt view each new time step as a
p.rename('p','p') # new variable
vtkfile_u << (u, t)
vtkfile_p << (p, t)
assign(du.sub(0), u)
`
