BME/ME 506
Computational Modeling
of Biological Tissues

Multilevel Homogenization Formulation for the Diffusion Problem

1.         Preliminaries: Governing Partial Differential Equation for the Diffusion Problem

            Diffusion of chemical species like oxygen plays a significant role in many physiological processes. If we consider a fixed volume, then the chemical species within that volume must be conserved. In other words, the rate of change of the chemical species u within the volume must equal the local production of u within the volume plus the next flux of u across the boundary. We can write this conservation law mathematically as:

where u is the chemical species concentration, the left hand side is the rate of change of u, f represents the local production of u, Ji is the flux of u, ni is the normal to the boundary S of the fixed volume V. If J is smooth and continuous, we can apply the divergence theorem to the flux Ji:

Applying this result, we can rewrite the basic conversation law using all volume integrals as:

We can move the time derivative inside the left hand side integral. We can then drop the integral in the conservation law and obtain the following partial differential equation for diffusion:

Like all physics problems, we have a relationship between the flux of a quantity and its spatial gradient through inherent material properties. This relationship is called a constitutive relationship. For diffusion equations, if the chemical species concentration is not too high, then we can specify a linear relationship between the chemical species flux and the chemical species gradient. For diffusion problems, this particular constitutive relationship is known as Fick's Law. Fick's law can be stated as follows:

If we substitute Fick's law into the governing diffusion differential equation, we obtain:

We note that the diffusion equation has a very similar form to the potential differential equations with the exception of the time dependent term.

2.         Homogenization Expansion of the Governing Diffusion Differential Equation

            The existence of multiple scales in diffusion problems can be traced to the different diffusion rates along cells, gap junctions between cells and tissue matrix. Chemical species can diffuse readily within cell cytoplasm, but not easily across gap junctions which are highly resistive to flow. Thus, diffusion at the cell and sub-cell level will very rapidly. At a more macroscopic level, we also have diffusion over the space of many cells and tissue matrix. To account for this feasibily, we must again turn to a multiple scale homogenization analysis. As with the potential equation homogenization formulation, we have three basic premises: 1) relationship between the structural scales, 2) the hierarchical gradient, and 3) asymptotic expansion of the chemical species concentration.

i.          The relationship between the different structural scales is written in homogenization theory as:

where xi is the more microscopic scale, xi-1 is the more macroscopic scale, and e is the ratio of characteristic lengths of the microscopic level to the macroscopic level as:

where |li| is the characteristic length scale of the microscopic level, |li-1| is the characteristic length scale of the macroscopic level, and e is the ratio of the length scales.

ii.         Hierarchical Gradient

            By using the chain rule and the above structural scale relationship, we can write the hierarchical gradient as:

iii.         Asymptotic Expansion of Chemical Species Concentration

            As with the potential equation, we use an asymptotic expansion of the chemical species concentration:

If we use the basic homogenization concepts i,ii, and iii in the governing diffusion partial differential equation we obtain:

We next assume that solving each separate equation related to order e will provide the solution to the above equation. When we separate equations of order e we obtain:

The equation associated with e-2 leads to the conclusion that u0 does not depend on the microscopic level, as shown for the potential problem. This leaves the differential equation associated with e-1 as:

We typically neglect terms of order e, so we neglect the time rate of change of u1. This gives the following differential equation associated with e0:

The above two differential equations represent the microscopic and macroscopic differential equations, respectively.

3.         Finite Element Formulation for the Homogenization Equations

            Using the above microscopic and macroscopic equations for the multiscale diffusion problem, we can develop a weak formulation and go on to write the finite element approximation. Let us first consider the microscopic differential equation. Since this equation has the same form as the potential microscopic equation, we used the same weak form:

where v is now the virtual chemical species concentration. If we make the same arguments for linear superposition as in the potential case, we end up with the same form of the auxillary equations:

             Auxillary Eq. 1

              Auxillary Eq. 2

               Auxillary Eq. 3

By making the appropriate substitution of shape functions as in the potential formulation, we derive the following finite element equations:

By solving the above equations, we can then compute an effective Diffusion coefficient as:

Using the above expression for the effective Diffusion tensor, we can then write the weak form for the global equation as:

This global equation is also similar to the global potential equation, with the exception that we have a time dependent term on the left hand side.   First we note that we can interpolate the time rate of change of u using the same shape functions with which we interpolate u0, the virtual concentration and the coordinates:

When we apply the interpolation or shape functions to the time rate of change term we have:

where the integral is performed over each element domain , after which the elements are summed. Since both v and du0/dt are constant over the element domain, they can be pulled out of the element integral. If we write the remaining terms inside the integral as:

By using the above matrix notation with following matrices for the effective diffusion and rate of production of a chemical species f:

We thus have the following finite element formulation for the time dependent diffusion problem:

where nelt is the total number of elements in the mesh, Me is the matrix associated with the rate of change of chemical species concentration, du0/dt is the rate of change of chemical species concentration, [Ke] is the diffusion coefficient matrix based on Fick's law, {u0} is the concentration of the chemical species, {fe} is the rate of local production of the chemical species and the boundary flux of the chemical species.

We now have to decide how to handle the time dependent term in the diffusion equation. The way time dependent terms are generally handled in finite element analysis is through finite difference approximation of the time derivative. Using finite difference, we can approximate the time derivative in a time dependent differential equation as:

where  is the value of u0 at time t+Dt,   is the value of u0 at t, Dt is the size of a time step, f is the remainder of the function in the time differential equation, and q is a parameter chosen for the time integration.   The above equation is a time marching algorithm to integrate the differential equation over time. We make the following substitutions into the time dependent finite element formulation for the diffusion equation:

we can thus write the finite element implementation using the above equations as:

In the above equation,  is the chemical species concentration that we need to calculate and  is the known chemical species concentration from the previous time step. Thus, we can place all terms associated with  on the right hand side of the equation because they are all known. This leaves:

We note that the above equation takes the familiar form of the standard linear algebraic equations we have obtained previously from finite element formulations if we make the following grouping of terms:

So if we make the following matrix notations:

we then have:

We note that these equations must be solved at each time step.

At this point, we have not said anything about the parameter q. This parameter will take on a value between 0 and 1. The choice of this parameter will affect the accuracy and stability of the time integration scheme. If we choose q = 0, then we have no diffusion cofficient term on the left hand side of the finite element equation. This time stepping algorithm is termed explicit, since we don't have to invert a diffusion coefficient matrix. The advantage of this method is that it is computationally less costly. However, this algorithm is conditionally stable, which means that we have to take a very small time step for the algorithm to work properly.

4.         Summary of the Homogenization Formulation for Time Dependent Diffusion

            In summary, we have the same four components for multiscale homogenization analysis of the time dependent diffusion equation as we do for the potential formulation:

i.           Calculate characteristic chemical concentration u* by solving local equation at RVE level:

            where the matrices are defined previously.

ii.          Compute the effective diffusion coefficients based on characteristic chemical concentration gradient (HOMOGENIZATION):

iii.         Solve the time dependent global diffusion problem using the effective diffusion coefficient:

iv.         Compute the local chemical concentration gradient for each time step, since u0 is time dependent (LOCALIZATION):

BME 506 Home