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

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

#=======================================================================
- (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
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