### (Deleted) Neumann boundary condition does't work

18
views
-1
9 weeks ago by
Hi!

I am trying to solve a stress-diffusion coupled problem with a certain inward flux on the boundary of a circular plate.

I am getting valid results when I use a concentration boundary condition (Dirichlet) but the flux boundary condition gives out 0 as the result, implying there is no inward flux.

Following is the code.
from __future__ import division
from fenics import *
from mshr import *
from math import sin, cos, pi
import numpy as np
from matplotlib import pyplot

#DEFINING CONSTANTS

nu = 0.28
R = 8.31
T = 273
R0 = 2.5
D = 1e-16
Cmax = 3.67e5
omg = 8.18e-6
E = 90.0e9*omg/(R*T)

T = 5e-1            # final time
num_steps = 100       # number of time steps (100)
dt = T / num_steps # time step size
print(1)

#DEFINING BOUDNARIES

sphere = Circle(Point(0,0),R0)
mesh = generate_mesh(sphere,30)
print(mesh.num_vertices())
print(mesh.topology().dim())

def centre(x,on_boundary):
return near(x[0],0) and near(x[1],0)

def boundary(x,on_boundary):
return on_boundary
print(2)

#DEFINING FUNCTION SPACES

P1 = VectorElement('CG', triangle, 2)
P2 = FiniteElement('CG', triangle, 2)
element = MixedElement([P1,P2])
W = FunctionSpace(mesh, element)

sol = Function(W)
(u, c) = split(sol)
(v, w) = TestFunctions(W)
print(3)

#DEFINING BOUNDARY CONDITIONS

bcc = DirichletBC(W.sub(0),(Constant(0.0),Constant(0.0)),centre,method='pointwise')
print(4)

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundary_markers.set_all(0)
class BoundaryX0(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 1)
meshfile=File('mesh.pvd')
meshfile << boundary_markers

bc = [bcc]

print(5)

#DEFINING TERMS

def F(u,c):

def Fi(c):
return ((1+c)**(1/3))*Identity(2)

def Fe(u,c):
return F(u,c)*inv(Fi(c))

def E_gl(u,c):
return 0.5*(F(u,c).T*F(u,c) - Identity(2))

def Ee_gl(u):
return 0.5*(Fe(u,c).T*Fe(u,c) - Identity(2))

def We(u,c):
return det(Fi(c))*E*0.5*(nu*(tr(Ee_gl(u))*tr(Ee_gl(u)))/\
(1-2*nu)+tr(Ee_gl(u)*Ee_gl(u)))/(1+nu)

def P(u):
return E*(nu*tr(Ee_gl(u))*Identity(2)/(1-2*nu)+Ee_gl(u))*Fe(u,c)/(1+nu)

def sigma(u,c):
return P(u)*F(u,c).T/det(F(u,c))

def sigmah(u,c):
return (1-nu)*(tr(sigma(u,c))/3)

k = Constant(dt)

kl = 1e-7
kt = (kl**2)/D
k3 = (inv(F(u,c).T))**2
k4 = 2*k3/3
kn = nu/(1-2*nu)
kj = kl/(kt*omg)

j0 = Constant(kj*1e-4)

#DEFINING INITIAL CONDITIONS

c_n=project(Constant(1.0),W.sub(1).collapse())

F2 = inner(-j0,w)*ds(1)+((c-c_n)/k)*w*dx
F = F1+F2+F3+F4

#FILE WRITIN
vtkfile_u = File('Trial_7/displacement.pvd')
vtkfile_c = File('Trial_7/concentration.pvd')
vtkfile_s = open('Trial_7/soc.txt','a+')

#TIME-STEPPING

t = 0
for n_step in range(num_steps):

# Update current time
t += dt
print(t)

jacobian = derivative(F,sol)

solve(F==0, sol, bc, J=jacobian)
print(10)

count = 0

print("Number of steps: ", n_step)

if (n_step%1) == 0:
(u,c) = sol.split()
assign(c_n,c)
vtkfile_u << sol
vtkfile_c << c
soc = assemble(c_n*dx)*100/(pi*R0**2)
print(soc)
vtkfile_s.write('%f '%soc)
vtkfile_s.write('%f \n'%t) ​

The first term of F2 is the Neumann boundary condition, where j0 is the magnitude of the flux.
I thought it might be leading to a trivial solution to the equation and tried with varying values of c_n (which is the initial concentration). But even then, it is solving only for 1 step and in the subsequent steps the Newton solver is running for 0 iterations.

Any help or pointers on this is appreciated.