1D Heat Equation that is Nonlinear in one subdomain and linear in another

variational formulation
1

Hi, I’m a new user. I want to model a 1D heat equation with time varying Dirichlet BC and two subdomains: one that has a temperature dependent thermal conductivity, the other has a constant thermal conductivity. I have tried modifying the example problems based on the heat equation with subdomains and nonlinear Poisson equation, but I haven’t been able to successfully combine them. Here is what I have so far:

from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *
import numpy as np
import sympy as sym

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 5         # parameter beta


tol = DOLFIN_EPS
def k_0(u):
    return 1+2*u

class K(Expression):
    def set_k_values(self, k_0, k_1):
        self.k_0, self.k_1 = k_0, k_1

    def eval(self, value, x, u):
        "Set value[0] to value at point x"
        tol = 1E-14
        if x[0] <= 0.5 + tol:
            value[0] = self.k_0
        else:
            value[0] = self.k_1

k_1 = 0.01

# Attempt at using C++ code snippet
#kappa = Expression('x[0] <= 0.5 + tol ? 1+1*u : k_1', degree=0,
#tol=tol, k_1=k_1, u=0)

kappa = K()
kappa.set_k_values(k_0,k_1)


# Create mesh and define function space
meshpoints = 50
my_mesh = UnitIntervalMesh ( meshpoints )
V = FunctionSpace(my_mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + beta*t',beta=beta, degree=1, t=0)

def boundary ( x ):
    value = x[0] < DOLFIN_EPS
    return value

bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + kappa*dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx

# Time-stepping

t = 0
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t
    # Compute solution
    solve(F==0, u, bc)

    # Plot solution
    plot(u)

    # Compute error at vertices
    # u_e = interpolate(u_D, V)
    # error = np.abs(u_e.vector().array() - u.vector().array()).max()
    # print('t = %.2f: error = %.3g' % (t, error))

    # Update previous solution
    u_n.assign(u)
2

Ok I think I have figured this out myself. I can just use a C++ expression after I have defined u as a function:

...
# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)
a = 0.01
b = 0
kappa = Expression('x[0] <= 0.5 + tol ? b+a*u : k_1', degree=0,
tol=tol, k_1=k_1, u=u, a=a, b=b)