Thermoelastic with NewtonMinimization

NGSolve
1

Hi there,

I am trying to solve a thermoelastic problem using NewtonMinimization, but I do not get any results. I have two materials (part1 and part2) with two different temperature gradients (dT1 and dT2). I am using NewtonMinimization to obtain the thermal expansions. Can someone please help me with this? Below is my code:

heat2.ipynb (3.7 KB)

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.csg import *

def thermal_stress(delta_T, E, nu):
cT = 2e-6
return -E * cT * delta_T / (1 - 2 * nu)

def stress(strain, mu, lam):
return 2 * mu * (strain) + lam * Trace(strain) * Id(3)

def mu(E, nu):
return E / 2 / (1 + nu)

def lam(E, nu):
return E * nu / ((1 + nu) * (1 - 2 * nu))

box1 = Box((0,0,0), (10,2.5,0.7))
box1.mat(‘part1’)
box1.faces.Min(Y).maxh = 0.1
box2 = Box((0,0,0), (10,-1,3))
box2.mat(‘part2’)
box2.faces.col = (0,0,1)

geo = OCCGeometry(Glue([box1,box2]))

mesh = Mesh(geo.GenerateMesh(maxh=0.5))

Draw(mesh)

Mechanical parameters:

E_p1, nu_p1 = 2.10, 0.1
E_p2, nu_p2 = 80, 0.1

fes = VectorH1(mesh, order=3)

u, v = fes.TnT()
gfu = GridFunction(fes)
dT1 = GridFunction(fes, definedon=mesh.Materials(“part1”))
dT2 = GridFunction(fes, definedon=mesh.Materials(“part2”))

a = BilinearForm(fes)
a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p1, nu_p1), lam=lam(E_p1, nu_p1)),
Sym(Grad(u))
), definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p2, nu_p2), lam=lam(E_p2, nu_p2)),
Sym(Grad(u))
), definedon=mesh.Materials(“part2”)
)

a += SymbolicEnergy(
thermal_stress(dT1, E_p1, nu_p1) *
u, definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
thermal_stress(dT2, E_p2, nu_p2) *
u, definedon=mesh.Materials(“part2”)
)

dT1.Set((-5e3x, -5e3y, -5e3z))
dT2.Set((-1e3x, -1e3y, -1e3z))

print(‘start’)

solvers.NewtonMinimization(a, gfu, printing=False, inverse=“sparsecholesky”)

print(‘end’)

2

I changed the mesh size, and NewtonMinimization works now but does not converge. Can anyone help me to understand what I am doing wrong, please? Below is the updated code:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.csg import *

def thermal_stress(delta_T, E, nu):
cT = 2e-6
return -E * cT * delta_T / (1 - 2 * nu)

def stress(strain, mu, lam):
return 2 * mu * (strain) + lam * Trace(strain) * Id(3)

def mu(E, nu):
return E / 2 / (1 + nu)

def lam(E, nu):
return E * nu / ((1 + nu) * (1 - 2 * nu))

box1 = Box((0,0,0), (10,2.5,0.7))
box1.mat(‘part1’)
#box1.faces.Min(Y).maxh = 0.1
box2 = Box((0,0,0), (10,-1,3))
box2.mat(‘part2’)
box2.faces.col = (0,0,1)

geo = OCCGeometry(Glue([box1,box2]))

mesh = Mesh(geo.GenerateMesh(maxh=0.5))

Draw(mesh)

Dumy mechanical parameters:

E_p1, nu_p1 = 56, 0.3
E_p2, nu_p2 = 80, 0.2

fes = VectorH1(mesh, order=3)

u, v = fes.TnT()
gfu = GridFunction(fes)
dT1 = GridFunction(fes, definedon=mesh.Materials(“part1”))
dT2 = GridFunction(fes, definedon=mesh.Materials(“part2”))

a = BilinearForm(fes)
a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p1, nu_p1), lam=lam(E_p1, nu_p1)),
Sym(Grad(u))
), definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p2, nu_p2), lam=lam(E_p2, nu_p2)),
Sym(Grad(u))
), definedon=mesh.Materials(“part2”)
)

a += SymbolicEnergy(
thermal_stress(dT1, E_p1, nu_p1) *
u, definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
thermal_stress(dT2, E_p2, nu_p2) *
u, definedon=mesh.Materials(“part2”)
)

Dumy temperature gradient

dT1.Set((-5e3x, -5e3y, -5e3z))
dT2.Set((-5e3x, -5e3y, -5e3z))

solvers.NewtonMinimization(a, gfu, printing=True, maxit=200, maxerr=1e-1, inverse=“sparsecholesky”)

3

I added a contact boundary, but still does not converge. Below is my code:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
from netgen.csg import *
from ngsolve.solvers import *

def thermal_stress(delta_T, E, nu):
cT = 10e-6
return -E * cT * delta_T / (1 - 2 * nu)

def stress(strain, mu, lam):
return 2 * mu * (strain) + lam * Trace(strain) * Id(3)

def mu(E, nu):
return E / 2 / (1 + nu)

def lam(E, nu):
return E * nu / ((1 + nu) * (1 - 2 * nu))

box1 = Box((0,0,0), (10,2.5,0.7))
box1.mat(‘part1’)
box1.faces.Min(Y).name = ‘buttom_contact’
#box1.faces.Min(Y).maxh = 0.1
box2 = Box((0,0,0), (10,-1,3))
box2.mat(‘part2’)
box2.faces.Max(Y).name = ‘top_contact’
box2.faces.col = (0,0,1)

geo = OCCGeometry(Glue([box1,box2]))

mesh = Mesh(geo.GenerateMesh(maxh=0.5))

Draw(mesh)

Dumy mechanical parameters:

E_p1, nu_p1 = 56, 0.3
E_p2, nu_p2 = 80, 0.2

fes = VectorH1(mesh, order=3)

u, v = fes.TnT()
gfu = GridFunction(fes)
dT1 = GridFunction(fes)
#dT2 = GridFunction(fes, definedon=mesh.Materials(“part2”))
dT1.Set((-5e3x, -5e3y, -5e3*z))

a = BilinearForm(fes)
a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p1, nu_p1), lam=lam(E_p1, nu_p1)),
Sym(Grad(u))
), definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
0.5 * InnerProduct(
stress(Sym(Grad(u)), mu=mu(E_p2, nu_p2), lam=lam(E_p2, nu_p2)),
Sym(Grad(u))
), definedon=mesh.Materials(“part2”)
)

a += SymbolicEnergy(
thermal_stress(dT1, E_p1, nu_p1) *
u, definedon=mesh.Materials(“part1”)
)

a += SymbolicEnergy(
thermal_stress(dT1, E_p2, nu_p2) *
u, definedon=mesh.Materials(“part2”)
)

contact = ContactBoundary(
mesh.Boundaries(“buttom_contact”),
mesh.Boundaries(“top_contact”)
)

X = CoefficientFunction((x, y, z))
cf = (X + u - (X.Other() + u.Other())) * contact.normal
contact.AddEnergy(IfPos(cf, 1e14cfcf, 0))

Dumy temperature gradient

with TaskManager():

contact.Update(gfu, a)
NewtonMinimization(a, gfu, printing=True, inverse="sparsecholesky")
4

I found out what was the problem. I did not define boundary conditions and the model was unstable. After defining the boundary condition it works now. Thanks