(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):
    return Identity(2) + grad(u)

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

F1 = -inner(P(u),grad(v))*dx
F2 = inner(-j0,w)*ds(1)+((c-c_n)/k)*w*dx
F3 = k3*inner(grad(c),grad(w))*dx
F4 =-k3*inner(grad(sigmah(u,c)),grad(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.
Thanks in advance.

~Raghunandan
Community: FEniCS Project
hi raghunandan,
check the code  once again,
Thanks
written 9 weeks ago by Md Ayub Ansari  
Please login to add an answer/comment or follow this question.
The thread is closed. No new answer/comment may be added.

Similar posts:
Search »
  • Nothing matches yet.