Rigid Plate Boundary Condition

15

Thank you so much for your detailed explanation!!!

I tried to implement, what you wrote, and used

Lag = inner(l,v)*dL + inner(u,e)*dL

for the implementation of the lagrangian multiplier like it was done here: Stress continuity using mixed-dimensional branch

But when I try to run it, I get the error, and I can’t figure out what I did wrong.

> Traceback (most recent call last):
>   File "mandel_v2.py", line 142, in <module>
>     solve(l == r, w, bc_ges)    
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
>     _solve_varproblem(*args, **kwargs)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 269, in _solve_varproblem
>     form_compiler_parameters=form_compiler_parameters)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/problem.py", line 59, in __init__
>     a = Form(a, form_compiler_parameters=form_compiler_parameters)
>   File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/form.py", line 23, in __init__
>     self.subdomains, = list(sd.values())  # Assuming single domain
> ValueError: too many values to unpack (expected 1)

I tried to create a “minimal” example by removing the inital conditions and the plotting etc.

from dolfin import *
import matplotlib.pyplot as plt

# Elasticity parameters 

E= 100*10**6		# E in Pa
nu = 0.25 		# Poissonnumber
k0= 1*10**(-14)
visk = 1 * 10**(-5)
k = k0/visk		# permeability	in m^4/Ns
b = 1           # Biot Koefficent
M = 10**14      # Biot Modul
F = 2        # Force at the top of the plate
a = 1           # length of the plate
G = Constant(E/(2*(1 + 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), 40,40)

# 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()

# Function spaces
LM = VectorElement('Lagrange', mesh.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)

# 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 = split(du)

# 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)) 

# 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(u0)) * sigma_s(u) -n_por(u0)*(p-p0)*I

# Measures
ds_test = ds(domain=mesh, subdomain_data=marker)
dL = Measure('dx', domain=submesh_top)

# Strss-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(l,v)*dL + inner(u,e)*dL
eqn = Imp + Fl + Lag

# free-slip BC
bcu_B  = DirichletBC(V.sub(0).sub(1), Constant((0)), bottom)     
bcu_L  = DirichletBC(V.sub(0).sub(0), Constant((0)), left)   
bcp_R  = DirichletBC(V.sub(1), Constant((0)), right)          

bc_ges = [bcu_L, bcu_B, bcp_R]

l = lhs(eqn)
r = rhs(eqn)

# Solve variational problem 
# Step in time
t = 0
      
for n in range(num_steps):
    # Update current time
    t += dt
    w = Function(V)

    solve(l == r, w, bc_ges)    
    (u,p,l) = w.split(True)

    assign(du.sub(0), u)