Homec4science
Diffusion Additive Manufacturing Work (mechanical_model_development) Browse

R11910/Thermal_model/superposition.py
8f5e03064641mechanical_model_develo…

Last Change
Raw File

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()
c4science · Help