# Problem with crank nicolson for coupled non linear pde

hi at all,

I started with FEA last December and my problem (thermoelectrcity) is not the best way to start with finite elements. But with the fenics and dolfin documentation and with some examples from the web I was able to solve the steady state problem.

Now I have no idea how I should implement the CN-method to solve the coupled problem in transient mode.

I used the cahn hilliard example for my problem but I do not get a correct time progress.

The solution is only a switching between the initial condition and the final time step.

Regardless of the time step size.

I think that the error is in the weak form and/or in the while loop.

But due to the lack of knowledge I find no way to solve my problem with the help of other documented examples.

I hope, somebody has time to look over my code and give me a right hint.

I had no idea how to reduce the problem into few lines.

I'm sorry. But maybe the problem is interesting enough to be used as a demo file. :)

best regards

Michael

Below the code:

import numpy

from dolfin import *

# Class representing the intial conditions

class InitialConditio

def __init_

def eval(self,vec,x):

def value_shape(self):

return (2,)

# Class for interfacing with the Newton solver

class ThermoelEquatio

def __init__(self, a, L,bc):

self.L = L

self.a = a

self.bc=bcs

def F(self, b, x):

for condition in self.bc:

def J(self, A, x):

for condition in self.bc:

# Model parameters

dt =5.0e-02 # time step

theta = 1# time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

Tinni = 300. # initial temperature

Vinni = 0.0 # initial electr. potential

# Form compiler options

parameters[

parameters[

parameters[

# Create mesh and define function spaces

# ===============

# Mesh of simple 2 dim. pellet

# ( 0.1 cm X 0.6 cm)

# ===============

mesh = RectangleMesh(

#######

#Subdomains

#three domains with different temperature depend. material properties

#######

class OmegaTop(

def inside(self, x, on_boundary):

return (between(x[1], (0.59, 0.6)) and between(x[0], (0.0, 0.1)))

class OmegaMiddle(

def inside(self, x, on_boundary):

return (between(x[1], (0.01, 0.59)) and between(x[0], (0.0, 0.1)))

class OmegaBottom(

def inside(self, x, on_boundary):

return (between(x[1], (0.0, 0.01)) and between(x[0], (0.0, 0.1)))

# Initialize sub-domain instances

omegatop = OmegaTop()

omegamiddle = OmegaMiddle()

omegabottom = OmegaBottom()

# Initialize mesh function for interior domains

subdomains = CellFunction(

subdomains.

omegatop.

omegamiddle.

omegabottom.

#plot(subdomain

#######

# boundary conditions

#######

class GammaBottom(

def inside(self, x, on_boundary):

return near(x[1], 0.0)

class GammaTop(

def inside(self, x, on_boundary):

return near(x[1], 0.6)

gammabottom =GammaBottom()

gammatop = GammaTop()

# Initialize mesh function for boundary domains

boundaries = FacetFunction(

boundaries.

gammabottom.

gammatop.

V1 = FunctionSpace(mesh, 'DG', 0)# temperatur depend. properties

P1 = FunctionSpace(mesh, "Lagrange", 1)#temperature

P2 = FunctionSpace(mesh, "Lagrange", 1)#electric potential

ME = P1*P2

# Define trial and test functions

du = TrialFunction(ME)

te, vo = TestFunctions(ME)

# Define functions

u = Function(ME) # current solution

u0 = Function(ME) # solution from previous converged step

# Split mixed functions

dc, dp = split(du)

temp, volt = split(u)

temp_prev, volt_prev = split(u0)#previous step

# Create intial conditions and interpolate

u_init = InitialConditio

u.interpolate(

u0.interpolate(

#I do not know if I'm right

temp_mid = (1.0-theta)

volt_mid = (1.0-theta)

#######

#Temperatur dependend material properties

#Coefficients of the polynomials A * temp² + B* temp + C

#######

#thermal conductivity

hcond1= Function(V1)

hcond2= Function(V1)

hcond3 = Function(V1)

#Seebeck coeff.

seeb1= Function(V1)

seeb2= Function(V1)

seeb3 = Function(V1)

#electr. conductivity

elcond1= Function(V1)

elcond2= Function(V1)

elcond3 = Function(V1)

k1_values = [0., 1.68181805028e-

k2_values = [0., -1.11448469117e

k3_values = [4., 0.0189634948598

a1_values = [0., -5.5012121918e-

a2_values = [0., 6.92678311478e-07 ,0.]

a3_values = [6.5e-6, -4.45128316556e

c1_values = [ 0.,0.0009090909

c2_values = [ 0.,-1.71636366069, 0.]

c3_values = [ 5.9e6,1603.

help = numpy.asarray(

hcond1.vector()[:] = numpy.choose(help, k1_values)

hcond2.vector()[:] = numpy.choose(help, k2_values)

hcond3.vector()[:] = numpy.choose(help, k3_values)

elcond1.vector()[:] = numpy.choose(help, c1_values)

elcond2.vector()[:] = numpy.choose(help, c2_values)

elcond3.vector()[:] = numpy.choose(help, c3_values)

seeb1.vector()[:] = numpy.choose(help, a1_values)

seeb2.vector()[:] = numpy.choose(help, a2_values)

seeb3.vector()[:] = numpy.choose(help, a3_values)

#######

rho=7.8 #Density g/cm³

cp=0.28 #heat capacity J/gK

Tc = Constant(300.15) # cold side temperature (bottom boundary)

V0 = Constant(0.0) # el. potential (bottom boundary)

Th = Constant(873.15) # hot side temperature (top boundary)

# Th could alsovary in time like:

#Th = Expression('400. +800.*sin(

#V1 = Constant(0.2) # el. potential (top boundary)

#### Define boundary conditions

bc0 = DirichletBC(

bc1 = DirichletBC(

bc2 = DirichletBC(

#addionally

#bc4 = DirichletBC(

#and/ or Neumann:

n = FacetNormal(mesh)

qc = as_vector(

jc = as_vector(

g_T = dot(n,qc)

g_R = dot(n,jc)

bcs = [bc0,bc1,

ds = Measure(

#electric current flux per unit area A/cm²:

j = -(elcond1*

- (seeb1*

elcond2*

#thermal flux per unit area W/cm²:

q = -((hcond1*

+ (seeb1*

#electric field V/cm

E = -nabla_

# generated heat W/cm³ (inner source)

Q = dot(j,E)

# Weak statement of the equations

#VOLT: -inner(

#TEMP: -inner(

L0 = -inner(

L1 = -(rho*cp)

L = L0 +L1

# Compute directional derivative about u in the direction of du (Jacobian)

a = derivative(L, u, du)

# Create nonlinear problem and Newton solver

problem = ThermoelEquation(a, L,bcs)

solver = NewtonSolver("lu")

solver.

solver.

# Output file

#fileU=

#fileT=

#######

# Step in time

t = 0.0

T = 20*dt

while (t < T):

print "======

print "Solving at t= "+str(t)+"s, tmax="+str(T)+\

"s, dt="+str(dt)+"s."

print "======

t += dt

u0.vector()[:] = u.vector()

solver.

#fileT << (u.split()[0], t)

#fileU << (u.split()[1], t)

plot(u.

plot(u.

interactive()

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- M Waste

- Solved:

- Last query:

- Last reply: