Fenisc fiber inclusion thermomechanical problem: two materials

gistfile1.py

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt


# https://comet-fenics.readthedocs.io/en/latest/demo/thermoelasticity/thermoelasticity_transient.html#References

L = 1.
R = 0.4
N = 50 # mesh density


domain = Rectangle(dolfin.Point(0., 0.), dolfin.Point(L, L))
# subdomain 0 = matrix
domain.set_subdomain(1, Rectangle(dolfin.Point(0., 0.), dolfin.Point(L, L)))
# subdomain 1 = fiber
domain.set_subdomain(2, Circle(Point(0., 0.), R) )

mesh = generate_mesh(domain, N)

dolfin.info(domain, True)
dolfin.plot(mesh, "2D mesh")
subdomains = dolfin.MeshFunction("size_t", mesh, 2, mesh.domains())
dolfin.plot(subdomains, "Subd")
plt.show()

T0 = Constant(70.)
#DThole = Constant(10.)
Em, num = 2500., 0.46
Ef, nuf = 500., 0.36 
alphaf, alpham = 129.0E-6, 61.2E-6
class Young(UserExpression):
    def __init__(self, subdomains, E_0, E_1, **kwargs):
        self.subdomains = subdomains
        self.E_0 = E_0
        self.E_1 = E_1
        super().__init__(**kwargs)
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.E_0
        else:
            values[0] = self.E_1
    def value_shape(self):
        return (2,)

class Poisson(UserExpression):
    def __init__(self, subdomains, E_0, E_1, **kwargs):
        self.subdomains = subdomains
        self.E_0 = E_0
        self.E_1 = E_1
        super().__init__(**kwargs)
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.E_0
        else:
            values[0] = self.E_1
    def value_shape(self):
        return (2,)

class CTE(UserExpression):
    def __init__(self, subdomains, E_0, E_1, **kwargs):
        self.subdomains = subdomains
        self.E_0 = E_0
        self.E_1 = E_1
        super().__init__(**kwargs)
    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.E_0
        else:
            values[0] = self.E_1
    def value_shape(self):
        return (2,)
E = Young(subdomains, Em, Ef, degree=0)
nu = Poisson(subdomains, num, nuf, degree=0)
alpha = CTE(subdomains, alphaf, alpham, degree=0 )

nu = 0.3
lmbda  = UserExpression(E*nu/((1+nu)*(1-2*nu)))
mu = UserExpression(E/2/(1+nu))
rho = 2700.     # density
kappa  = UserExpression(alpha*(2*mu + 3*lmbda))
cV = Constant(910e-6)*rho # specific heat per unit volume at constant strain
k = Constant(237e-6)  # thermal conductivity

d = 1 # interpolation degree
Vue = VectorElement('CG', mesh.ufl_cell(), d) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), d) # temperature finite element
V = FunctionSpace(mesh, MixedElement([Vue, Vte]))

def inner_boundary(x, on_boundary):
    return near(x[0]**2+x[1]**2, R**2, 1e-3) 
def bottom(x, on_boundary):
    return near(x[1], 0) and on_boundary
def left(x, on_boundary):
    return near(x[0], 0) and on_boundary

bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), bottom)
bc2 = DirichletBC(V.sub(0).sub(0), Constant(0.), left)
bc3 = DirichletBC(V.sub(1), 10, inner_boundary)
bcs = [bc1, bc2, bc3]


# Variational formulation and time discretization
U_ = TestFunction(V)
(u_, T_) = split(U_)
dU = TrialFunction(V)
(du, dT) = split(dU)
Uold = Function(V)
(uold, Told) = split(Uold)

def eps(v):
    return sym(grad(v))
def sigma(v, dT):
    return (lmbda*tr(eps(v)) - kappa*dT)*Identity(2) + 2*mu*eps(v)

dt = Constant(0.)
mech_form = inner(sigma(du, dT), eps(u_))*dx
therm_form = (cV*(dT-Told)/dt*T_ + kappa*T0*tr(eps(du-uold))/dt*T_ + dot(k*grad(dT), grad(T_)))*dx
form = mech_form + therm_form

# Resolution
Nincr = 1
t = np.logspace(1, 4, Nincr+1)
Nx = 100
x = np.linspace(R, L, Nx)
T_res = np.zeros((Nx, Nincr+1))
U = Function(V)
(u, T) = split(U)
for (i, dti) in enumerate(np.diff(t)):
    print("Increment " + str(i+1))
    dt.assign(dti)
    solve(lhs(form) == rhs(form), U, bcs)
    Uold.assign(U)
    T_res[:, i+1] = [U(xi, 0.)[2] for xi in x]

plt.figure()
p = plot(sigma(u, T)[1,1], title="$\sigma_{yy}$ stress near the hole")
plt.xlim((0, L))
plt.ylim((0, L))
plt.colorbar(p)
plt.show()