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


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
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)
#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). 

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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »