From fcb9460ea6dc65b3075d0730341f09a3e7832c61 Mon Sep 17 00:00:00 2001 From: Saswata Hier-Majumder Date: Tue, 12 Nov 2019 18:26:32 +0000 Subject: [PATCH] 3 component advection diffusion reaction solver added to darcy advection class and a new script LVL_3comp added --- src/modules/mupopp.py | 112 ++++++++++++++- src/run/LVL_3comp/LVL_3comp.py | 246 +++++++++++++++++++++++++++++++++ 2 files changed, 357 insertions(+), 1 deletion(-) create mode 100644 src/run/LVL_3comp/LVL_3comp.py diff --git a/src/modules/mupopp.py b/src/modules/mupopp.py index 839c055..033112b 100644 --- a/src/modules/mupopp.py +++ b/src/modules/mupopp.py @@ -848,7 +848,117 @@ def advection_diffusion_two_component(self,W,mesh,sol_prev,dt,f1,K=1.0,zh=Consta term_SUPG = tau_SUPG*dot(u, grad(vc))*r*dx F += term_SUPG return lhs(F), rhs(F) - + def advection_diffusion_three_component(self,W,mesh,sol_prev,dt,f1,K=1.0,zh=Constant((0.0,1.0)),SUPG=1,gam_rho=-1.0,rho_0=1.0): + """ + This function returns the bilinear form for the system of PDES governing + an advection-reaction-two-component flow, the governing + PDEs for Darcy flow are: + div(u) = 0 (1) + phi*u = -k*(grad(p)-rho*zhat) (2) + dc0/dt + dot(u,grad(c0)) = div(grad(c0))/Pe - Da*c0*c1**2/phi + beta*f (3) + and + dc1/dt = - 2*Da*c0*c1**2/phi (4) + where + drho = 1 + gam_rho c0 (5) + where c0 and c1 are concentrations of the reactants in the liquid + and solid, u is the fluid velocity, p is pressure, k is permeability + drho=difference between liquid and solid densities, zhat is a unit + vecotr in vertically upward direction, Pe is Peclet number, Da is + the Dahmkoler number, beta is source strength, f is a function + for lateral variations in source of c0, and phi is the constant porosity + Input: + W : Mixed function space containing velocity, pressure and two + concentrations + mesh : Fenics mesh on which W is defined + sol_prev : All unknowns from the previous time step + dt : time step + f1 : Spatially variable scalar function for the source term + K : Constant permeability + zh : Unit vector in the vertical direction + SUPG : A counter for chosing between different SUPG terms + gam_rho : The coefficient of concentration in the expression for drho + rho_0 : The intercept in the expression for drho + Output: + lhs(F) : Left hand side of the combined bilinear form + rhs(F) : Right hand side of the bilinear formulation + The bilinear form uses Crank-Nicholson time stepping and gives + several options for SUPG stabilization. We recommend using Sendur 2018 formulation + for the SUPG term. + The bilinear form, combined together, becomes + phi*dot(u,v)-k*p*div(v)+div(u)*q + eta*(c0n-c0nm1)+eta*dot(u,grad(c_mid))*dt + +dt*dot(grad(c_mid),grad(eta))/Pe+(c1_n-c1_nm1)*omega + = K*drho*dot(zhat,v)-Da*c0*c1*dt*eta/phi + - beta*f1*dt*eta + Da*c0*c1*dt*omega/(1-phi) + In this code c0 is uc, c1 is cc, eta is vc, omega is qc + We recommend using this function preferentially over the other functions as + this combines all four equations into one bilinear form. If the density contrast + is constant, just replace drho with a constant value. + """ + h = CellDiameter(mesh) + # TrialFunctions and TestFunctions + U = TrialFunction(W) + (v, q, vc, qc,qc1) = TestFunctions(W) + u, p, uc, cc,cc1 = split(U) + zhat=zh + # Density contrast as a function of concentration + deltarho=rho_0-gam_rho*uc + + # Define the variational form + F = (inner(self.phi*u,v) - K*div(v)*p+div(u)*q)*dx\ + - K*deltarho*inner(v,zhat)*dx + # uc and cc are the trial functions for the next time step + # uc for comp cc and d comp1 + # u0 (component 0), c0(component 1), and c1(component 3) + # are known values from the previous time step + u,p,u0 ,c0,c1 = split(sol_prev) + # Mid-point solution for comp 0 + u_mid = 0.5*(u0 + uc) + # First order reaction term + # For chemical reaction + # MgCO3+2Fe = C + 3 Fe(2/3)Mg(1/3)O + # Assuming forward reaction controls the rate + # Total molar mass is calculated from + # MgCO3 = 24.0+12.0+3.0*16.0 + # 2Fe = 2.0*56 + # C = 12.0 + # 3 Fe(2/3)Mg(1/3)O = 3.0*(0.67*56.0+0.33*24.0+16.0) + # M = MgCO3+2*Fe+C+3*Oxide = 392.32 + Total_molar_mass = 392.32 + Fe_frac = 2.0*56.0/Total_molar_mass + carbonate_frac = 84.0/Total_molar_mass + C_frac = 12.0/Total_molar_mass + f21 = Fe_frac*self.Da*u0*c0*c0 #reaction rate for Fe + f11 = carbonate_frac*self.Da*u0*c0*c0 #reaction rate for carbonate + f31 = C_frac*self.Da*u0*c0*c0 + + F += vc*(uc - u0)*dx + dt*(vc*dot(u, grad(u_mid))*dx\ + + dot(grad(vc), grad(u_mid)/self.Pe)*dx) \ + + dt*f11/self.phi*vc*dx - self.beta*dt*f1*vc*dx \ + + qc*(cc - c0)*dx + dt*f21/(1-self.phi)*qc*dx \ + + qc1*(cc1 - c1)*dx - dt*(f31/(1-self.phi))*qc1*dx + # Residual + h = CellDiameter(mesh) + r = uc - u0 + dt*(dot(u, grad(u_mid)) - div(grad(u_mid))/self.Pe+f11/self.phi)\ + - self.beta*dt*f1 + cc-c0 + dt*f21/(1.0-self.phi) + cc1-c1 - dt*f31/(1.0-self.phi) + # Add SUPG stabilisation terms + # Default is Sendur 2018, a modification of Codina, 1997 also works + vnorm = sqrt(dot(u, u)) + + if SUPG==1: + #Sendur 2018 + tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h) + elif SUPG==2: + #Codina 1997 eq. 114 + tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h+self.Da) + #tau_SUPG = 1.0/(4.0/(self.Pe*h*h)+2.0*vnorm/h+self.Da*np.max(c0)/self.phi) + else: + alpha_SUPG=self.Pe*vnorm*h/2.0 + #Brookes and Hughes + coth = (np.e**(2.0*alpha_SUPG)+1.0)/(np.e**(2.0*alpha_SUPG)-1.0) + tau_SUPG=0.5*h*(coth-1.0/alpha_SUPG)/vnorm + term_SUPG = tau_SUPG*dot(u, grad(vc))*r*dx + F += term_SUPG + return lhs(F), rhs(F) ################################################################### ### Advection diffusion equation in Stokes flow diff --git a/src/run/LVL_3comp/LVL_3comp.py b/src/run/LVL_3comp/LVL_3comp.py new file mode 100644 index 0000000..e0f4244 --- /dev/null +++ b/src/run/LVL_3comp/LVL_3comp.py @@ -0,0 +1,246 @@ +################################################### +## This is a part of Mupopp +## Copyright Saswata Hier-Majumder, July, 2016 +## Modified on August 2018 +## This program solves an advection diffusion problem +## with Darcy flow, using Dirichlet boundary conditions +## for velocity and concentration +## Modified by Joe Sun, February 2019 +#################################################### + +from fenics import * +from mshr import* +import numpy, scipy, sys, math +#import matplotlib.pyplot as plt +#Add the path to Mupopp module to the code +sys.path.insert(0, '../../modules/') +from mupopp import * + + +##################################################### +parameters["std_out_all_processes"]=False +#set_log_level(30) #Critical info only +#################################### + +# Parameters for initializing the object +Da0 = 100.0 +Pe0 = 100.0 +alpha0= 0.01 # beta = alpha0/phi +phi0 = 0.01 +beta=alpha0/phi0 +Fe=0.01 +cfl0 = 0.01 +# Parameters for iteration +T0 = 2.0 +dt0 = 0.1 +out_freq0 = 1 + +# Parameters for mesh +mesh_density = 40 + +# Output files for quick visualisation + +file_name = "Da_%3.2f_Pe_%.1E_beta_%3.2f_Fe_%3.2f"%(Da0,Pe0,beta,Fe) +output_dir = "output/" + +extension = "pvd" # "xdmf" or "pvd" + +velocity_out = File(output_dir + file_name + "_velocity." + extension, "compressed") +pressure_out = File(output_dir + file_name + "_pressure." + extension, "compressed") +c0_out = File(output_dir + file_name + "_concentration0." + extension, "compressed") +c1_out = File(output_dir + file_name + "_concentration1." + extension, "compressed") +c2_out = File(output_dir + file_name + "_concentration2." + extension, "compressed") +initial_c0_out = File(output_dir + file_name + "_initial_c0." + extension, "compressed") +initial_c1_out = File(output_dir + file_name + "_initial_c1." + extension, "compressed") + +# Output parameters +def output_write(mesh_density,Da,phi,Pe,alpha,cfl,fname= output_dir + "/a_parameters.out"): + """This function saves the output of iterations""" + file=open(fname,"a") + file.write("####################################") + file.write("\n") + file.write("Mesh density: %g\n" %mesh_density) + file.write("Da: %g\n" %Da0) + file.write("phi: %g\n" %phi0) + file.write("Pe: %g\n" %Pe0) + file.write("alpha: %g\n" %alpha0) + file.write("cfl: %g\n" %cfl0) + file.write("####################################") + file.close + +output_write(mesh_density,Da0,phi0,Pe0,alpha0,cfl0) + +#Define function for source term in Governing equations +class SourceTerm(Expression): + """ Creates an expression for the source term + in the advection reaction equations. + The source term consists of a series of + sine waves. + """ + def __init__(self, mesh,element): + self.mesh = mesh + self.element=element + def eval(self, values, x): + g1=x[0]*0.0 + for ii in range(0,20): + g1+=0.1*np.abs(np.sin(ii*x[0]*np.pi)) + + g = (1.0-tanh(x[1]/0.01))*g1 + values[0] = g + def value_shape(self): + return (1,) + +# Define function for BC +class BoundarySource(Expression): + def __init__(self, mesh,element): + self.mesh = mesh + self.element=element + def eval_cell(self, values, x, ufl_cell): + cell = Cell(self.mesh, ufl_cell.index) + n = cell.normal(ufl_cell.local_facet) + g1=x[0]*0.0 + for ii in range(0,20): + g1+=0.01*np.abs(np.sin(ii*x[0]*np.pi)) + g = -g1 + values[0] = g*n[0] + values[1] = g*n[1] + def value_shape(self): + return (2,) + + +############################ +## Numerical solution +############################ + +# Define the mesh +xmin = 0.0 +xmax = 4.0 +ymin = 0.0 +ymax = 1.0 +domain = Rectangle(Point(xmin,ymin),Point(xmax,ymax)) +mesh = generate_mesh(domain,mesh_density) +#mesh = RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), 100, 50) +##################### + +# Mark facets +## Define subdomain for flux calculation +class Plane(SubDomain): + def inside(self, x, on_boundary): + return x[1] > ymax - DOLFIN_EPS +facets = FacetFunction("size_t", mesh) +Plane().mark(facets, 1) +ds = Measure("ds")[facets] +n = FacetNormal(mesh) +#################### +# Define essential boundary +def top_bottom(x): + return x[1] < DOLFIN_EPS or x[1] > ymax - DOLFIN_EPS +def bottom(x): + return x[1] < DOLFIN_EPS +def top(x): + return x[1] > ymax - DOLFIN_EPS + +# Sub domain for Periodic boundary condition +class PeriodicBoundary(SubDomain): + # Bottom boundary is "target domain" G + def inside(self, x, on_boundary): + return bool(near(x[0], xmin) and on_boundary) + # Map right boundary (H) to left boundary (G) + def map(self, x, y): + if near(x[0], xmax): + y[0] = x[0] - xmax + y[1] = x[1] + else: + y[0] = -1000 + y[1] = -1000 +# Create periodic boundary condition +pbc = PeriodicBoundary() + +############################ +## Darcy velocity +############################ +# Create FunctionSpaces +# Velocity +V = VectorElement("Lagrange", mesh.ufl_cell(), 2) +# Pressure +Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1) +#Concentration +Qc = FiniteElement("Lagrange",mesh.ufl_cell(), 1) +# Make a mixed space +W = dolfin.FunctionSpace(mesh, MixedElement([V,Q,Qc,Qc,Qc]), constrained_domain=pbc) +X = FunctionSpace(mesh,"CG",1, constrained_domain=pbc) + + +# Define boundary conditions +G=BoundarySource(mesh,element=V) +#bc1 = DirichletBC(W.sub(0), G, bottom) +bc1 = DirichletBC(W.sub(0), Constant((0.0,0.1)), bottom) +bc = [bc1] + +########################### +## Create an object +########################### +darcy = DarcyAdvection(Da=Da0,phi=phi0,Pe=Pe0,alpha=alpha0,cfl=cfl0) + +# Define initial conditions +sol_0 = Function(W) +temp2 = Function(X) +c01_temp=Expression("0.01",element=Q) #Fe +temp2.interpolate(c01_temp) +c01 = temp2 +assign(sol_0.sub(3), c01) + +#c0_initial,c1_initial=c0.split() +#initial_c0_out << c0_initial +#initial_c1_out << c1_initial + +########################### +## Solve for Darcy velocity +########################### +sol = Function(W) +# Parameters for iteration +T = T0 +dt = dt0 +t = dt +flux=np.array([]) +carbon=np.array([]) +time_array=np.array([]) +i = 1 +out_freq = out_freq0 +S=SourceTerm(mesh,element=Qc) + +while t - T < DOLFIN_EPS: + # Update the concentration of component 0 + a,L = darcy.advection_diffusion_three_component(W,mesh,sol_0,dt,f1=S,K=0.1) + solve(a==L,sol,bc) + sol_0 = sol + u0,p0,c00,c01,c02 = sol.split() + if i % out_freq == 0: + u0.rename("velocity","") + velocity_out << u0 + p0.rename("pressure","") + pressure_out << p0 + c00.rename("[CO3]","") + c0_out << c00 + c01.rename("[Fe]","") + c1_out << c01 + c02.rename("[C]","") + c2_out << c02 + ## Calculate flux + flux1 = assemble(c00*phi0*dot(u0, n)*ds(1)) + flux= np.append(flux,flux1) + time_array=np.append(time_array,t) + #calculate mass of carbon deposited + mass1=assemble(c02*(1.0-phi0)*dx) + carbon = np.append(carbon,mass1) + #print "flux 1: ", flux_1 + # Move to next interval and adjust boundary condition + info("time t =%g\n" %t) + info("iteration =%g\n" %i) + #print 'iteration',i + t += dt + i += 1 +flux_file=output_dir + file_name + "_flux.csv" +np.savetxt(flux_file,(time_array,flux),delimiter=',') +mass_file=output_dir + file_name + "_mass.csv" +np.savetxt(mass_file,(time_array,carbon),delimiter=',')