Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F96117434
fem_solution
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, Dec 22, 19:36
Size
6 KB
Mime Type
text/x-python
Expires
Tue, Dec 24, 19:36 (2 d)
Engine
blob
Format
Raw Data
Handle
23124589
Attached To
R11910 Additive Manufacturing Work
fem_solution
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
#define simulation parameters
rho = 4420 #kg/m^3
c_p = 750 #J/kgK
k = 27 #W/mK
alpha = k/(rho*c_p) #m^2/s
time = 0.009 #s
num_steps = 120
dt = time/num_steps #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, 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
#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 initial condition
T_d = Constant(0)
T_sides = Constant(0)
#define Neumann boundary conditions
g_top = Constant(0.0)
#list of boundary conditions for convenience
boundary_conditions = {0: {'Dirichlet': T_d},
1: {'Neumann': g_top},
2: {'Dirichlet': T_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]:
if boundary_conditions[i]['Neumann'] != 0:
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'):
os.remove(file.path)
#define vtk file to save the data at each iteration
vtkfile_analytical = File('results/analytic_temperature.pvd')
vtkfile_fem = File('results/fem_temperature.pvd')
#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()
#solve the weak form for each timestep
for n in range(num_steps):
#assemble the equations
A, b = assemble_system(a, L)
#determine the location of the source in given time
source = source_loc + velocity*t
#apply the source if it is on
if t < 0.002984:
delta = PointSource(V, Point(source), 5.27e-5*dt)
delta.apply(b)
#solve the system
solve(A, T.vector(), b, "cg", "jacobi")
#assign the temperature as the initial condition for the next step
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", "FEM Temperature")
vtkfile_fem << (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()
"""
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
Log In to Comment