Evaluating discontinuous derivatives in internal boundary (electric stress)


255
views
0
6 months ago by
Hello, I want to solve a multiphysics problem.

The first problem I want to solve is the Poisson equation for the electric field in a 2D domain formed by two subdomains  $\Omega_1$Ω1,  $\Omega_2$Ω2  separated by an internal interface  $\Gamma$Γ . The electric permeabilities are  $\epsilon_1$ϵ1  and  $\epsilon_2$ϵ2  respectively. The only sources of charge are a linear density of charge  $\sigma$σ  at the boundary  $\Gamma$Γ . Therefore, I will solve the problem using the laplace equation in the interiors of the domain  $\Omega_1$Ω1 and  $\Omega_2$Ω2  That is:
\(\nabla ^2 V = 0\)
The potential field is continuous across the domain, but the electric field \(\nabla V\) has a discontinuity at \(\Gamma\)

The boundary conditions at the interface are something like \(\epsilon_1 \nabla V_1 · \vec{n} - \epsilon_2 \nabla V_2 · \vec{n} = \sigma\) and \(\nabla V_1 · \vec{t} =  \nabla V_2 · \vec{t}\)

The second problem I want to solve are the Stokes fluid equations in the \(\Omega_1\) domain. The boundary conditions will depend on the electric stress at the boundary. The electric stress depends on \(E^n_{\Omega_1}, E^n_{\Omega_2}\).

My question is how can I compute this discontinuous values of the electric field at each points of the boundary, that is, evaluating \(\nabla V · \vec{n}\) at each side of the boundary. I have tried to solve the problem. My code is shown below. The problem is that I don't know how to tell FEniCS to compute the derivatives at both of the parts of the boundary. The code I did so far gives NaN for the flux ..., also, is there a way to calculate the individual values of the electric field at the boundaries and not the flux?

    V = FunctionSpace(mesh, "CG", 1)
    dS = Measure("dS", subdomain_data=boundaries)
    dx = Measure("dx", subdomain_data=domains)
    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    sigma = Constant(4.0)
    a = e1*dot(grad(u), grad(v)) * dx(1) + e2*dot(grad(u), grad(v)) * dx(2)
    L = - sigma * v('+') * dS(1)
    # Define boundary condition values
    u1 = Constant(0)
    u2 = Constant(-14)

    # Define boundary conditions
    # bc1 = DirichletBC(V, u1, boundaries, 0) 
    bc1 = DirichletBC(V, u1, boundaries, 2)
    bc2 = DirichletBC(V, u2, boundaries, 5)

    # Compute solution
    phi = Function(V)
    solve(a == L, phi, [bc1, bc2])

#Evaluate derivatives
    w = TrialFunction(W)
    v = TestFunction(W)
    grad_phi1 = Function(W)
    grad_phi2 = Function(W)
    a1 = inner(w,v)*dx(1)
    L1 = inner(grad(phi),v)*dx(1)
    a2 = inner(w,v)*dx(2)
    L2 = inner(grad(phi),v)*dx(2)
    solve(a1 == L1, grad_phi1)
    solve(a2 == L2, grad_phi2)​
#Evaluate normals at mesh
    n = FacetNormal(mesh)
    dS = Measure("dS", domain = mesh, subdomain_data = boundaries)
    flux1 = assemble(dot(grad_phi1,n)*dS(1))
    flux2 = assemble(dot(grad_phi2,n)*dS(1))
Community: FEniCS Project

1 Answer


2
6 months ago by
Hi,
1. Have you marked your boundaries and your subdomains?

2.
 assemble(dot(grad_phi1,n)*dS(1))

is the integral of something over dS(1). So you are computing the total flux coming through the boundary, which is a number. Not a variable function over the boundary. I guess you want to compute dot( grad_phi1, n ) on every node and not the integral of it.

3. Why do you solve the a1==L1 and a2==L2 PDEs? You are projecting grad(phi) onto the same functionspace that it already lives. With this you are good to go.
grad(phi)*dx(1) or dx(2)​


4. Your boundary is one. There is no such thing as the two sides of the boundary. You want to compute the values of this derivative of phi on the nodes that live on the boundary. Now the only thing that changes is the normal vector according to your geometry and your subdomains. Maybe you should try to define the normal vector of the subdomain, because the way you are using it is pointless.

Consider this:

from fenics import *
mesh    = UnitSquareMesh(20, 20)

V       = VectorFunctionSpace(mesh, 'P', 1)
n       = FacetNormal(mesh)

u       = TrialFunction(V)
v       = TestFunction(V)

a = inner(u,v) * ds
l = inner(n,v) * ds

A = assemble(a, keep_diagonal=True)
A.ident_zeros() 
b = assemble(l)

n_ = Function(V)

solve(A, n_.vector(), b)


This is a way to obtain the normal vector on any measured boundary ds(0,1....) and not on the whole boundary as FacetNormal(mesh) will do . I guess this should solve your problems.

Hi Minas,

Thank you very much for your response! Let me ask you a few questions more:

1. Have you marked your boundaries and your subdomains?

Yes, I did. I marked the interface boundary as 1. That is dS(1).

2. assemble(dot(grad_phi1,n)*dS(1))

is the integral of something over dS(1). So you are computing the total flux coming through the boundary, which is a number. Not a variable function over the boundary. I guess you want to compute dot( grad_phi1, n ) on every node and not the integral of it.

---------------Yes! Do you know how can I evaluate this function at a Point = {x, f(x)}, where f is the interface boundary? The boundary is defined as a single line, and query Points for evaluating the normal vector could not lie exactly at the boundary. How can I "extrapolate" ?

3. Why do you solve the a1==L1 and a2==L2 PDEs? You are projecting grad(phi) onto the same functionspace that it already lives. With this you are good to go.grad(phi)*dx(1) or dx(2)​

---------------Yes you are right!  At the beginning I was doing exactly what you said, the problem is that I am using Paraview to visualize and I cannot open the results with integrals performed over subdomains, I had to do it separately.

4. Your boundary is one. There is no such thing as the two sides of the boundary. You want to compute the values of this derivative of phi on the nodes that live on the boundary. Now the only thing that changes is the normal vector according to your geometry and your subdomains. Maybe you should try to define the normal vector of the subdomain, because the way you are using it is pointless.

---------------The problem is that the gradient is discontinuous just at the interface boundary. I want to compute the values of the derivatives of Phi at the nodes Immediately above the interface boundary (region 1) and immediately below the interface boundary (region 2) because at the interface boundary, the gradents don't exist (there is a jump). Once I compute the normal vectors using the code you posted, how can I obtain the different values of the gradients across the interface boundary?
written 6 months ago by Ximo GC  
Hi,

I would like to see how you have defined the subdomains and the interface and also how they are updated in time as the interface moves.

written 6 months ago by Minas Mattheakis  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »