Page MenuHomec4science
Files F83489838

box_visuals.py
No OneTemporary
Actions

File Metadata

Created
Tue, Sep 17, 10:29

box_visuals.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 22 16:08:08 2021
@author: kubilay
"""
from fenics import *
import os
import numpy as np
from helpers import *
import time as timer
import matplotlib.pyplot as plt
from fem_scan_strategy2 import *
#define simulation parameters
rho = 4420 #kg/m^3
c_p = 750 #J/kgK
k = 27 #W/mK
rho = 8000 #kg/m^3
c_p = 500 #J/kgK
k = 19.2 #W/mK
dt = 0.0005
store = fem_scan_strategy(dt, 60)
time = store[-1,-1]
num_steps = len(store)
print('Simulation will be {} steps'.format(num_steps))
alpha = k/(rho*c_p) #m^2/s
#time = 0.0075 #s
#t_off = 0.002984 #s
power = 175 #W
tol = 1e-14
velocity = 0.8*np.array([1, 0, 0]) #m/s
velocity_mag = np.linalg.norm(velocity) #m/s
source_loc = np.array([0, -0.225, 1.8])*1e-3-np.array([0, 0, 1])*1e-9 #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
mesh = Mesh()
filename = '../../Part_geometry/mesh/layer_060.xdmf'
f = XDMFFile(MPI.comm_world, filename)
f.read(mesh)
MeshTransformation.rescale(mesh, 1e-3, Point(0,0,0)) #in m
#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)
#redefine d in therm of the boundary markers
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
#define Dirichlet initial condition
T_d = Constant(0)
#define Neumann boundary conditions
g_top = Constant(0.0)
g_sides = Constant(0)
#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 = 0
#create directory to save files
try:
os.mkdir('results')
except FileExistsError:
for file in os.scandir('results'):
pass
#define vtk file to save the data at each iteration
vtkfile_fem = File('results/fem_temperature.pvd')
xdmf_file = XDMFFile(MPI.comm_world, "results/fem_temperature.xdmf")
xdmf_file.parameters["flush_output"] = True
xdmf_file.parameters["functions_share_mesh"] = True
xdmf_file.parameters["rewrite_function_mesh"] = False
#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, 1.8])*1e-3-np.array([0, 0, 300])*1e-6
mid_point = np.array([1.19375, -0.225, 1.8])*1e-3-np.array([0, 0, 300])*1e-6
end_point = np.array([2.3875, -0.225, 1.8])*1e-3-np.array([0, 0, 300])*1e-6
midpoint = np.array([-1.19375, -0.225, 1.8])*1e-3
endpoint = np.array([-2.3875, -0.225, 1.8])*1e-3
T_store = np.zeros([num_steps, 10])
points = np.array([midpoint-np.array([0, 0, 100])*1e-6, midpoint-np.array([0, 0, 2*100])*1e-6, midpoint-np.array([0, 0, 3*100])*1e-6,
midpoint-np.array([0, 0, 4*100])*1e-6, midpoint-np.array([0, 0, 5*100])*1e-6, endpoint-np.array([0, 0, 50])*1e-6,
endpoint-np.array([0, 0, 100])*1e-6, endpoint-np.array([0, 0, 150])*1e-6,
endpoint-np.array([0, 0, 200])*1e-6, endpoint-np.array([0, 0, 500])*1e-6])
#calculate time
#start = timer.time()
A, b = assemble_system(a, L, bcs)
#solve the weak form for each timestep
for n in range(num_steps):
#assemble the equations
b = assemble(L)
#determine the location of the source in given time
source = source_loc + velocity*t
#apply the source if it is on
delta = PointSource(V, Point(store[n][0], store[n][1], z_max-tol), dt*power/(rho*c_p))
delta.apply(b)
#solve the system
solve(A, T.vector(), b, "cg", "sor")
#assign the temperature as the initial condition for the next step
T_0.assign(T)
#store the temperature at points of interests
for i,p in enumerate(points):
T_store[n,i] = T(p)
#write temperature data
T.rename("Temperature", "FEM Temperature")
vtkfile_fem << (T,t)
xdmf_file.write(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)))
#output total simulation time
#print('it took:', timer.time()-start)
#print temperature evalution at points of interest
# 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.show()
#store the temperature at points of interests
for i,p in enumerate(points[1:5]):
plt.plot(np.linspace(0,time,num_steps), T_store[:,i+2], marker='.')
plt.show()
#store the temperature at points of interests
plt.figure(figsize=(12, 6))
for i,p in enumerate(points[5:]):
plt.plot(np.linspace(0,time,num_steps), T_store[:,i+5], marker='.')
plt.ylim([-50, 1650])
plt.xlim([-0.01, 0.35])
plt.xlabel('Time (s)')
plt.ylabel('$\Delta$T (K)')
plt.show()
#store the temperature at points of interests
plt.figure(figsize=(12, 5))
plt.rcParams['font.size'] = '16'
plt.plot(np.linspace(0,time,num_steps), T_store[:,0], label='-100', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,1], label='-200', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,2], label='-300', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,3], label='-400', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,4], label='-500', linewidth=4)
plt.ylim([-50, 1650])
plt.xlim([-0.01, 0.3])
plt.xlabel('Time (s)')
plt.ylabel('$\Delta$T (K)')
plt.title('Point A (-1.19375, -0.225) mm')
plt.legend(title='$\Delta z (\mu$m)')
plt.savefig('point_A.png')
plt.show()
#store the temperature at points of interests
plt.figure(figsize=(12, 5))
plt.rcParams['font.size'] = '16'
plt.plot(np.linspace(0,time,num_steps), T_store[:,5], label='-50', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,6], label='-100', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,7], label='-150', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,8], label='-200', linewidth=4)
plt.plot(np.linspace(0,time,num_steps), T_store[:,9], label='-500', linewidth=4)
plt.ylim([-50, 1650])
plt.xlim([-0.01, 0.25])
plt.xlabel('Time (s)')
plt.ylabel('$\Delta$T (K)')
plt.title('Point B (-2.3875, -0.225) mm')
plt.legend(title='$\Delta z (\mu$m)')
plt.savefig('point_B.png')
plt.show()
"""
np.savetxt('Figures/fem_start.txt', T_start/power)
np.savetxt('Figures/fem_middle.txt', T_middle/power)
np.savetxt('Figures/fem_end.txt', T_end/power)
np.savetxt('Figures/fem_time.txt', np.linspace(0,t,num_steps))
"""
"""
def run_analytic_analysis(t):
analytical_temperature = MyExpression(time=t, velocity=velocity, source_loc=source_loc+velocity*t*1e3, alpha=alpha)
T = interpolate(analytical_temperature, V)
#print(np.argmax(T.vector()[:]), np.max(T.vector()[:]))
#print(T(0,0,0.5))
T_start[n] = T(0*5, 2.45*5, 0.2*5)
T_middle[n] = T(1.19375*5, 2.45*5, 0.2*5)
T_end[n] = T(2.3875*5, 2.45*5, 0.2*5)
T.rename("Temperature", "Analytic Temperature")
vtkfile_analytical << (T,t)
def run_fem_analysis(t, T):
source= np.array([0, 0, 0.5]) + velocity*t*1e3
if t < 0.001:
delta = PointSource(V, Point(source), 175)
else:
delta = PointSource(V, Point(source), 0)
delta.apply(b)
solve(A, T.vector(), b)
T_start[n] = T(0, 0, 0.34)
T_middle[n] = T(0.5, 0, 0.34)
T_end[n] = T(0.9, 0, 0.34)
T.rename("Temperature", "FEM Temperature")
vtkfile_fem << (T,t)
return T
"""

Event Timeline

c4science · Help