Skip to content

Instantly share code, notes, and snippets.

@ilyasst
Last active October 21, 2021 21:43
Show Gist options
  • Save ilyasst/7fcafc757e1f045845ecc9e914a3f79e to your computer and use it in GitHub Desktop.
Save ilyasst/7fcafc757e1f045845ecc9e914a3f79e to your computer and use it in GitHub Desktop.
Fenisc fiber inclusion thermomechanical problem: two materials
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment