Hi,
I would like to make a very basic thermal expansion in 3d linear analysis (with tet mesh).
I didn’t found anything in the documentation and either in the forum (simple version).
Can someone help me ?
Thanks.
you just need to add the thermal stress term, similar to what is done here (of course without the whole contact stuff, and couple with a thermal simulation):
I don’t succeed to apply thermal_stress on my object, i tried with :
[…]
dt = GridFunction(fes, "dt")
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
a += SymbolicEnergy(
thermal_stress(dt) * u)
pre = Preconditioner(a, "bddc")
a.Assemble()
[…]
Thanks for your answer.
You need to set your grid function to the desired thermal gradient, e.g. with
dt.Set()
Hi,
I tried to do as you said, but my results are not as expected.
To add context, I would like to make a 3D bending beam with sample uniform thermal.
First of all, I made a code to do the linear analysis without thermal, and it work perfectly as you can see on the first picture.
For the material properties, it is also working, but, when I had the thermal expension I obtain this (with the same force on Z, the same boundary etc…):
Here is my code :
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)
def thermal_stress(dt):
global E, nu
cT = 2e-6
return -E * cT * dt / (1 - 2 * nu)
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u,v = fes.TnT()
gfu = GridFunction(fes)
dt = GridFunction(fes, "dt")
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
a += SymbolicEnergy(thermal_stress(dt) * u)
dt.Set((-5e3 * x, -5e3 * y, -5e3 * z)) # Set dummy temperature gradient
pre = Preconditioner(a, "bddc")
a.Assemble()
Do you have any idea of what I did wrong ?
Thanks a lot for your help.
can you send a minimal example which we can run ?
Yes of course, here is my code :
from netgen.occ import *
from ngsolve import *
# Variable
E = 210000
nu = 0.3
# Dirichlet Boundaries
Face_Dirichletx = 'Face4'
Face_Dirichlety = 'Face4'
Face_Dirichletz = 'Face4'
# Loads parameter
FaceForce = 'Face2'
ForceVal = 200
# Degree of Element
EleDeg = 2
# Mesh the geometrie and get the boundaries --------------------------------------------------
geo = OCCGeometry('Geometrie.stp')
mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
boundaries = list(mesh.GetBoundaries())
# Stress Function ----------------------------------------------------------------------------
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)
def thermal_stress(dt):
global E, nu
cT = 2e-6
return -E * cT * dt / (1 - 2 * nu)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u,v = fes.TnT()
gfu = GridFunction(fes)
dt = GridFunction(fes, "dt")
# Assemble of a ----------------------------------------------------------------------------
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
a += SymbolicEnergy(thermal_stress(dt) * u)
dt.Set((-5e3 * x, -5e3 * y, -5e3 * z)) # I would like to make uniform temperature on all the beam, so I guess it is not like that.
pre = Preconditioner(a, "bddc")
a.Assemble()
# Creation of the force value (applied on a face) -----------------------------------------
Dom = mesh.Boundaries(FaceForce)
area = Integrate(1, Dom)
print("Face Area : ", area)
Force_Value = ForceVal/area
force = CF((0, 0, Force_Value))
f = LinearForm(force * v * ds(FaceForce)).Assemble()
# Run the Solver --------------------------------------------------------------------------
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-6, maxiter=500)
gfu.vec.data = inv * f.vec
# Print Result --------------------------------------------------------------------------
with TaskManager():
EFstress = MatrixValued(H1(mesh,order=EleDeg), symmetric=True)
VMstress = GridFunction(EFstress)
VMstress.Interpolate (Stress(Sym(Grad(gfu))))
SymbolicEnergy requires to call AssembleLinearization, like
a.AssembleLinearization(gfu.vec)
otherwise, one gets the error message:
SymbolicEnergy :: CalcMatrix not implemented
Now I get from your code this result:
Joachim
This solve one problem, but the thermal expansion is still not taked into account.
When I change the cT factor, or the value in “dt.set” the result is still the same.
When you put the source term (your thermal force) into the BilinearForm, you have to update with the residual (instead of gfu.vec.data = inv * f.vec):
gfu.vec.data += inv * (f.vec - a.Apply(gfu.vec))
Joachim
I tried with your update, but i still have strange results.
I delete the force to isolate the thermal effect, and I obtain this :
My new code updated :
from netgen.occ import *
from ngsolve import *
# Variable
E = 210000
nu = 0.3
# Dirichlet Boundaries
Face_Dirichletx = 'Face4'
Face_Dirichlety = 'Face4'
Face_Dirichletz = 'Face4'
# Loads parameter
FaceForce = 'Face2'
ForceVal = 0
# Degree of Element
EleDeg = 2
# Mesh the geometrie and get the boundaries --------------------------------------------------
geo = OCCGeometry('Geometrie.stp')
mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
boundaries = list(mesh.GetBoundaries())
# Stress Function ----------------------------------------------------------------------------
mu = E / 2 / (1+nu)
lam = E * nu / ((1+nu)*(1-2*nu))
def Stress(strain):
return 2*mu*strain + lam*Trace(strain)*Id(3)
def thermal_stress(dt):
global E, nu
cT = 2e-12
return -E * cT * dt / (1 - 2 * nu)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u,v = fes.TnT()
gfu = GridFunction(fes)
dt = GridFunction(fes, "dt")
# Assemble of a ----------------------------------------------------------------------------
with TaskManager():
a = BilinearForm(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(v))).Compile()*dx)
a += SymbolicEnergy(thermal_stress(dt) * u)
dt.Set((5e3, 5e3, 5e3))
pre = Preconditioner(a, "bddc")
a.AssembleLinearization(gfu.vec)
# Creation of the force value (applied on a face) -----------------------------------------
Dom = mesh.Boundaries(FaceForce)
area = Integrate(1, Dom)
print("Face Area : ", area)
Force_Value = ForceVal/area
force = CF((0, 0, 0))
# f = LinearForm(force * v * ds(FaceForce)).Assemble()
# Run the Solver --------------------------------------------------------------------------
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-6, maxiter=500)
gfu.vec.data += inv * (a.Apply(dt.vec))
# Print Result --------------------------------------------------------------------------
with TaskManager():
EFstress = MatrixValued(H1(mesh,order=EleDeg), symmetric=True)
VMstress = GridFunction(EFstress)
VMstress.Interpolate (Stress(Sym(Grad(gfu))))
Hi,
I’ve tried many other things, but i didn’t found anything to make it work. I think it’s about the bilinear formulation. Can someone have a solution ?
My latest test code :
from netgen.occ import *
from ngsolve import *
# Variable
E = 210000
nu = 0.3
# Dirichlet Boundaries
Face_Dirichletx = 'Face4'
Face_Dirichlety = 'Face4'
Face_Dirichletz = 'Face4'
# Loads parameter
FaceForce = 'Face2'
ForceVal = 200
# Degree of Element
EleDeg = 2
# Mesh the geometrie and get the boundaries --------------------------------------------------
geo = OCCGeometry('Geometrie.stp')
mesh = Mesh(geo.GenerateMesh(maxh=6)).Curve(3)
boundaries = list(mesh.GetBoundaries())
# Stress Function ----------------------------------------------------------------------------
mu = E / 2 / (1 + nu)
lam = E * nu / ((1 + nu) * (1 - 2 * nu))
def Stress(strain):
return 2 * mu * strain + lam * Trace(strain) * Id(3)
def thermal_stress(gfu):
global E, nu
cT = 2e-6
return -E * cT * gfu / (1 - 2 * nu)
# FE Space ---------------------------------------------------------------------------------
fes = VectorH1(mesh, order=EleDeg, dirichletx=Face_Dirichletx, dirichlety=Face_Dirichlety, dirichletz=Face_Dirichletz)
u, v = fes.TnT()
gfu = GridFunction(fes)
# Assemble of a ----------------------------------------------------------------------------
with TaskManager():
a = BilinearForm(fes)
a += SymbolicEnergy(InnerProduct(Stress(Sym(Grad(u))), Sym(Grad(u))))
a += SymbolicEnergy(thermal_stress(gfu) * u)
pre = Preconditioner(a, "bddc")
a.AssembleLinearization(gfu.vec)
# Creation of the force value (applied on a face) -----------------------------------------
Dom = mesh.Boundaries(FaceForce)
area = Integrate(1, Dom)
print("Face Area : ", area)
#Force_Value = ForceVal / area
#force = CF((0, 0, Force_Value))
# f = LinearForm(force * v * ds(FaceForce)).Assemble()
gfu.Set((10,10,10))
# Run the Solver --------------------------------------------------------------------------
from ngsolve.krylovspace import CGSolver
inv = CGSolver(a.mat, pre, printrates='\r', tol=1e-6, maxiter=500)
gfu.vec.data += inv * (a.Apply(gfu.vec))
# Print Result --------------------------------------------------------------------------
with TaskManager():
EFstress = MatrixValued(H1(mesh, order=EleDeg), symmetric=True)
VMstress = GridFunction(EFstress)
VMstress.Interpolate(Stress(Sym(Grad(gfu))))
Thanks,
I sent you a snippet two weeks ago, it should work.
Can you post a complete (!) example using that, and the result you get form it ?
Yes the snippet that you sent me is integrated in the 2 previous code that I post.
In the latest code that i post, i’ve made some simplification;
(I made only 1 GridFunction (gfu), and I deleted the force apply on face 2)
This is actually my complete code because I open the results with another app (using VTKOutput) :
vtk = VTKOutput(ma=mesh,
coefs=[VMstress,gfu],
names = [“Stress”,“Displacement”],
filename=“C:/desktop”)
vtk.Do()
If you need anything else, please ask.
Thanks a lot for your help.
here is your notebook, with my proposed modification. Isn’t that what you want ?
thermal.ipynb (4.3 KB)
Thanks for what you sent me.
I’ve worked on it, but it still don’t work. There are still things happening that I cannot explain.
First of all, my objective is to apply on a simple object (geometrie.step) a thermal load and see the thermal expension in the linear solver. In order to do this, I use the Duhamel analogy:
X = -D * H*α*ΔT
with respectively:
X the volumic load to apply (acceleration)
D the operator of derivation connecting the displacement to the strain.
α the thermal expension coefficient
H the Hook matrix
End ΔT the temperature difference. (ΔT=T2-T1)
In order to test the program a use a simple beam, but it can be anything else instead.
To show you the result that I would like to obtain, I made a simulation on another FE code.
This calculation was made with the exact same geometrie and material.
With T1=0 and T2=500 K.
The code that you sent me work better but there are still some point to light.
First of all, the result that I obtain with the original value are completly out of range, like 10e4 on X axis. So in order to see what happen, I change the value of ct (α) from 2e-6 to 2e-9 (not physical at all). Then I obtain this:
First weird thing, why there is contraction instead of an expension on the X axis ?
Second thing, the Y and Z displacement look better, but it look like there kind of a temperature grandiant. And I don’t succeed to make it uniform on all the object.
After many test I still don’t understand what happen here.
During these tests I retry with my step beam, and surprise, the results are completly different even if it is the exact same geometrie.
I turn that it is because the geometrie of step file is not properly aligne with the origin.
I don’t know how to solve this problem for the moment.
If you want, I can send you the last code that i’ve done, but there is major change since the last that you send.
Please find here, the geometrie of the beam in step file.
Geometrie.stp (7.8 KB)
Thanks again for your precious help.
Thomas
Hi,
I’m still working on the problem and I’ve maybe found some interesting way.
First of all, the volumic load are implemented by:
dt.set((500*x,500*y,500*z))
This is probably the cause of the problem, because the acceleration is not uniformly apply on the geometry.
In order to apply it uniformly, I think, I need to creat a new coefficient function.
If I understood correctly, I currently use the green function, but to apply the load uniformly, I need to apply the sign function in blue.
I tried this, of course it don’t work but It should do something like that:
dt.Set((500*math.copysign(1,x), 500*math.copysign(1,y),500*math.copysign(1,z)))
Secondly, the origin has to be on the center of mass of the geometry.
So I need to find how to determine the coordinate of the the new origin, in order to use it for the coefficient function. I have no idea, how to do this.
Is there a function to do this automaticcally ?
If you have any idea of solution.
Thanks a lot for your help.
Thomas