Hi,
I am pretty new to FEniCS and I have a problem creating a neumann boundary condition on a curved boundary.
So that you understand my problem, I am going to explain it in more detail:
I have a coupled field problem - a solid mechanical field (with the displacement u) and a fluid mechanical field (with the fluid pressure p). The theory is known as the Theory of porous media.
Now I know the force which points in y-direction and which should be applied to the upper boundary (that is in reality a 3D-geometry, of course). That can be converted into a stress vector. If the stress vector is applied on the upper boundary, I expect that the geometry deforms in a way, that the upper bundary shifts downwards and the right boundary shifts sideways and a little bit downwards (like a gas pedal). The lower boundary is fixed.
The results I get are strange. The area near the right boundary deforms. However, the upper boundary doesn’t deform anyway. I don’t know exactly, whether the implementation of the traction or stress vector is false, because it is a curved boundary. Can anyone help me with this?
This is my code:
''' coupled field problem '''
from fenics import *
import matplotlib.pyplot as plt
import sys
import numpy as np
from numpy import array
from ufl import nabla_div
from numpy import cos , pi
import time
import pylab as lab
start = time.time()
Titel = 'Mechanical'
##############################################################################
# Definieren von Parametern #
##############################################################################
## Materialeigenschaften
#Lame-Konstanten
mu = 0.05 # Paper Travascio
lambda_ = 0.05 # Paper Travascio
#Permeabilität
K = 1.26e-3; #[mm^4/Ns)]
# allgemeiner initialer Druck im System
p_0 = 0.5 # MPa, wird auf 0 gesetzt (Meniskus hat in Anfangsposition nur durch den Fluiddruck diese Form)
## Definieren von Zeitparametern
T = 10.0 # final time
num_steps = 100 # number of time steps
dt = T / num_steps # time step size
## Größen, die im Variationsproblem benötigt werden und als Konstanten gesetzt werden müssen
p_0 = Constant(p_0)
K = Constant(K)
mu = Constant(mu)
lambda_ = Constant(lambda_)
##############################################################################
# Netz definieren #
##############################################################################
Netz = Mesh("Meniskus.xml")
##############################################################################
# Definition mixed elements #
##############################################################################
#Mixed Elements
P2 = VectorElement('Lagrange', Netz.ufl_cell(), 2) # quadratische Ansätze für Verschiebung
P1 = FiniteElement('Lagrange', Netz.ufl_cell(), 1) # lineare Ansätze für Druck
element = P2 * P1 #mixed element
# Define FunctionSpace
V = FunctionSpace(Netz, element) #mixed FunctionSpace für mixed element
##############################################################################
# Ränder und Randbedinungen #
##############################################################################
# Toleranz
tol = 1E-14
## gesamte Ränder
def boundary(x, on_boundary): # für Druck-RB
return on_boundary
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", Netz, "Meniskus_facet_region.xml")
# Nummer 1: unterer Rand
# Nummer 2: rechter Rand
# Nummer 3: oberer Rand
# für Randintegral
ds_test = Measure('ds', domain=Netz, subdomain_data=boundaries)
# für Normalenvektor
n = FacetNormal(Netz)
# Verschiebungsrandbedingung
bc_u1 = DirichletBC(V.sub(0), Constant((0.0, 0.0)), boundaries, 1)
# Druckrandbedingungen
bc_p1 = DirichletBC(V.sub(1), Constant(p_0), boundaries, 1)
bc_p3 = DirichletBC(V.sub(1), Constant(p_0), boundaries, 3)
# Class representing the intial conditions/Anfangsbedingungen
class InitialConditions(UserExpression):
def eval(self, values, x):
values[0] = 0 # ux
values[1] = 0 # uy
values[2] = p_0 # p; ist überhaupt Anfangsdruck des Fluids vorhanden, wenn es in Ruhe ist?
def value_shape(self):
return (3,)
t = Constant((0.0000 , 0.0005))
##############################################################################
# Definition wichtiger Größen #
##############################################################################
#Test- und Trial-Functions
(u, p) = TrialFunctions(V) # Größen u und p zum neuen Zeitschritt n+1
(v_1, v_2) = TestFunctions(V) # Wichtungsfunktionen
x_n = Function(V) # x: Unbekannte zum Zeitpunkt n; Hilfsgröße zum Splitten
u_n, p_n = split(x_n) # u, p zum Zeitpunkt n
# Dimension (1D/2D/3D)
d = u.geometric_dimension() # geometrische Dimension
I = Identity(d) # Einheitstensor
## Definieren von Spannung und Verzerrung für Festkörperanteil
# Kleine Verzerrungen
def epsilon(u):
return sym(nabla_grad(u))
# Gesamtspannung aus Anteilen des Festkörpers und des Fluids
def sigma_ges(u,p):
return lambda_*nabla_div(u)*I + 2*mu*epsilon(u) - p*I
## Variationsproblem
# Festkörperfeld: Impulsbilanz für das Gemisch
# Fluidfeld: Volumenbilanz für das Gemisch
Imp = inner(sigma_ges(u,p), nabla_grad(v_1)) * dx \
- dot(t, v_1) * ds_test(3)
Vol = - dot(((u-u_n)/dt), nabla_grad(v_2)) * dx \
+ K * dot(nabla_grad(p), nabla_grad(v_2))*dx
# Zusammgesetzte Gleichung
F = Imp + Vol
# Aufteilen der Gleichung in linke und rechte Seite
l = lhs(F)
r = rhs(F)
# Erstellen von VTK-files
vtkfile_u = File("Ergebnisse/u_{}.pvd".format(Titel))
vtkfile_p = File("Ergebnisse/p_{}.pvd".format(Titel))
##############################################################################
# Lösungsschleife #
##############################################################################
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
#gesamten Randbedingungen
bc_ges = [bc_u1, bc_p1, bc_p3]
# Solve variational problem for time step
w = Function(V)
solve(l == r, w, bc_ges)
(u,p) = w.split(True)
# für Paraview
u.rename('u','u') # Das ist notwendig, damit Paraview nicht jeden
p.rename('p','p') # Zeitschritt als eine neue Variable ansieht.
# Save solution to file (VTK)
vtkfile_u << (u, t)
vtkfile_p << (p, t)
# Update previous solution
assign(x_n, [u,p])
end = time.time()
ellapsed = end-start
