Stress vector on a curved boundary in a 2D coupled field problem

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

Hi,

I’m quite new to fenics as well, and therefore I probably can’t help your problem. Maybe it’s due to the initital pressure. I had some trouble with that as well.
But I would like to ask you about your variational problem. I also work in the framework of the Theory of Porous Media and I miss the volume fractions in your equations. I would write the variational problem for a two phase media somehow like this:

Imp = inner(sigma_ges(u,p), nabla_grad(du_s)) * dx - dot(b, du_s) * dx - dot(t, du_s) * ds_test(3)

Vol = (nS-nS_n) / dt + dot(((u-u_n)/dt), nabla_grad(dpI)) * dx - K * dot(nabla_grad( p ), nabla_grad(dpI))*dx #- dot(tilde_v, nabla_grad(dpI)) * ds_test(3)

In the beginning you find the volume fraction of the fluid constituent.

Hi mrBrown,

In the picture you can see the field equation and the way to the variational form of the second equation. I don’t see a missing volume fraction of the fluid constituent.

Please let me know, if I did any mistake. Thank you!

Hi Supporter53,

thanks, that looks fine, I think we just ran in a missunderstanding. When I read “Theory of Porous Media”, I thought of the approach discussed in here: Theory of Porous Media - Highlights in Historical Development and Current State | Reint de Boer | Springer.
Have you been able to get over your problem? Maybe someone else can help?

Hi mrBrown,

yes, I think it depends on the model that you use for your calculations - so it’s probably a different approach :slight_smile:
I was able to get over this, it was a problem in connection with the boundary conditions. Now it’s working, thank you!