Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76523907
layer_scan.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
Thu, Aug 8, 11:10
Size
9 KB
Mime Type
text/x-python
Expires
Sat, Aug 10, 11:10 (2 d)
Engine
blob
Format
Raw Data
Handle
19724435
Attached To
R11910 Additive Manufacturing Work
layer_scan.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
def
parse_args
():
"""
Function to handle arguments
"""
parser
=
argparse
.
ArgumentParser
(
description
=
'Run FEM of single scan on a simple box'
)
parser
.
add_argument
(
"rho"
,
type
=
float
,
help
=
'density in kg/m^3'
)
parser
.
add_argument
(
"c_p"
,
type
=
float
,
help
=
'specific heat in J/kgK'
)
parser
.
add_argument
(
"k"
,
type
=
float
,
help
=
'thermal conductivity in W/Km'
)
parser
.
add_argument
(
"power"
,
type
=
float
,
help
=
'laser power in W'
)
parser
.
add_argument
(
"vel_mag"
,
type
=
float
,
help
=
'velocity magnitude in m/s'
)
parser
.
add_argument
(
"time_step"
,
type
=
float
,
help
=
'time step for simulation (s)'
)
args
=
parser
.
parse_args
()
return
args
def
main
():
#get the arguments
args
=
parse_args
()
#define simulation parameters
rho
=
args
.
rho
#kg/m^3
c_p
=
args
.
c_p
#J/kgK
k
=
args
.
k
#W/mK
power
=
args
.
power
#W
velocity_mag
=
args
.
vel_mag
dt
=
args
.
time_step
alpha
=
k
/
(
rho
*
c_p
)
#m^2/s
tol
=
1e-14
velocity
=
velocity_mag
*
np
.
array
([
1
,
0
,
0
])
#m/s
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
(
197
,
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
)
#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_scan.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
([
0
,
0
,
13
-
0.4
])
*
1e-3
p_2
=
np
.
array
([
0
,
0
,
13
-
0.4
])
*
1e-3
p_3
=
np
.
array
([
0
,
0
,
13
-
0.3
])
*
1e-3
p_4
=
np
.
array
([
0
,
0
,
13
-
0.2
])
*
1e-3
p_5
=
np
.
array
([
0
,
0
,
13
-
0.1
])
*
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
(
num_steps
):
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
,
tensor
=
b
)
#solve the system
solve
(
A
,
T
.
vector
(),
b
,
"cg"
,
"sor"
)
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
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
()[:]
+
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
"""
T_1[n] = T_superposition(p_1)
T_2[n] = T_superposition(p_2)
T_3[n] = T_superposition(p_3)
T_4[n] = T_superposition(p_4)
T_5[n] = T_superposition(p_5)
"""
#output progress
#print('{}% Complete'.format(round(100*n/num_steps,2)))
"""
#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
(
'it took:'
,
timer
.
time
()
-
start
)
if
__name__
==
'__main__'
:
main
()
Event Timeline
Log In to Comment