From cb995c7d0049c83e199a9ae307979197f068af54 Mon Sep 17 00:00:00 2001 From: Zhe Rex Li Date: Thu, 12 May 2022 15:39:50 +1000 Subject: [PATCH] add sample files for plane membrane --- example/PlaneMembrane/plane_membrane.db | 117 ++++++++++++++++++++++++ example/PlaneMembrane/plane_membrane.py | 90 ++++++++++++++++++ 2 files changed, 207 insertions(+) create mode 100644 example/PlaneMembrane/plane_membrane.db create mode 100755 example/PlaneMembrane/plane_membrane.py diff --git a/example/PlaneMembrane/plane_membrane.db b/example/PlaneMembrane/plane_membrane.db new file mode 100644 index 00000000..3414af24 --- /dev/null +++ b/example/PlaneMembrane/plane_membrane.db @@ -0,0 +1,117 @@ +MultiphysController { + timestepMax = 20000 + visualization_interval = 1000 // Frequency to write visualization data + analysis_interval = 20 // Frequency to perform analysis +} +Stokes { + epsilonR = 78.5 //fluid dielectric constant [dimensionless] + tau = 1.0 + F = 0, 0, 0 + rho_phys = 998.2 + nu_phys = 1.003e-6 //fluid kinematic viscosity; user-input unit: [m^2/sec] + BC = 3 // Pressure constant BC + din = 1.0 // Inlet pressure + dout = 1.0 // Outlet pressure + UseElectroosmoticVelocityBC = true + SolidLabels = 0, -1 + ZetaPotentialSolidList = -0.005, -0.03 // unit [v] +} +Ions { + temperature = 310.15 //unit [K] + //number_ion_species = 5 //number of ions + //tauList = 1.0, 1.0, 1.0, 1.0, 1.0 // H+, OH-, Na+, Cl-, Fe3+ + //IonDiffusivityList = 9.3e-9, 5.3e-9, 1.3e-9, 2.0e-9, 0.604e-9 //user-input unit: [m^2/sec] + //IonValenceList = 1, -1, 1, -1, 3 //valence charge of ions; dimensionless; positive/negative integer + //IonConcentrationList = 1.0e-4, 1.0e-4, 100, 100, 0 //user-input unit: [mol/m^3] + number_ion_species = 2 //number of ions + //IonConcentrationFile = "Pseudo3D_plane_membrane_concentration_Na_z192_xy64.raw", "double", "Pseudo3D_plane_membrane_concentration_Na_z192_xy64.raw", "double" + tauList = 1.0,1.0 // Na+, anion + IonDiffusivityList = 1e-9,1e-9 //user-input unit: [m^2/sec] + IonValenceList = 1,-1 //valence charge of ions; dimensionless; positive/negative integer + IonConcentrationList = 145e-3,145e-3 //user-input unit: [mol/m^3] + MembraneIonConcentrationList = 15e-3, 15e-3 + BC_InletList = 0,0 //boundary condition for inlet; 0=periodic; 1=ion concentration; 2=ion flux + BC_OutletList = 0,0 //boundary condition for outlet; 0=periodic; 1=ion concentration; 2=ion flux + InletValueList = 15e-3, 15e-3 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] + OutletValueList = 145e-3, 145e-3 //if ion concentration unit=[mol/m^3]; if flux (inward) unit=[mol/m^2/sec] + BC_Solid = 0 //solid boundary condition; 0=non-flux BC; 1=surface ion concentration + //SolidLabels = 0 olid labels for assigning solid boundary condition; ONLY for BC_Solid=1 + //SolidValues = 1.0e-5 // user-input surface ion concentration unit: [mol/m^2]; ONLY for BC_Solid=1 + FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging +} +Poisson { + epsilonR = 80.4 //fluid dielectric constant [dimensionless] + tau = 4.5 + BC_Inlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential + BC_Outlet = 0 // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential + InitialValueLabels = 1,2//a list of labels of fluid nodes + InitialValues = 60.6e-3, 0 //unit: [V] + //------- Boundary Voltage for BC = 1 (Inlet & Outlet) --------------------- + Vin = 60.6e-3 //ONLY for BC_Inlet = 1; electrical potential at inlet + Vout = 0 //ONLY for BC_Outlet = 1; electrical potential at outlet + //-------------------------------------------------------------------------- + //------- Boundary Voltage for BC = 2 (Inlet & Outlet) --------------------- + //Vin0 = 0.01 //(ONLY for BC_Inlet = 2); unit:[Volt] + //freqIn = 1.0 //(ONLY for BC_Inlet = 2); unit:[Hz] + //t0_In = 0.0 //(ONLY for BC_Inlet = 1); unit:[sec] + //Vin_Type = 1 //(ONLY for BC_Inlet = 1); 1->sin(); 2->cos() + //Vout0 = 0.01 //(ONLY for BC_Outlet = 1); unit:[Volt] + //freqOut = 1.0 //(ONLY for BC_Outlet = 1); unit:[Hz] + //t0_Out = 0.0 //(ONLY for BC_Outlet = 1); unit:[sec] + //Vout_Type = 1 //(ONLY for BC_Outlet = 1); 1->sin(); 2->cos() + //-------------------------------------------------------------------------- + BC_SolidList = 1 //solid boundary condition; 1=surface potential; 2=surface charge density + SolidLabels = 0 //solid labels for assigning solid boundary condition + SolidValues = -0.001 //if surface potential, unit=[V]; if surface charge density, unit=[C/m^2] + WriteLog = true //write convergence log for LB-Poisson solver + // ------------------------------- Testing Utilities ---------------------------------------- + // ONLY for code debugging; the followings test sine/cosine voltage BCs; disabled by default + TestPeriodic = false + TestPeriodicTime = 1.0 //unit:[sec] + TestPeriodicTimeConv = 0.01 //unit:[sec] + TestPeriodicSaveInterval = 0.2 //unit:[sec] + //------------------------------ advanced setting ------------------------------------ + timestepMax = 10000 //max timestep for obtaining steady-state electrical potential + analysis_interval = 200 //timestep checking steady-state convergence + tolerance = 1.0e-6 //stopping criterion for steady-state solution +} +Membrane { + MembraneLabels = 1 + VoltageThreshold = 100.0, 100.0 + MassFractionIn = 1,0 + MassFractionOut = 1,0 + ThresholdMassFractionIn = 1, 0 + ThresholdMassFractionOut = 1, 0 + +} +Domain { + Filename = "Pseudo3D_double_plane_membrane_z192_xy64_InsideLabel1_OutsideLabel2.raw" + nproc = 1, 1, 3 // Number of processors (Npx,Npy,Npz) + n = 64, 64, 64 // Size of local domain (Nx,Ny,Nz) + N = 64, 64, 192 // size of the input image + voxel_length = 0.01 //resolution; user-input unit: [um] + BC = 0 // Boundary condition type0 + ReadType = "8bit" + ReadValues = 2, 1 + WriteValues = 2, 1 + //InletLayers = 0, 0, 1 + //OutletLayers = 0, 0, 1 + //InletLayersPhase = 1 + //OutletLayersPhase = 1 + //checkerSize = 3 // size of the checker to use +} +Analysis { +} +Visualization { + save_electric_potential = true + save_concentration = true + #save_velocity = true + #save_pressure = true + save_8bit_raw = true +} + + + + + + diff --git a/example/PlaneMembrane/plane_membrane.py b/example/PlaneMembrane/plane_membrane.py new file mode 100755 index 00000000..a525dda0 --- /dev/null +++ b/example/PlaneMembrane/plane_membrane.py @@ -0,0 +1,90 @@ +import numpy as np +import math +import matplotlib.pyplot as plt + +#physical constant +k_B_const = 1.380649e-23 #[J/K] +N_A_const = 6.02214076e23 #[1/mol] +e_const = 1.602176634e-19 #[C] +epsilon0_const = 8.85418782e-12 #[C/V/m] + +#other material property parameters +epsilonr_water = 80.4 +T=310.15 #[K] + +#input ion concentration +C_Na_in = 15e-3 #[mol/m^3] +C_Na_out = 145e-3 #[mol/m^3] +C_K_in = 150e-3 #[mol/m^3] +C_K_out = 4e-3 #[mol/m^3] +C_Cl_in = 10e-3 #[mol/m^3] +C_Cl_out = 110e-3 #[mol/m^3] + +#calculating Debye length +#For the definition of Debye lenght in electrolyte solution, see: +#DOI:10.1016/j.cnsns.2014.03.005 +#Eq(42) in Yoshida etal., Coupled LB method for simulator electrokinetic flows +prefactor= math.sqrt(epsilonr_water*epsilon0_const*k_B_const*T/2.0/N_A_const/e_const**2) +debye_length_in = prefactor*np.sqrt(np.array([1.0/C_Na_in,1.0/C_K_in,1.0/C_Cl_in])) +debye_length_out = prefactor*np.sqrt(np.array([1.0/C_Na_out,1.0/C_K_out,1.0/C_Cl_out])) +print("Debye length inside membrane in [m]") +print(debye_length_in) +print("Debye length outside membrane in [m]") +print(debye_length_out) + + +#setup domain +cube_length_z = 192 +cube_length_xy = 64 +#set LBPM domain resoluiton +h=0.01 #[um] +print("Image resolution = %.6g [um] (= %.6g [m])"%(h,h*1e-6)) + +domain=2*np.ones((cube_length_z,cube_length_xy,cube_length_xy),dtype=np.int8) +zgrid,ygrid,xgrid=np.meshgrid(np.arange(cube_length_z),np.arange(cube_length_xy),np.arange(cube_length_xy),indexing='ij') +domain_centre=cube_length_xy/2 +make_bubble = np.logical_and(zgrid>=cube_length_z/4,zgrid<=cube_length_z*0.75) +domain[make_bubble]=1 + +##save domain +file_name= "Pseudo3D_double_plane_membrane_z192_xy64_InsideLabel1_OutsideLabel2.raw" +domain.tofile(file_name) +print("save file: "+file_name) + +#debug plot +#plt.figure(1) +#plt.pcolormesh(domain[:,int(domain_centre),:]) +#plt.colorbar() +#plt.axis("equal") +#plt.show() + +##generate initial ion concentration - 3D +#domain_Na = C_Na_out*np.ones_like(domain,dtype=np.float64) +#domain_Na[make_bubble] = C_Na_in +#domain_K = C_K_out*np.ones_like(domain,dtype=np.float64) +#domain_K[make_bubble] = C_K_in +#domain_Cl = C_Cl_out*np.ones_like(domain,dtype=np.float64) +#domain_Cl[make_bubble] = C_Cl_in +# +#domain_Na.tofile("Pseudo3D_plane_membrane_concentration_Na_z192_xy64.raw") +#domain_K.tofile("Pseudo3D_plane_membrane_concentration_K_z192_xy64.raw") +#domain_Cl.tofile("Pseudo3D_plane_membrane_concentration_Cl_z192_xy64.raw") + +##debug plot +#plt.figure(2) +#plt.subplot(1,3,1) +#plt.title("Na concentration") +#plt.pcolormesh(domain_Na[:,int(bubble_centre),:]) +#plt.colorbar() +#plt.axis("equal") +#plt.subplot(1,3,2) +#plt.title("K concentration") +#plt.pcolormesh(domain_K[:,int(bubble_centre),:]) +#plt.colorbar() +#plt.axis("equal") +#plt.subplot(1,3,3) +#plt.title("Cl concentration") +#plt.pcolormesh(domain_Cl[:,int(bubble_centre),:]) +#plt.colorbar() +#plt.axis("equal") +#plt.show()