Diffusion Additive Manufacturing Work (scratch_cluster_run) Browse

tutorial.py
676f9f89ff91scratch_cluster_run

tutorial.py
View Options

			#!/usr/bin/env python3
			# -- coding: utf-8 --
			"""
			Created on Fri Dec 3 12:35:17 2021

			@author: ekinkubilay
			"""

			#from dolfin import *
			#from dolfin_utils import meshconvert
			from fenics import *
			import os
			import numpy as np

			#define simulation parameters
			alpha = 3
			beta = 1.2
			time = 10
			num_steps = 40
			dt = time/num_steps
			tol = 1e-8
			velocity = 0.36

			#read the mesh
			mesh = Mesh('../Part_geometry/mesh/layer_004.xml')

			#determine the limits of the domain
			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 function space for v
			V = FunctionSpace(mesh, 'P', 1)

			#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)

			#initial source location
			class BoundarySource(SubDomain):
			def inside(self, x, on_boundary):
			return on_boundary and (near(x[2], z_max, tol) and ((x[0]-1.8)2 + x[1]2 < 0.14))
			b_source = BoundarySource()
			b_source.mark(boundary_markers, 3)

			#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 = Expression('10(-4+(x[1]x[1]-2t)+alphax[0]*x[0])', degree=2, alpha=alpha, t=0)
			g_top = Constant(0.0)
			g_sides = Constant(0)
			source = Constant(-500)

			#list of boundary conditions for convenience
			boundary_conditions = {0: {'Dirichlet': T_d},
			1: {'Neumann': g_top},
			2: {'Neumann': g_sides},
			3: {'Neumann': source}}


			#project the initial condition
			T_0 = project(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(gvds(i))



			#write the weak form and seperate into right- and left-hand sides
			F = Tvdx + dtdot(grad(T), grad(v))dx - T_0vdx - dtfvdx + dtsum(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 = File('results/heat_equation_tutorial.pvd')

			#solve the weak form for each timestep
			for n in range(num_steps):

			#increment time and recalculate time dependent functions
			t += dt
			T_d.t = t
			g_top.t = t
			g_sides.t = t

			#solve the weak form
			solve(a == L, T, bcs)

			#save the temperature data the the .vtk file
			vtkfile << (T,t)

			#assign T as the initial condition for the next iteration
			T_0.assign(T)

			#remove Source Boundary from markers and RHS
			b_source.mark(boundary_markers, 4)
			L -= dtgv*ds(3)
			#redefine Source Boundary with new position and mark
			class BoundarySource(SubDomain):
			def inside(self, x, on_boundary):
			return on_boundary and (near(x[2], z_max, tol) and ((x[0]-1.8+velocityt)2 + x[1]*2 < 0.14))
			b_source = BoundarySource()
			b_source.mark(boundary_markers, 3)
			#redefine d in therm of the new boundaries
			ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
			#add new Neumann boundary (Source Boundary) to RHS
			L += dtgv*ds(3)

R11910/Thermal_model/drafts/tutorial.py676f9f89ff91scratch_cluster_run

tutorial.pyView Options

R11910/Thermal_model/drafts/tutorial.py
676f9f89ff91scratch_cluster_run

tutorial.py
View Options