Splitting Mixed Function Space to Individual Function Spaces- Thermoelasticity

variational formulation
1

Hello to All,

I wanted to separate the function space within Thermo-elastic evolution problem (full coupling) https://comet-fenics.readthedocs.io/en/latest/demo/thermoelasticity/thermoelasticity_transient.html into separate Function Space for a staggered algorithm.

The following is the modified code of the thermoelasticity_transient.py found for the above problem

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

from __future__ import print_function
from dolfin import *
from mshr import *
import numpy as np

L = 1.
R = 0.1
N = 50

domain = Rectangle(Point(0.,0.), Point(L, L)) - Circle(Point(0., 0.), R, 100)
mesh = generate_mesh(domain, N)

T0 = Constant(293.)  
DThole = Constant(10.)
E = 70e3
nu = 0.3
lmbda  = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu)) 
rho = 2700. 
alpha = 2.31e-5 
kappa  = Constant(alpha*(2*mu + 3*lmbda)) 
cV = Constant(910e-6)*rho 
k = Constant(237e-6)  

d = 1 # interpolation degree
Vue = VectorElement('CG', mesh.ufl_cell(), d)
Vte = FiniteElement('CG', mesh.ufl_cell(), d) 
MixedSpace = FunctionSpace(mesh, MixedElement([Vue, Vte]))
VectorSpace = FunctionSpace(mesh,Vue)
ScalarSpace = FunctionSpace(mesh,Vte)

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

bc1 = DirichletBC(VectorSpace.sub(0), Constant(0.), left)
bc2 = DirichletBC(VectorSpace.sub(1), Constant(0.), bottom) 
bcT = DirichletBC(ScalarSpace, DThole, inner_boundary) 
bcU = [bc1, bc2]

u =  Function(VectorSpace)
T =  Function(ScalarSpace)

vu =  TestFunction(VectorSpace)
vT =  TestFunction(ScalarSpace)

un = Function(VectorSpace)
Tn = Function(ScalarSpace)

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

dt = Constant(0.)

u =  TrialFunction(VectorSpace)
mech_form = (inner(sigma(u, Tn), eps(vu)))*dx
au, Lu = lhs(mech_form), rhs(mech_form)
u = Function(VectorSpace)
mechproblem = LinearVariationalProblem(au, Lu, u, bcU)
mechsolver = LinearVariationalSolver(mechproblem)

T =  TrialFunction(ScalarSpace)
therm_form = (cV*(T-Tn)/dt*vT + kappa*T0*tr(eps(u-un))/dt*vT + dot(k*grad(T), grad(vT)))*dx
aT, LT = lhs(therm_form), rhs(therm_form)
T = Function(ScalarSpace)
thermalproblem = LinearVariationalProblem(aT, LT, T, bcT)
thermsolver = LinearVariationalSolver(thermalproblem)

fileResults = XDMFFile('FENICS_FILE/Dynamic.xdmf')
fileResults.parameters["flush_output"] = True
fileResults.parameters["functions_share_mesh"] = True
tt = 0

Nincr = 10
t = np.logspace(1, 4, Nincr+1)
Nx = 100

u.rename("Displacement", "label")
T.rename("Temperature", "label")
fileResults.write(u,tt)
fileResults.write(T,tt)

for (i, dti) in enumerate(np.diff(t)):
     dt.assign(dti)
    print("Increment " + str(i+1),float(dt))
    tt+=float(dt)

    thermsolver.solve()
    Tn.assign(T)

    mechsolver.solve()
    un.assign(u)
    
    fileResults.write(u,tt)
    fileResults.write(T,tt)

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

Why is this code giving a different result to the original code?

2

Hi,
I haven’t tried to implement the staggered approach for the same problem but I imagine that the difference comes from the fact that you use very large time increments and only one pass of staggered iterations
The monolithic approach finds the solution for a given time step in a single solve but the staggered approach should be iterated until convergence on this time step to yield the same solution.
However, both approaches should yield the same solution for very small time steps.

2 Likes
3

Thank you Sir, The reason seems apt.

This staggering was done in order to implement the FAR91 Paper, which eventually led to the following question: