# Non-linear coupled equations are not converging

350
views
1
12 months ago by
Hello everyone
I am trying to simulate non-linear coupled equations and this is my code:

# -*- coding: utf8 -*-

from __future__ import print_function
from fenics import *
import numpy as np
from dolfin import *
from mshr import *
tol = 1E-14
t = 0.0

class InitialConditions(Expression):
def eval(self, values, x):
values[0] = 0.0
values[1] = 315
if  (near(x[0], 2.0, 0.2) or near(x[0], 13.0, 0.2)) and \
near(x[1], 1.0, 0.1):
values[0] = 5000.0
values[1] = 315.01
def value_shape(self):
return (2,)

T = 1.0e-1    # final time
num_steps = 100    # number of time steps
dt = T / num_steps # time step size
D = 1.000e-10    # Diffusion Coefficient
K = 8.000e-1       # Thermometric diffusion ratio
c = 4.000e3            # Heat Capacity
X = 2.000e-7          # Thermoemtric diffusion coefficient

# Create mesh and define function space

mesh = RectangleMesh(Point(0., 0.),Point(15.0, 2.0), 30, 30)
P1 = FiniteElement("P", triangle,2)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define test functions
v_1, v_2 = TestFunctions(V)
# Define functions for concentration and Temperature
u = Function(V)
u_n = Function(V)

# Split system functions to access components
u_1, u_2 = split(u)
# Define boundary condition

z_1 = Expression(('0.0','315.0'), degree = 2)
z_2 = Expression(('0.0','315.0'), degree = 2)

#tol = 1E-14
def boundary_D(x, on_boundary):
return on_boundary and not ((near(x[0], 2.0, 0.5) or near(x[0], 13.0, 0.5)) and\
near(x[1], 2.0, tol)) and not (near(x[0], 0.0, tol) or\
near(x[0], 15.0, tol))
def boundary_D_source(x, on_boundary):
return on_boundary and (near(x[0], 2.0, 0.5) or near(x[0], 13.0, 0.5)) and \
near(x[1], 2.0, tol)
bc_1= DirichletBC(V, z_1, boundary_D_source)
bc_2 = DirichletBC(V, z_2, boundary_D)

bcs = [bc_1, bc_2]

#bcs = bc_2
# Define initial value
f = InitialConditions(degree = 1)
u_n = interpolate(f, V)
S_1 = Source(degree = 2)
u_n1, u_n2 = u_n.split()

# Define variational problem
v_1, v_2= TestFunctions(V)
k = Constant(dt)
K = Constant(K)
D = Constant(D)
c = Constant(c)
X = Constant(X)

g = Constant(-1.0)
r = Constant(-1.00e-6)

+ u_1*((u_2 - u_n2))*v_2*dx  - (K/c)*(u_2)*((u_1 - u_n1))*v_2*dx\

vtkfile_1 = File("Project/u1.pvd")
vtkfile_2 = File("Project/u2.pvd")
du = TrialFunction(V)

# Time-stepping
for n in range(num_steps):

_u_1, _u_2 = u_n.split()
vtkfile_2 << (_u_1,t)

# Compute solution
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)

solver.solve()
# Update previous solution
u_n.assign(u)
t += dt

​

The problem is with changing coefficients or time, either is not converging or it is converging to non-sense(value of u_1 for example becomes zero in whole domain).

What is problem in your opinion?!
Initial condition seems reasonable. Sometimes I think problem is variational formulation, but I am not sure.

These are the equations:
$\frac{\partial u_1}{\partial t}=D\times\left(\Delta u_1+\frac{K}{u_2}\times\Delta u_2\right)\backslash\frac{\partial u_2}{\partial t}-\left(\frac{K}{c}\right)\times\left(\frac{u_2}{u_1}\right)\times\frac{\partial u_1}{\partial t}=X\times\Delta u_2$u1t =D×(Δu1+Ku2 ×Δu2)\u2t (Kc )×(u2u1 )×u1t =X×Δu2
I cannot reproduce your problem since I get a different error in line 70 NameError: name 'Source' is not defined.
written 11 months ago by Eleonora Piersanti
Oh I am sorry. Delete that line. I forgot to delete it, it does not affect the code. Thank you for checking my code by the way.
written 11 months ago by Morteza

0
11 months ago by
When dealing with such low values for coefficients such as the diffusion and thermometric diffusion maybe it would be a good idea to nondimensionalize your equation. We need to be careful with too little and too big values as we are dealing with float points arithmetic.