Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88763313
superposition.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Oct 20, 13:52
Size
8 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 13:52 (2 d)
Engine
blob
Format
Raw Data
Handle
21807939
Attached To
R11910 Additive Manufacturing Work
superposition.py
View Options
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Feb 3 18:53:00 2022
@author: kubilay
"""
from
fenics
import
*
import
os
import
numpy
as
np
from
helpers
import
*
import
time
as
timer
import
matplotlib.pyplot
as
plt
import
argparse
#define simulation parameters
rho
=
4420
#kg/m^3
c_p
=
750
#J/kgK
k
=
27
#W/mK
power
=
175
#W
velocity_mag
=
0.8
dt
=
0.0005
alpha
=
k
/
(
rho
*
c_p
)
#m^2/s
tol
=
1e-14
t
=
0.00001
#read the mesh
mesh
=
Mesh
()
filename
=
'../Part_geometry/mesh/layer_001.xdmf'
f
=
XDMFFile
(
MPI
.
comm_world
,
filename
)
f
.
read
(
mesh
)
MeshTransformation
.
rescale
(
mesh
,
1e-3
,
Point
(
0
,
0
,
0
))
#in m
boundary_mesh
=
BoundaryMesh
(
mesh
,
'exterior'
)
#determine the limits of the domain in m
x_max
=
np
.
max
(
mesh
.
coordinates
()[:,
0
])
x_min
=
np
.
min
(
mesh
.
coordinates
()[:,
0
])
y_max
=
np
.
max
(
mesh
.
coordinates
()[:,
1
])
y_min
=
np
.
min
(
mesh
.
coordinates
()[:,
1
])
z_max
=
np
.
max
(
mesh
.
coordinates
()[:,
2
])
z_min
=
np
.
min
(
mesh
.
coordinates
()[:,
2
])
factor_flux
=
power
/
(
8
*
rho
*
c_p
*
alpha
*
(
np
.
pi
*
alpha
)
**
1.5
)
factor_analytical
=
power
/
(
4
*
rho
*
c_p
*
(
np
.
pi
*
alpha
)
**
1.5
)
infill_velocity_mag
=
velocity_mag
contour_velocity_mag
=
0.9
infill_scan_limits
=
get_interpolation_limits
(
0.00377
/
2
,
infill_velocity_mag
,
alpha
)
contour_scan_limits
=
get_interpolation_limits
(
0.00377
/
2
,
contour_velocity_mag
,
alpha
)
layer
=
get_scan_strategy
(
23
,
contour_velocity_mag
,
infill_velocity_mag
,
contour_scan_limits
,
infill_scan_limits
,
0.01
)
time
=
1.1
*
layer
.
totalLayerTime
()
num_steps
=
int
(
np
.
ceil
(
time
/
dt
))
print
(
'Simulation will have {} steps'
.
format
(
num_steps
))
print
(
'{} seconds will be simulated'
.
format
(
time
))
#define markers for the boundary faces
boundary_markers
=
MeshFunction
(
'size_t'
,
mesh
,
mesh
.
topology
()
.
dim
()
-
1
,
4
)
#define boundary subclasses for bottom, top and side surfaces, mark the boundaries
class
BoundaryBottom
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
near
(
x
[
2
],
z_min
,
tol
)
b_bottom
=
BoundaryBottom
()
b_bottom
.
mark
(
boundary_markers
,
0
)
class
BoundaryTop
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
near
(
x
[
2
],
z_max
,
tol
)
b_top
=
BoundaryTop
()
b_top
.
mark
(
boundary_markers
,
1
)
class
BoundarySides
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
(
near
(
x
[
0
],
x_min
,
tol
)
or
near
(
x
[
0
],
x_max
,
tol
)
or
near
(
x
[
1
],
y_min
,
tol
)
or
near
(
x
[
1
],
y_max
,
tol
))
b_sides
=
BoundarySides
()
b_sides
.
mark
(
boundary_markers
,
2
)
#define function space for v
V
=
FunctionSpace
(
mesh
,
"CG"
,
1
)
V_b
=
FunctionSpace
(
boundary_mesh
,
"CG"
,
1
)
#redefine d in therm of the boundary markers
ds
=
Measure
(
'ds'
,
domain
=
mesh
,
subdomain_data
=
boundary_markers
)
#define initial condition
T_d
=
Constant
(
0
)
#define Neumann boundary conditions
g_top
=
Constant
(
0.0
)
g_sides
=
interpolate
(
Constant
(
0.0
),
V_b
)
#list of boundary conditions for convenience
boundary_conditions
=
{
0
:
{
'Dirichlet'
:
T_d
},
1
:
{
'Neumann'
:
g_top
},
2
:
{
'Neumann'
:
g_sides
}}
#interpolate the initial condition
T_0
=
interpolate
(
Constant
(
0.0
),
V
)
#define test and trial function and the source term
T
=
TrialFunction
(
V
)
T_superposition
=
interpolate
(
Constant
(
0.0
),
V
)
v
=
TestFunction
(
V
)
#f = Constant(0)
#gather all Dirichlet boundary conditions in one list
bcs
=
[]
for
i
in
boundary_conditions
:
if
'Dirichlet'
in
boundary_conditions
[
i
]:
bc
=
DirichletBC
(
V
,
boundary_conditions
[
i
][
'Dirichlet'
],
boundary_markers
,
i
)
bcs
.
append
(
bc
)
#gather all Neumann boundary conditions in weak form
integrals_N
=
[]
for
i
in
boundary_conditions
:
if
'Neumann'
in
boundary_conditions
[
i
]:
g
=
boundary_conditions
[
i
][
'Neumann'
]
integrals_N
.
append
(
g
*
v
*
ds
(
i
))
#write the weak form and seperate into right- and left-hand sides
F
=
T
*
v
*
dx
+
alpha
*
dt
*
dot
(
grad
(
T
),
grad
(
v
))
*
dx
-
T_0
*
v
*
dx
-
alpha
*
dt
*
sum
(
integrals_N
)
a
,
L
=
lhs
(
F
),
rhs
(
F
)
#define the function and time
T
=
Function
(
V
)
T_0
=
interpolate
(
Constant
(
0.0
),
V
)
#create directory to save files
try
:
os
.
mkdir
(
'results/layer_scan'
)
except
FileExistsError
:
pass
#define vtk file to save the data at each iteration
xdmf_file
=
XDMFFile
(
MPI
.
comm_world
,
"results/layer_scan/layer_scan2.xdmf"
)
xdmf_file
.
parameters
[
"flush_output"
]
=
True
xdmf_file
.
parameters
[
"functions_share_mesh"
]
=
True
xdmf_file
.
parameters
[
"rewrite_function_mesh"
]
=
False
start
=
timer
.
time
()
#calculate time
start
=
timer
.
time
()
flux_vector
=
interpolate
(
Constant
(
0.0
),
V_b
)
#define arrays to store data on the points of interest
T_1
=
np
.
zeros
(
num_steps
)
T_2
=
np
.
zeros
(
num_steps
)
T_3
=
np
.
zeros
(
num_steps
)
T_4
=
np
.
zeros
(
num_steps
)
T_5
=
np
.
zeros
(
num_steps
)
p_1
=
np
.
array
([
4.9
/
2
,
4.9
/
2
,
2
-
0.1
])
*
1e-3
p_2
=
np
.
array
([
4.9
/
2
,
4.9
/
2
,
2
-
0.2
])
*
1e-3
p_3
=
np
.
array
([
4.9
/
2
,
4.9
/
2
,
2
-
0.3
])
*
1e-3
p_4
=
np
.
array
([
4.9
/
2
,
4.9
/
2
,
2
-
0.4
])
*
1e-3
p_5
=
np
.
array
([
4.9
/
2
,
4.9
/
2
,
2
-
0.5
])
*
1e-3
boundary_normals
=
get_boundary_normals
(
mesh
)
normal_vectors
=
np
.
zeros
(
np
.
shape
(
boundary_mesh
.
coordinates
()[
dof_to_vertex_map
(
V_b
)]))
for
i
,
v
in
enumerate
(
boundary_mesh
.
coordinates
()[
dof_to_vertex_map
(
V_b
)]):
normal
=
boundary_normals
(
v
)
normal
/=
np
.
max
(
normal
)
normal_vectors
[
i
]
=
normal
/
np
.
linalg
.
norm
(
normal
)
A
,
b
=
assemble_system
(
a
,
L
,
bcs
)
dim_len
=
len
(
T
.
vector
()
.
get_local
())
#calculate anaytical temperature at each time step
for
n
in
range
(
200
):
g_sides
.
vector
()
.
set_local
(
run_correction_fields
(
boundary_mesh
.
coordinates
()[
dof_to_vertex_map
(
V_b
)],
normal_vectors
,
layer
,
t
,
alpha
,
factor_flux
))
as_backend_type
(
g_sides
.
vector
())
.
vec
()
.
ghostUpdate
()
#A, b = assemble_system(a, L, bcs)
b
=
assemble
(
L
)
#solve the system
solve
(
A
,
T
.
vector
(),
b
,
"cg"
,
"jacobi"
)
as_backend_type
(
T
.
vector
())
.
vec
()
.
ghostUpdate
()
#T_superposition.vector().set_local(T.vector()[:] + run_analytical_fields(mesh.coordinates()[dof_to_vertex_map(V)], layer, t, alpha, factor_analytical))
#as_backend_type(T_superposition.vector()).vec().ghostUpdate()
T_superposition
.
vector
()
.
set_local
(
run_analytical_fields
(
mesh
.
coordinates
()[
dof_to_vertex_map
(
V
)],
layer
,
t
,
alpha
,
factor_analytical
))
as_backend_type
(
T_superposition
.
vector
())
.
vec
()
.
ghostUpdate
()
T_superposition
.
vector
()
.
set_local
(
T
.
vector
()[:]
+
T_superposition
.
vector
()[:])
as_backend_type
(
T_superposition
.
vector
())
.
vec
()
.
ghostUpdate
()
#write temperature data
if
n
%
10
==
0
:
#print('step {}'.format(n))
#write temperature data
T_superposition
.
rename
(
"Temperature"
,
"Superposition Temperature"
)
xdmf_file
.
write
(
T_superposition
,
t
)
T_0
.
assign
(
T
)
for
scan
in
layer
:
if
t
>
1.1
*
scan
.
t_end
:
layer_off
=
Layer
()
layer_off
.
addScan
(
scan
)
T_0
.
vector
()
.
set_local
(
T
.
vector
()
.
get_local
()
+
run_analytical_fields
(
mesh
.
coordinates
()[
dof_to_vertex_map
(
V
)[
0
:
dim_len
]],
layer_off
,
t
,
alpha
,
factor_analytical
))
as_backend_type
(
T
.
vector
())
.
vec
()
.
ghostUpdate
()
layer
.
removeScan
(
scan
)
#increment time and recalculate time dependent functions
t
+=
dt
"""
try:
T_1[n] = T(p_1)
except RuntimeError:
pass
try:
T_2[n] = T(p_2)
except RuntimeError:
pass
try:
T_3[n] = T(p_3)
except RuntimeError:
pass
try:
T_4[n] = T(p_4)
except RuntimeError:
pass
try:
T_5[n] = T(p_5)
except RuntimeError:
pass
"""
#output progress
print
(
'{}% Complete'
.
format
(
round
(
100
*
n
/
num_steps
,
2
)))
print
(
'it took:'
,
timer
.
time
()
-
start
)
"""
#print temperature evalution at points of interest
plt.plot(np.linspace(0,time,num_steps), T_1, label='p_1', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_2, label='p_2', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_3, label='p_3', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_4, label='p_4', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_5, label='p_5', marker='.')
#plt.ylim([-20, 2000])
plt.xlabel(r'$time (s)$')
plt.ylabel("$\Delta T (K)$")
plt.legend()
plt.savefig("Figures/layer_scan_temperature_evolution.jpg")
"""
#print temperature evalution at points of interest
plt
.
plot
(
T_1
[
0
:
200
],
label
=
'p_1'
,
marker
=
'.'
)
plt
.
plot
(
T_2
[
0
:
200
],
label
=
'p_2'
,
marker
=
'.'
)
plt
.
plot
(
T_3
[
0
:
200
],
label
=
'p_3'
,
marker
=
'.'
)
plt
.
plot
(
T_4
[
0
:
200
],
label
=
'p_4'
,
marker
=
'.'
)
plt
.
plot
(
T_5
[
0
:
200
],
label
=
'p_5'
,
marker
=
'.'
)
#plt.ylim([-20, 2000])
plt
.
xlabel
(
r'$time (s)$'
)
plt
.
ylabel
(
"$\Delta T (K)$"
)
plt
.
legend
()
plt
.
show
()
Event Timeline
Log In to Comment