### Neumann BC depending on solution (time dependent)

83
views
0
8 weeks ago by
Hi all,

Here's the problem I want to solve :

dC/dt=D*lap(C) in the domain
dC/dx=g=7e-6 C_surf ^0,74 on the boundary

The weak form goes like this :
F=c*vc*dx + dt*D*dot(grad(c), grad(vc))*dx - c_n*vc*dx + g*vc*ds​

The issue I'm facing is in the updating of the Neumann condition in the time stepping. I found ways to update this equation when g is an Expression depending on t. The thing is it's actually depending on the solution itself.

g = Function(V)
g = conditional(gt(c_n, 0), 7e-6*c_n**0.74, Constant(0.0))​

Rem

Here's my code

from fenics import *
from dolfin import *
import numpy as np

print('Defining the solving parameters')
Time = 6000000.0         # final time
num_steps = 10    # number of time steps
dt = Time / num_steps # time step size
t=0 #Initialising time to 0s

print('Defining mesh')
# Create mesh and define function space
mesh = Mesh('geo/mesh_rcb_3_layers_not_refined.xml')

concrete_thickness=240e-3
steel_thickness=2e-3
polymer_thickness=20e-3
internal_cavity_dimension=1738e-3

print('Defining Functionspaces')
V = FunctionSpace(mesh, 'P', 1) #FunctionSpace of the solution c
V0 = FunctionSpace(mesh, 'DG', 0) #FunctionSpace of the materials properties

### Define initial values
print('Defining initial values')
##Tritium concentration
iniC = Expression('0',degree=1)
c_n = interpolate(iniC, V)

### Boundary Conditions
print('Defining boundary conditions')

inside =  CompiledSubDomain("(fabs(x[0])<=side && fabs(x[1])<side && fabs(x[2])<side) && on_boundary", side = (internal_cavity_dimension+steel_thickness)/2)
outside = CompiledSubDomain("(fabs(x[0])>=side || fabs(x[1])>=side || fabs(x[2])>=side) && on_boundary", side = (internal_cavity_dimension+concrete_thickness+polymer_thickness+steel_thickness)/2)

##Tritium concentration
inside_bc_c=1
outside_bc_c=0
bci_c=DirichletBC(V,inside_bc_c,inside)
bco_c=DirichletBC(V,outside_bc_c,outside)
g = Function(V)
g = conditional(gt(c_n, 0), 7e-6*c_n**0.74, Constant(0.0))#
#g = Expression('-t', t=0,degree=2)
bcs_c=[bci_c]

##Temperature
inside_bc_T=1
outside_bc_T=Expression('14+cos(2*3.14*t/325/24/3600)+cos(2*3.14*t/24/3600)', t=0, degree=2)
bci_T=DirichletBC(V,inside_bc_T,inside)
bco_T=DirichletBC(V,outside_bc_T,outside)
bcs_T=[bci_T]

###Defining materials properties
print('Defining the materials properties')

D  =1e-6 #Diffusion coefficient

### Define variational problem
print('Defining the variational problem')

c = TrialFunction(V)
vc = TestFunction(V)
ac, Lc = lhs(F), rhs(F) #Rearranging the equation
AC=assemble(ac)

### Time-stepping
c = Function(V)

fileC = File("solutions/solutionCpourri/solutionC.pvd")

for n in range(num_steps):
# Update current time
print("t= "+str(t)+" s")
print(str(100*t/Time)+" %")
t += dt
# Compute solution concentration

#solve(ac == Lc, c, bcs_c)
#g.t = t + dt
LC=assemble(Lc)
bci_c.apply(AC,LC)
solve(AC,c.vector(),LC)
fileC << (c,t)

# Update previous solution
c_n.assign(c)


Community: FEniCS Project

3
8 weeks ago by
If I understand correctly, it looks like the code already implements a "lagged" version of the desired BC, in which $g$ depends on the previous time step's solution of $c$ (although I believe there should be a factor of $\Delta t$ multiplying $g$ somewhere).  It is possible to include the current time step's $c$ in the residual instead, for a fully-implicit approach, but this makes the problem nonlinear.  You can still try to solve it with something like
### Define variational problem
print('Defining the variational problem')

#c = TrialFunction(V)
c = Function(V)
vc = TestFunction(V)
#ac, Lc = lhs(F), rhs(F) #Rearranging the equation
#AC=assemble(ac)

# (Note: Grouped time derivative terms and effectively multiplied g by dt
# relative to F defined in comments above.)
g = conditional(gt(c, 0), 7e-6*c**0.74, Constant(0.0))

### Time-stepping
#c = Function(V)

fileC = File("solutions/solutionCpourri/solutionC.pvd")

for n in range(num_steps):
# Update current time
print("t= "+str(t)+" s")
print(str(100*t/Time)+" %")
t += dt
# Compute solution concentration

#solve(ac == Lc, c, bcs_c)
#g.t = t + dt
#LC=assemble(Lc)
#bci_c.apply(AC,LC)
#solve(AC,c.vector(),LC)

# Notation to invoke nonlinear solver:
solve(F==0,c,J=derivative(F,c),bcs=[bci_c,])

fileC << (c,t)

# Update previous solution
c_n.assign(c)
​

However, convergence of Newton's method is not guaranteed (especially for a non-smooth residual).

Hi,

Thank you very much for your answer. It seems to work, though I know have to find a way to make it converge ! :)
written 8 weeks ago by Rémi Delaporte-Mathurin