Homec4science
Diffusion Additive Manufacturing Work (scratch_cluster_run) Browse

R11910/Thermal_model/single_scan/analytical_flux_final.py
bd0232b0e72dscratch_cluster_run

Last Change
Raw File

analytical_flux_final.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jan 29 19:20:47 2022
@author: kubilay
"""
from fenics import *
import os
import numpy as np
from helpers import *
import time as timer
import matplotlib.pyplot as plt
from functools import partial
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("time", type=float, help='total simulation time in s')
parser.add_argument("t_off", type=float, help='time when the scan is turned off in s')
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("num_steps", type=int, help='total number of steps')
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
time = args.time #s
t_off = args.t_off #s
power = args.power #W
velocity_mag = args.vel_mag
num_steps = args.num_steps
alpha = k/(rho*c_p) #m^2/s
dt = time/num_steps #s
tol = 1e-14
velocity = velocity_mag*np.array([1, 0, 0]) #m/s
source_loc = np.array([0, -0.225, 2])*1e-3-np.array([0, 0, 1])*1e-6 #m
#read the mesh
mesh = Mesh('../Part_geometry/mesh/layer_005.xml') #in mm
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])
#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)
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_off = Function(V)
t = 0.00001
#create directory to save files
try:
os.mkdir('flux')
except FileExistsError:
for file in os.scandir('flux'):
os.remove(file.path)
#define vtk file to save the data at each iteration
vtkfile_flux = File('flux/flux_temperature.pvd')
start = timer.time()
#T_final = Function(V)
#define arrays to store data on the points of interest
T_start = np.zeros(num_steps)
T_middle = np.zeros(num_steps)
T_end = np.zeros(num_steps)
start_point = np.array([0, -0.225, 2])*1e-3-np.array([0, 0, 300])*1e-6
mid_point = np.array([1.19375, -0.225, 2])*1e-3-np.array([0, 0, 300])*1e-6
end_point = np.array([2.3875, -0.225, 2])*1e-3-np.array([0, 0, 300])*1e-6
#calculate time
start = timer.time()
flux_vector = interpolate(Constant(0.0), V_b)
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)
factor = power/(8*rho*c_p*alpha*(np.pi*alpha)**1.5)
layer = Layer(0.008)
scan = Scan(tol, t_off, source_loc, velocity)
#layer.addScan([scan, Scan(0.0031, 0.006084, source_loc + t_off*velocity + np.array([0,0.00007, 0]), -velocity)])
layer.addScan([scan])
#calculate anaytical temperature at each time step
for n in range(num_steps):
g_sides.vector()[:] = run_correction_fields(boundary_mesh.coordinates()[dof_to_vertex_map(V_b)], normal_vectors, layer, t, alpha, factor)
b = assemble(L)
#solve the system
solve(A, T.vector(), b, "cg", "jacobi")
T_0.assign(T)
#store the temperature at points of interests
T_start[n] = T(start_point)
T_middle[n] = T(mid_point)
T_end[n] = T(end_point)
#write temperature data
T.rename("Temperature", "Flux Temperature")
vtkfile_flux << (T,t)
#increment time and recalculate time dependent functions
t += dt
T_d.t = t
g_top.t = t
#output progress
print('{}% Complete'.format(round(100*n/num_steps,2)))
print('it took:', timer.time()-start)
plt.plot(np.linspace(0,time,num_steps), T_start/power, label='start point', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_middle/power, label='middle point', marker='.')
plt.plot(np.linspace(0,time,num_steps), T_end/power, label='end point', marker='.')
plt.xlabel(r'$time (s)$')
plt.ylabel("$\Delta T/P (K/W)$")
plt.legend()
plt.savefig("Figures/flux_temperature_evolution.jpg")
np.savetxt('Figures/flux_start.txt', T_start/power)
np.savetxt('Figures/flux_middle.txt', T_middle/power)
np.savetxt('Figures/flux_end.txt', T_end/power)
np.savetxt('Figures/flux_time.txt', np.linspace(0,t,num_steps))
if __name__=='__main__':
main()
c4science · Help