How to change sub-domains and their property by "MeshFunction" marker?


109
views
0
4 months ago by
Hi,
I am doing 3D simulation of heat equation in a cube, where in the beginning at the bottom of the cube there is solid powder and the the rest is air. After some time steps I need to add another powder layer. The materials property is depend on the phases(solid,liquid,powder) So I decided to use the linear solver and change materials properties after each solve by changing the sub_domains.
The main code and the Function are bellow.
It works fine. But the problem is:
If instead of adding powder before defining the equation(F), it is added after that(commenting the line 41 and un-comment line 67) result shows that powder dose not added at all.
So my question is: is there any way to have linear solver and change material properties at different time steps?

many thanks in advance.
from fenics import *
from Functions import Materials, Laser,Boundary, Refine
#=======================================================================
laser_power = 0.2 * 400	    		# [j/(s)] = [W]  
scan_speed = 0.5                	# [m/s]
laser_radius = 400 * 1e-06		    # [m]
lw =0.5 * 1e-03
powderlayer = 1

#=======================================================================
x0 = y0 = z0 = -2.5 * 1e-03		# first point in each direction
x1 = y1 = z1 = +2.5 * 1e-03		# last point in each direction

nx = ny = nz = 10
deltax = (x1-x0) / nx
deltay = (y1-y0) / ny
deltaz = (z1-z0) / nz
dt = deltax / (scan_speed)
powder_height = z0 + powderlayer * lw
laser_coor_x = x0 
laser_coor_y = y0 
laser_coor_z = z0 + powderlayer * lw - 0.5 * deltaz
#=======================================================================
f = Laser(laser_power,laser_radius,laser_coor_x,laser_coor_y,laser_coor_z,powder_height) 	# [w/m3]
#=======================================================================
original_mesh = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), nx, ny, nz)
refine_criteria = deltaz/2
refine_level = 2
mesh = Refine(laser_coor_x,laser_coor_y,laser_coor_z,original_mesh,refine_criteria,refine_level,powder_height)
V = FunctionSpace(mesh, 'CG', 2)
T_0 = Constant(298.15)	 #[K] = 25 c
T_n = interpolate(T_0, V)
T = TrialFunction(V)
v = TestFunction(V)
#=======================================================================
T_D = Constant(1073.15)	#[k] = 800 c
dT_Neumann = Constant(0.0)
bc, dT_dn = Boundary(V,T_D,dT_Neumann,z0)
#=======================================================================
k,ru,cp = Materials(mesh, z0) # without any powder
k,ru,cp = Materials(mesh, powder_height ) # add powder

#=======================================================================
F = T*v*dx + (dt*k/(ru*cp))*dot(grad(T), grad(v))*dx	\
		- (T_n*v*dx) - (dt/(ru*cp))*f*v*dx				\
		- (dt*k/(ru*cp))*dT_dn*v*ds			# due to neumann condition

a, L = lhs(F), rhs(F)
T = Function(V)

problem = LinearVariationalProblem(a, L, T, bc)
solver = LinearVariationalSolver(problem)
#=======================================================================
vtkfile_T = File('ssolution/T1.pvd')

tsteps = 6
t = 0

f_power = f.Volumetric_laser_power
f.Volumetric_laser_power = 0.0	#turn laser off
t += dt
solver.solve()
vtkfile_T << (T, t)
T_n.assign(T)
f.Volumetric_laser_power = f_power #turn laser on

#k,ru,cp = Materials(mesh, powder_height ) # add powder

t += dt
solver.solve()
vtkfile_T << (T, t)
T_n.assign(T)

for n in range(tsteps):
	t += dt 
	solver.solve()
	vtkfile_T << (T, t)
	T_n.assign(T)
	f.laser_coor_x += deltax

​

Materials: Two sub-domains air and powder are defined and marked:


def Materials(mesh,powder_height):
	# define a meshfunction for numbering subdomains
	subdomains = MeshFunction("size_t",mesh,3)

	# define the subdomains	
	class Air(SubDomain):
		def inside(self, x,on_boundary):
			return (x[2] >= powder_height)
	
	class Powder(SubDomain):
		def inside(self, x,on_boundary):
			return (x[2] <= powder_height)
		

	# mark the subdomains with numbers
	subdomain_air = Air()
	subdomain_powder= Powder()
	
	subdomains.set_all(1)
	subdomain_air.mark(subdomains,0)

	V0 = FunctionSpace(mesh, "DG", 0)
	k = Function(V0)
	ru = Function(V0)
	cp = Function(V0)
	
	help = numpy.asarray(subdomains.array(), dtype=numpy.int64)
	k.vector()[:] = numpy.choose(help, k_values)
	ru.vector()[:] = numpy.choose(help, ru_values)
	cp.vector()[:] = numpy.choose(help, cp_values)

	return(k,ru,cp)​

Boundary: There is Dirichlet boundary condition at the botomn and other edge has homogenous Neumann condition:

def Boundary(V,u_Dirichlet,du_dn_Neumann,z0):
	# Define boundaries
	tol=1e-14
	bot   = 'near(x[2], -2.5 * 1e-03, tol)'
	top   = 'near(x[2], +2.5 * 1e-03, tol)' 
	left  = 'near(x[0], -2.5 * 1e-03, tol)' 
	right = 'near(x[0], +2.5 * 1e-03, tol)' 

	front = 'near(x[1], -2.5 * 1e-03, tol)' 
	back  = 'near(x[1], +2.5 * 1e-03, tol)' 

	def boundary_D(x, on_boundary):
		if on_boundary:
			if near(x[2], z0, tol): #bot
				#print('z=%g' % x[2])
				return True
			else:
				return False
		else:
			return False

	#Dirichlet condition
	u_D = u_Dirichlet  
	bc = DirichletBC(V, Constant(u_D), boundary_D)
	#Neuman condition
	du_dn = Constant(du_dn_Neumann)
	return(bc,du_dn)

laser: laser will be applied to the powder

def Laser(laser_power,laser_radius,xf,yf,zf,zp):
	# laser applied to midlle of elements
	# in gaussian distribution 99.73% data are within mu-3sigma & mu+3sigma
	#       ---> laser_radius = 3*sigma
	tol = 1e-7
	stdev = laser_radius / 3.0      # [m] 
	a = -0.5 / (stdev * stdev) 		
	sqrt_det_covariance_matrix = stdev * stdev * stdev 		# [m3]
	Volumetric_laser_power = laser_power / (sqrt_det_covariance_matrix * pow(2*pi,1.5)) # [w/m3]

	f = Expression('(x[2]>=powder_height-tol ? 0:2) * Volumetric_laser_power \
					* exp( a * (										\
						+ (x[0]-laser_coor_x) * (x[0]-laser_coor_x) 	\
						+ (x[1]-laser_coor_y) * (x[1]-laser_coor_y)		\
						+ (x[2]-laser_coor_z) * (x[2]-laser_coor_z)))'	\
						 , degree=2,  a=a, laser_coor_x=xf,				\
						 laser_coor_y=yf, laser_coor_z=zf,				\
						 Volumetric_laser_power=Volumetric_laser_power,	\
						 powder_height=zp, tol=tol)
	return(f)

Refine: mesh is locally refined by CellFunction marker:

def Refine(x,y,z,mesh,refine_criteria,repeat,powder_height):
    #Local refinement
	refined_mesh = mesh
	for n in range(repeat):
		cell_markers = CellFunction("bool", mesh)
		cell_markers.set_all(False)
		origin = Point(x, y, z)
		for cell in cells(mesh):
			center = cell.midpoint()
			if abs(center.z()-z) < refine_criteria:
				cell_markers[cell] = True
		refined_mesh = refine(mesh, cell_markers,False)
		mesh = refined_mesh
		refine_criteria = refine_criteria / 2.0
	if(rank==0):
		print('refine is done....powder_height:%g    x:%g  y:%g  z:%g ' %(powder_height, x, y, z))
	return(refined_mesh)
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »