**2. Multi-Physics Problem** =========================================== .. note:: The reader is encouraged to look at `this paper `_ to make a sense out of numerical solution of a coupled physics problem: 2.1 Coupled System of Equations ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ In this section, we solve a coupled system of partial differential equations. To be more specific, we look into the **Poisson-Nernst-Planck** equations representing coupled physics including the chemical and electrical fields. The strong from of the Nernst-Planck equations is expressed as: .. math:: :name: eq.79 \frac{\partial c}{\partial t}+ \nabla . (-D_i \nabla c_i- z_i \mu_i F c_i \nabla \phi)=0,\quad i=Na^+,Cl^- Where :math:`c` is concentration of the mobile species, :math:`D` is diffusion constant of the mobile species, :math:`\mu = \frac{D}{RT}` is mobility of the ionic species ,:math:`R` is universal gas constant , :math:`T` is temperature, :math:`F` is Faraday constant and :math:`\phi` is Electrical potential. .. math:: :name: eq.80 \nabla^2 \phi = -\frac{F}{\varepsilon_0 \varepsilon_r}(z_{Na} c_{Na} + z_{Cl} c_{Cl} + z_{fixed} c_{fixed}) The approximate solution is sought on scalar spaces including :math:`S^{(c)}` and :math:`S^{(\phi)}` corresponding to the function spaces representing the concentration of the ions and electric potential respectively. .. math:: :name: eq.81 S^{c}=span \{\omega_0^{c} ,...,\omega_N^{c} \} S^{\phi}=span \{\lambda_0^{\phi} ,...,\lambda_N^{\phi} \} The :math:`\omega_i` and :math:`\lambda_i` are the basis functions which are linearly independent. This approximation is used to discretize the problem in space. In this regard the approximate solution and can be expressed as a linear combination of the basis functions: .. math:: :name: eq.82 c= \sum_{j=0}^N \alpha_j \omega_j (x,y) \phi= \sum_{j=0}^N \beta_j \lambda_j (x,y) The coefficients :math:`\alpha_i` and :math:`\beta_i` should be computed. A mixed scalar space is then defined as: :math:`\textbf{S} = S^{(c)} \times S^{(\phi)}` along with a test function :math:`\textbf{v} = (v_c,v_{\phi})` .The test functions are chosen as :math:`\forall v_c , v_{\phi} \in \textbf{S}`, So we can define the test function: :math:`v_c = \omega_i` and :math:`v_{\phi} = \lambda_i` . In order to derive the weak form of the NP equation, we should multiply it by :math:`v_c`: .. math:: :name: eq.83 \int_{\Omega} v_c [\frac{\partial c}{\partial t}+ \nabla . (-D \nabla c- z \mu F c \nabla \phi)]dx Then the above equation could be expanded as following: .. math:: :name: eq.84 \int_{\Omega} v_c \frac{\partial{c}}{\partial t} dx+\int_{\Omega} v_c D {\nabla}^2{c} dx - \int_{\Omega} v_c \mu F z \nabla{c} \nabla{\phi} dx -\int_{\Omega} v_c \mu F z {\nabla}^2{\phi} dx = 0 By integration by part and applying the divergence theorem, it yields to: .. math:: :name: eq.85 \int_{\Omega} v_c \frac{\partial c}{\partial t} dx - \int_{\partial \Omega} v_c D \nabla{c}.n.ds + \int_{\Omega} D(\nabla c.\nabla \phi) dx-\int_{\Omega} v_c \mu Fz \nabla{c} \nabla{\phi} dx - \int_{\partial \Omega} v_c \mu Fzc \nabla{\phi}.n.ds+\int_{\Omega} \mu Fz \nabla(cv_c).\nabla \phi.dx=0 In the above equation, :math:`n` is the outward normal vector on the boundaries. The terms including :math:`ds` correspond to the Neumann boundary condition that vanish as they are equal to zero: .. math:: :name: eq.86 \int_{\Omega} v_c \frac{\partial c}{\partial t} dx + \int_{\Omega} D(\nabla c.\nabla \phi) dx- \int_{\Omega} v_c \mu Fz \nabla{c} \nabla{\phi} dx +\int_{\Omega} v_c \mu F(\nabla c.\nabla\phi)dx+ \int_{\Omega} \mu Fcz \nabla(v_c).\nabla\phi.dx=0 In the above equation, the third and fourth terms cancel out: .. math:: :name: eq.87 \int_{\Omega} v_c \frac{\partial c}{\partial t}dx+\int_{\Omega}D(\nabla c.\nabla v_c)dx+\int_{\Omega} \mu zcF(\nabla{\phi}.\nabla{v_c})dx = 0 The derivative in time could be written according to the backward-Euler scheme: .. math:: :name: eq.88 \frac{\partial c}{\partial t} = \frac{c-c^0}{\Delta t} Then after plugging into the :ref:`equation.87 `: .. math:: :name: eq.89 \int_{\Omega} v_c c dx-\int_{\Omega}v_c c^0 dx+\Delta t \int_{\Omega}D(\nabla c.\nabla{v_c})dx+\Delta t \int_{\Omega}\mu z Fc (\nabla \phi.\nabla{v_c})dx=0 The same approach could be taken by deriving the weak form of the Poisson equation after multiplying it into a test function :math:`v_{\phi}`: .. math:: :name: eq.90 \int_{\Omega} v_{\phi}\nabla^2 \phi dx = -\int_{\Omega} v_{\phi} [-\frac{F}{\varepsilon_0 \varepsilon_r}(z_{Na} c_{Na} + z_{Cl} c_{Cl} + z_{fixed} c_{fixed})]dx After integration by part and applying the divergence theorem, it yields to: .. math:: :name: eq.91 \int_{\partial \Omega} v_{\phi}\nabla \phi n.ds -\int_{\Omega} \nabla v_{\phi}.\nabla \phi dx = \int_{\Omega} v_{\phi} [-\frac{F}{\varepsilon_0 \varepsilon_r}(z_{Na} c_{Na} + z_{Cl} c_{Cl} + z_{fixed} c_{fixed})]dx The first term in the left-hand side of the above equation, corresponds to the Neumann boundary condition that is equal to zero. So, this could be rewritten by expanding the right-hand side to find the final weak from of the Poisson equation: .. math:: :name: eq.92 \int_{\Omega} \nabla v_{\phi} \nabla \phi dx - \int_{\Omega} v_{\phi} \frac{F}{\varepsilon_0 \varepsilon_r}(z_{Na} c_{Na})dx- \int_{\Omega} v_{\phi} \frac{F}{\varepsilon_0 \varepsilon_r}(z_{Cl} c_{Cl})dx - \int_{\Omega} v_{\phi} \frac{F}{\varepsilon_0 \varepsilon_r}(z_{fixed} c_{fixed})dx = 0 The constants used in the simulation as well as the initial conditions are presented in below tables: .. csv-table:: Constants :widths: 5, 5, 5,5,5,5,5,5,5,5 :align: center :math:`F(\frac{c}{mol})`, :math:`\epsilon_0(\frac{F}{m})`,:math:`\epsilon_r`,:math:`z_{Na}`,:math:`z_{Cl}`,:math:`z_{fixed}`,:math:`R(\frac{J}{mol.K})``,:math:`T(K)` ,:math:`D_{Na}(\frac{m^2}{s})`,:math:`D_{Cl}(\frac{m^2}{s})` 96487, :math:`8.85 \times 10^{-12}` , 100,1,-1,-1,8.3114,293,:math:`10^{-7}`,:math:`10^{-7}` .. csv-table:: Initial Condition :widths: 5,5 :align: center Solution, Gel :math:`c_{Na}=1mM`,:math:`c_{Na}=5.193mM` :math:`c_{Cl}=1mM`,:math:`c_{Cl}=0.193mM` ,:math:`c_{fixed}=1mM` 2.2 Finite Element Implementation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. note:: In this section we replicate the results in `this paper `_ The summary of the FEM code implemented in FeniCS is presented here: .. code-block:: python from dolfin import * Then we should define the constants: .. code-block:: python # Diffusion constant for Na+ cations GEL D_Na = 1.0 * (10 ** (-7.)) # Diffusion constant for Cl- anioins GEL D_Cl = 1.0 * (10 ** (-7.)) # Diffusion constant for Na+ cations SOLUTION D_Na_s = 1.0 * (10 ** (-7.)) # Diffusion constant for Cl- anioins SOLUTION D_Cl_s = 1.0 * (10 ** (-7.)) # The gas constant R = 8.31 # Initial value for Na+ in the solution c_init_Na_sol = 1 # Initial value for Cl- in the solution c_init_Cl_sol = 1 # Initial value for fixed anioin (e.g. bound charge) in the gel c_init_Fixed_gel = 5 c_init_Na_gel = ((c_init_Fixed_gel + pow(c_init_Fixed_gel ** 2 + 4 * pow(c_init_Cl_sol, 2), 0.5))) / 2. c_init_Cl_gel = ((-c_init_Fixed_gel + pow(c_init_Fixed_gel ** 2 + 4 * pow(c_init_Cl_sol, 2), 0.5))) / 2. # Valence of Na+ z_Na = 1.0 # Valence of Cl- z_Cl = -1.0 # Valence of Fixed anioin z_Fixed = -1.0 # Anion concentration # C_0 = Constant(500.0) eps_vacume = 8.85 * pow(10, -12) eps_water = 100. eps_sclera = 100. # Dielectric permittivity epsilon = eps_vacume * eps_water epsilon_p = eps_vacume * eps_sclera # Temperature Temp = 293.0 # Electrical mobility for Na+ cations mu_Na = D_Na / (R * Temp) # Electrical mobility for Cl- cations mu_Cl = D_Cl / (R * Temp) mu_Na_s = D_Na_s / (R * Temp) # Electrical mobility for Cl- cations mu_Cl_s = D_Cl_s / (R * Temp) # Faraday number Faraday = 96485.34 Next we should present the mesh in the simulation. The domains include a hydrogel (Shown in red color with the dimensions of :math:`4 \times 10 mm^2`) placed in the middle of a solution domain (Shown in green color with the dimensions of :math:`50 \times 50 mm^2`). The mesh was created in the GMSH. .. note:: The density of the element should be significantly increased near the boundaries of the hydrogel and solution to capture the steep gradient of the concentration of the mobile ions The domains and the mesh are shown here: .. figure:: PNG/11.png :align: center Computational Domains (Left) and Mesh (Right) .. note:: The mesh file is available in the github repository in the folder **PNP-MESH**. .. code-block:: python mesh = Mesh() hdf = HDF5File(mesh.mpi_comm(), "0file.h5", "r") hdf.read(mesh, "/mesh", False) subdomains = MeshFunction('size_t', mesh, mesh.topology().dim()) hdf.read(subdomains, "/subdomains") boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1) hdf.read(boundaries, "/boundaries") The integration symbol :math:`dx` is defined as: .. code-block:: python dx = Measure('dx', domain=mesh, subdomain_data=subdomains) Next we should define the appropriate element and then form the mixed-function space containing all unknowns: .. code-block:: python # Defining element for scalar variables (e.g. concentrations and voltage) Element1 = FiniteElement("CG", mesh.ufl_cell(), 1) # Defining the mixed function space W_elem = MixedElement([Element1, Element1, Element1]) W = FunctionSpace(mesh, W_elem) The summary of the FEniCS code is presented here: .. code-block:: python from dolfin import * # Diffusion constant for Na+ cations GEL D_Na = 1.0 * (10 ** (-7.)) # Diffusion constant for Cl- anioins GEL D_Cl = 1.0 * (10 ** (-7.)) # Diffusion constant for Na+ cations SOLUTION D_Na_s = 1.0 * (10 ** (-7.)) # Diffusion constant for Cl- anioins SOLUTION D_Cl_s = 1.0 * (10 ** (-7.)) # The gas constant R = 8.31 # Initial value for Na+ in the solution c_init_Na_sol = 1 # Initial value for Cl- in the solution c_init_Cl_sol = 1 # Initial value for fixed anioin (e.g. bound charge) in the gel c_init_Fixed_gel = 5 ############################ ############################ c_init_Na_gel = ((c_init_Fixed_gel + pow(c_init_Fixed_gel ** 2 + 4 * pow(c_init_Cl_sol, 2), 0.5))) / 2. c_init_Cl_gel = ((-c_init_Fixed_gel + pow(c_init_Fixed_gel ** 2 + 4 * pow(c_init_Cl_sol, 2), 0.5))) / 2. ############################## ############################## # Valence of Na+ z_Na = 1.0 # Valence of Cl- z_Cl = -1.0 # Valence of Fixed anioin z_Fixed = -1.0 # Anion concentration # C_0 = Constant(500.0) eps_vacume = 8.85 * pow(10, -12) eps_water = 100. eps_sclera = 100. # Dielectric permittivity epsilon = eps_vacume * eps_water epsilon_p = eps_vacume * eps_sclera # Temperature Temp = 293.0 # Electrical mobility for Na+ cations mu_Na = D_Na / (R * Temp) # Electrical mobility for Cl- cations mu_Cl = D_Cl / (R * Temp) mu_Na_s = D_Na_s / (R * Temp) # Electrical mobility for Cl- cations mu_Cl_s = D_Cl_s / (R * Temp) # Faraday number Faraday = 96485.34 mesh = Mesh() hdf = HDF5File(mesh.mpi_comm(), "0file.h5", "r") hdf.read(mesh, "/mesh", False) subdomains = MeshFunction('size_t', mesh, mesh.topology().dim()) hdf.read(subdomains, "/subdomains") boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1) hdf.read(boundaries, "/boundaries") # File("domains.pvd") << subdomains # lable domain dx = Measure('dx', domain=mesh, subdomain_data=subdomains) . . . # Varitional form for the Nernst-Planck equation in the hydrogel domain for Na+ and Cl- Weak_NP_Na = ... Weak_NP_Cl = ... # Varitional form for the Poisson equation in the hydrogel domain for Na+ and Cl- and fixed charge Weak_Poisson = ... ################################## # Varitional form for the Nernst-Planck equation in the solution domain for Na+ and Cl- Weak_NP_Na_sol = ... Weak_NP_Cl_sol = ... # Varitional form for the Poisson equation in the solution domain for Na+ and Cl- and fixed charge Weak_Poisson_sol = ... # Summing up variational forms F = Weak_NP_Na + Weak_NP_Cl + Weak_Poisson + Weak_NP_Cl_sol + Weak_NP_Na_sol + Weak_Poisson_sol . . . Na_Ion = File("Na.pvd") Cl_Ion = File("Cl.pvd") Voltage = File("phi.pvd") while t <= T: bc_left_volt = DirichletBC(W.sub(2), -0.1, boundaries, 35) bc_right_volt = DirichletBC(W.sub(2), 0.1, boundaries, 36) bcs = [bc_left_volt, bc_right_volt, bc_left_c, bc_right_c, bc_left_c_1, bc_right_c_1] J = derivative(F, z, dz) problem = NonlinearVariationalProblem(F, z, bcs, J) solver = NonlinearVariationalSolver(problem) solver.solve() (c_Na, c_Cl, phi) = z.split(True) C_previous_Na.assign(c_Na) C_previous_Cl.assign(c_Cl) t += dt Na_Ion << c_Na Cl_Ion << c_Cl Voltage << phi The results could be visualized in Paraview. Here are the results after 10s for the concentration of the sodium ion, chloride ion and electrical voltage: .. figure:: PNG/12.png :align: center Initial and stationary solutions The figures (a), (c) and (e) correspond to the initial values of the :math:`c_{Na}` , :math:`c_{Cl}` and :math:`\phi` respectively. The figures (b), (d) and (f) correspond to the stationary values of the :math:`c_{Na}` , :math:`c_{Cl}` and :math:`\phi` respectively. .. note:: The above code has been written with the capability for parallel computation. You can significantly increase the speed if running it in parallel. In order to run the code in parallel, all you need to do is just opening the terminal and then run this command: .. code-block:: python mpirun -np 4 python3 code.py Where 4 is the number of the cores you want to use to run the code. You can increase it based on number of the cores you have access to. .. note:: The complete FEniCS code is available upon request (Please refer to the section : **About the Author**)