How to calculate a particular subdomain in thermo-elasticity and save both temperature and displacement results

General
1

Hi, everybody! I’m doing a thermos-elasticity analysis of an object made up of concentric circles. The inner one has a subdomain marker of 1 and the outer on has a subdomain marker of 2. I have so far completed the heat transfer calculation of the whole object. Now I want to calculate the thermal deformation and thermal stress of subdomain 2. I use the “interpolate” function and the code goes like:

from dolfin import *

from subprocess import call

import os

import numpy as np

import vtk

from vtk.util.numpy_support import vtk_to_numpy

import math

#################### Thermal of whole body #####################
Lmbda1 = 0.2
Lmbda2 = 1.2

call("dolfin-convert whole.msh whole.xml;gzip whole *.xml",shell=True)

mesh = Mesh("whole.xml.gz")

boundaries = MeshFunction("size_t", mesh, "whole_facet_region.xml.gz")

subdomains = MeshFunction("size_t", mesh, "whole_physical_region.xml.gz")


VT = FunctionSpace(mesh, "Lagrange", 1)

dT = TrialFunction(VT)

T_ = TestFunction(VT)


bc = DirichletBC(VT, wall_temp , boundaries, wall_id)

integrals_N = g*T_*ds(heat_wall_id)


ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)


aT = (inner(Lmbda1*grad(dT), grad(T_))*dx(1)+inner(Lmbda2*grad(dT), grad(T_))*dx(2)

LT = integrals_N

Delta_T = Function(VT, name="Temperature")

Delta_T.set_allow_extrapolation(True)


solve(aT == LT, Delta_T , bc, solver_parameters={'linear_solver': 'bicgstab','preconditioner': 'hypre_amg'})

#################### elasticity of subdomain 2#####################

E = 11.5e3

nu = 1.5

alpha = 8e-6


call("dolfin-convert subdomain2.msh subdomain2.xml;gzip subdomain2 *.xml",shell=True)

subdomain_mesh = Mesh("subdomain2.xml.gz")

subdomain_boundaries = MeshFunction("size_t", subdomain_mesh, "subdomain2_facet_region.xml.gz")

subdomain_subdomains = MeshFunction("size_t", subdomain_mesh, "subdomain2_physical_region.xml.gz")

subdomain_VT = FunctionSpace(subdomain_mesh, "Lagrange", 1)


subdomain_Delta_T = Function(subdomain_VT,name = " subdomain_temperature")

subdomain_Delta_T.set_allow_extrapolation(True)

subdomain_Delta_T.interpolate(Delta_T)


def m_mu(m_E,m_nu):

return(m_E/2/(1+m_nu))


def m_lmbda(m_E,m_nu):

return(m_E*m_nu/(1+m_nu)/(1-2*m_nu))


def eps(v):

return 0.5*(grad(v) + grad(v).T)


def sigma(v, dT,m_E,m_nu,m_alpha):

return (m_lmbda(m_E,m_nu)*tr(eps(v))- m_alpha*(3*m_lmbda(m_E,m_nu)+2*m_mu(m_E,m_nu))*dT)*Identity(3) + 2.0*m_mu(m_E,m_nu)*eps(v)


Vu = VectorFunctionSpace(subdomain_mesh, 'CG', 1)

du = TrialFunction(Vu)

u_ = TestFunction(Vu)


bc = DirichletBC(Vu.sub(2), Constant(0.0), boundaries, bottom_facet_id)


ds = Measure('ds', domain= subdomain_mesh, subdomain_data= subdomain_boundaries)

dx = Measure('dx', domain= subdomain_mesh, subdomain_data= subdomain_subdomains)


Integrals_Wint= inner(sigma(du, subdomain_Delta_T, E, nu, alpha), eps(u_))*dx

aM = lhs(Integrals_Wint)

LM = rhs(Integrals_Wint)

u = Function(Vu, name="Displacement")


solve(aM == LM, u, bc, solver_parameters={'linear_solver': 'mumps'})


pvd_file_u = File("Subdomain_displacement.pvd")

pvd_file_u << u

pvd_file_T = File("Subdomain_temperature.pvd")

pvd_file_T << subdomain_Delta_T
  1. Is that right to get the displacement of subdomain?
  2. How to output the result of temperature of subdomain? Cuz here encounters error of the last line, which is:
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [12], in <cell line: 143>()
    142 pvd_file_T = File("Subdomain_temperature.pvd")
--> 143 pvd_file_T << subdomain_Delta_T

File ~/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/io/__init__.py:21, in __lshift__(self, u)
     19         self.write(u[0], u[1])
     20 else:
---> 21     self.write(u)

AttributeError: 'Sum' object has no attribute '_cpp_object'