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

83

views

0

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 :

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.

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

Thank you very much in advance for your help !

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)
F=c*vc*dx + dt*D*dot(grad(c), grad(vc))*dx - c_n*vc*dx + g*vc*ds
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

### 1 Answer

3

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)
#F=c*vc*dx + dt*D*dot(grad(c), grad(vc))*dx - c_n*vc*dx + g*vc*ds
#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))
F=((c-c_n)/dt)*vc*dx + D*dot(grad(c), grad(vc))*dx + g*vc*ds
### 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).

Please login to add an answer/comment or follow this question.

Thank you very much for your answer. It seems to work, though I know have to find a way to make it converge ! :)