BME/ME 506
Computational Modeling
of Biological Tissues
Homogenization Formulation for the Potential Problem with Multiple Scales
1. Homogenization Theory Preliminaries
a. Length Scales
In homogenization analysis, at least two different coordinate scales are assumed to exist. The 0th level scale, or macroscopic scale, is denoted as x0. The 1st level scale, or microscopic scale, is denoted as x1. It is important to note that this notation differs from that traditionally used in the homogenization literature. In that literature, the macroscopic scale is denoted by x and the microscopic scale is denoted by y. However, since biological tissues almost inevitably have more than two scales, I have adopted the nomenclature using 0 to represent the most macroscopic scale going to n to represent the most microscopic scale of a tissue with n levels. This allows the formulation to be systematic and when coupled with a specific tissue morphology allows the user to know what scale the are dealing with and where they are within the hierarchy. We assume that the two scales within the hierarchy are related as shown below:
Figure 1. Relationship between steps on microscopic scale x1 and macroscopic scale x0. Note that 10 steps are taken on the microscopic scale for each 1 step on the macroscopic scale.
As stated in the figure legend, we take 10 steps on the microscopic scale for
each single step taken on the macroscopic scale. This means that: 
In terms of length scales, however, the length scale l1 on the microscopic level is 1/10 of the length scale l0 on the macroscopic scale. This can be written as:
were the parameter e in homogenization theory denotes the ratio of characteristic length scales between levels of hierarchy. By inspection, we can easily see that the ratio in steps between the microscopic coordinates and the macroscopic coordinates in eq. 1 is 1/(1/10) or 1/e. Therefore, we can write the formal relationship between the microscopic coordinate position x1 and the macroscopic coordinate position as:
If we assume that their are n levels in a tissue hierarchy, then we can write the general relationship between the ith level (more microscopic level) and the (i-1)th level (more macroscopic level) as:
where e is the previously defined ratio of length scales.
b. Asymptotic Expansion of Potential
The second preliminary step in creating the homogenization formulation for the potential problem is to write an asymptotic expansion of the field variable. With the asymptotic expansion, we approximate the complete value of the potential over all scales using an asymptotic expansion that is the addition of the potential at the macroscopic level plus peturbations due to the existence of microstructure. We first define the components of the asymptotic expansion:
The asymptotic expansion then becomes:
where the supersciprts refer to the quantities defined above and higher order terms are multiplied by the appropriate order of e. In other words, perturbations due to the first level of microstructure are multipled by the first order of e, perturbations due to the second level of microstructure are multiplied by e squared, and so on. We retain as many orders of terms in the asymptotic expansion as we have levels of microstructure in the tissue we are analyzing.
We next have to develop a gradient that accounts for the multiple levels of structure. To do this, we first recognize that the potential is now a function at several levels of structural coordinates as defined below:
In addition, the terms in the asymptotic expansion are also assumed to depend on all the structural levels as:
We know turn to the definition of the gradient. The total gradient can be taken with respect to any scale. We will take the gradient with respect in general to x0. We then use the chain rule to define the hierarchical gradient:
If we then utilize the relationship between scales in the above hierarchical gradient we obtain:
If we then apply the above hierarchical gradient with the asymptotic expansion of the hierarchical potential we obtain:
This is the final form of the hierarchical gradient potential that we will use in the potential homogenization formulation.
c. Periodicity of Potential Perturbations
The final preliminary step before we create the potential homogenization formulation is our assumption about the behavior of the perturbations terms. Before plunging into this, we need to have a conceptual idea of what constitutes hierarchical analysis. The reason that homogenization theory, or any other microstructural analysis theory, exists is that we cannot perform an analysis of all levels of scale in hierarchical materials like biological tissues in one model. This is really based on two practical limitations. First, it is not possible to image all levels of scale and therefore create a computational model that incorporates all levels of hierarchy in one dataset. Think about the analogy of computed tomography scanning. Typical patient CT scanners have a resolution of 0.3 to 0.5mm. Yet if a significant portion of the patients body is scanned, say a 30cm x 24 cm x 15 cm, we will have 600 x 480 x 300 voxels assuming isotropic voxels (equal resolution in all directions). Assuming byte data, this dataset will be 86.4 MB. A large dataset, but not unreasonable with todays computers and storage. Micro-CT scanners produce voxels typically at a resolutionof 50 microns (this is an upper bound, many such scanners can produce data at 10 micron resolution). Note that if we could even scan the same volume at a resolution of 50 microns, we would produce 10x the number of voxels in each dimenion! This means that a micro-CT data set would be 1000 times larger for the same volume than a patient CT dataset. Thus, a micro-CT dataset would be 86,400 MB or 86.4 GB. This is clearly an unreasonable dataset to handle and is larger than most computer hard drives. When we note that most biological cells have a diameter of 10 to 20 microns, requiring a 10 to 125 large dataset to resolve their placement in tissue matrices, you can see that such model creation is not feasible with current technology. Second, even if we could handle and create models from such large datasets, it would not be computationally possible to analyze such large models.
Based on these limitations, many researchers have turned to microstructural or hierarchical analysis methods. Even this case, rather than analyzing the whole hierarchy at once, we analyze pieces of the hierarchy, and assume that the pieces we analyze are representative of the complete hierarchy. This is illustrated conceptually below in Figure 2:
Figure 2: Concept of using a Representative Volume Element (RVE) to perform analysis of a hierarchical microstructured materials. The RVE should be statistically representative of the microstructure. The RVE is then analyzed under assumed boundary conditions to calculate a typical microstructural response to a global potential gradient.
From this illustration you can see that we need to perform analyses on each hierarchical level. In fact, that is the purpose of microstructural analysis to break up the modeling task into smaller analyses separated at each level rather than one large all encompassing analysis. You can also see that the analysis on the more microscopic level is performed on a representative piece of the microstructure, denoted as a Representative Volume Element, or RVE. We thus call this type of hierarchical analysis a Representative Volume Element or RVE analysis. The critical question is what boundary conditions should be assigned to the RVE when performing the analysis. We do not a priori know the RVE boundary conditions because this is really a chicken and egg question. To determine the RVE boundary conditions we would have to know results of the global analysis, but to perform the global analysis we would have to first perform the RVE analysis. Furthermore, due to the nature of the RVE model, we typically don't even know the structure of the neighboring material to hazard a guess at boundary conditions. Therefore, we must assume boundary conditions for the RVE, leading to an approximation in the analysis. The two most assumed boundary condtions are uniform values of the potential, or periodic values of the potential. In homogenization theory, the latter, periodic boundary condtions are assumed. Periodic boundary conditions state that the value of the potential on one side of the RVE must equal the value of the potential on the opposite side of the RVE. Periodic boundary conditions on Figure 3:
would require the following:
d. Homogenization Preliminaries - Summary
In summary, there are three major preliminary concepts in homogenization theory:
i. Relationship
between coordinate levels determined by e, the ratio of length scales, as: 
ii. Asymptotic
expansion of the field variable, in this case the potential as:
which
leads to a hierarchical gradient defined as:
iii. The
potential on the RVE boundary is assumed periodic, which in the 2D case leads
to:
2. Homogenization Formulation of the Potential Problem
a. Formulation of the governing Potential PDE using Homogenization Theory
To begin the homogenization formulation for the potential problem, we recall the governing PDE for the Potential problem:
To simplify our derivation without losing generality, we will assume we have a two level hierarchical structure. In this case, our hierarchical gradient becomes:
We can use the above hierarchical gradient to rewrite the governing Potential PDE as:
We require that the above PDE be satisfied if all equations associated with a given magnitude of e are satisfied independently. This gives us the following separate equations:
Note that in this case we have neglected terms of order e. This is because we generally assume that the ratio between hierarchical length scales is large, meaning that e will be very small, generally less than 0.1. Thus, we assume the influence of terms associated with e will be small. Let us examine the first PDE, that associated with e-2. We first note that we can integrate this equation over x1 to obtain the following:
where C(x0) is a constant over x1. In this case, we assume that sij cannot be zero, which is physically reasonable since the conductivity or permeability cannot be zero. Otherwise, we would have no need to solve the problem! Next, we integrate our first solution above again over x1. Noting that f is assumed to be periodic, we obtain the following:
This means that C(x0) must be zero. Using the fact with our first solution, we show that:
This result shows that f0 depends only on the macroscopic level x0. This result also means that the first terms of the equation associated with e-1 is automatically satisfied. This leaves us with with two governing equations. These equations turn out to be the governing PDE's for the Potential problem under the homogenization formulation:
We will denote these as the microscopic and macroscopic governing Potential PDE's.
b. Finite Element Approximation for the Homogenization Form of the Potential Problem
As noted above, the homogenization approximation for the hierarchical Potential problem leaves us with two governing PDE's, which we label appropriately enough the microscopic and macroscopic Potential PDE's. To begin the finite element approximation to these PDE's, we will write them in the matrix vector form we used for the finite element approximation to the single scale potential problem:
Let us begin with the microscopic equilibrium equation. To create the weak form, we multiply the microscopic PDE by a trial function or virtual potential y. We then integrate over the volume of the domain, which in this case is both the microscopic and macroscopic volume domains:
At this point, we note an important property about periodic functions. In the limit as e -> 0, we can break the total integral over the volume domains into a double integral with one integral over the microscopic volume domain divided by the total volume of the domain, followed by a second integral over the global domain. Denoting the microscopic volume as Vrve (for the representative volume element of the microscopic domain) and the macroscopic volume as Vg, we obtain:
since y is an arbitrary trial function, we will assume that it varies only over the microstructure. This will allow us to focus on the inner integral over microscopic domain. Performing integration by parts as in the previous potential finite element formulation, and assuming that no flux is applied at the microscopic level, we obtain:
We have no boundary flux term in the above equation because we assumed at the beginning of our derivation that no flux is applied on the microscopic level. The above equation indicates that the perturbation in the potential due to the microstructure F1, is a function of the global potential F0. This again sounds like the chicken and the egg problem, namely that we must know the global potential F0 to calculate F1, but we must also know F1 to calculate F0. However, if we examine the right hand side of the above equation, we will note that the gradient of the global potential is first of all constant over the microscopic domain. We also note that since the problem is linear, we can represent any arbitary global potential gradient as the superposition of a unit global potential gradients and the actual global potential gradient itself. This is just multiplying a vector by the identity matrix to obtain the vector back again. We can write this superposition as:
Noting that the three columns of the identity matrix represent three separate unit global potential gradients, we can then solve three auxillary problems where we substitute these unit global potential gradients for the actual global potential gradient. This gives us the following three auxillary problems to solve:
Auxillary
Eq. 1
Auxillary
Eq. 2
Auxillary
Eq. 3
In auxillary eqs. 1-3, the unit global potential gradient has
been taken to be negative, replacing the negative sign on the right hand side
of the equation. Since the solution to the auxillary equation is not the same
as the solution to the original global potential gradient, we denote this potential
as the characteristic potential F*. If we
multiply the right hand side out, we see that each right hand side is actually
equal to a column of the conductivity or permeability matrix:
Aux eq. 1
Aux Eq. 2
Aux Eq. 3
We can now note that since the problem is linear, we can superpose the solutions to the auxillary problem to find the solution to our original microscopic equilibrium problem in the same way that we superpose the unit global potenitals to obtain the actual global potentials -- namely by multiplying a matrix with the characteristic solutions by the original global potential:
At this point, we are ready to generate the finite element formulation for the auxillary equations. This is done by first interpolating the trial potential, the characteristic potential, and the microscopic level coordinates using shape functions:
Note that the left hand side of the equation is the same as the traditional potential problem. Therefore, would go through the exact same finite element approximations to obtain:
For the right hand side, we substitute for the trial potential gradient, and generate the appropriate matrices to obtain:
If we identity the following matrices:
We then have the finite element formulation for the local potential problem:
Once we solve this, recall that the total potential gradient is written as:
If we substitute our expression for the gradient of F1 into the above equation we obtain:
Now, let us turn our attention to the global potential equation. We have the following weak form for the global equation:
We note that we can substitute for the gradient in this equation using the above expression to obtain:
We have pulled both the global potential and the trial potential gradient out of the RVE integral because they are by definition constant over the RVE. This means that all the terms inside the RVE integral are known and we can compute these before we perform the global analysis. In fact, if we examine the terms inside the integral, we can see that if we define these as a conductivity or permeability, that the above equation has the same form has our original single scale potential equation. We thus can define the effective conductivity as:
Thus, the above equation is a precise mathematical description of a structure function relationship since it shows how an effective material propety depends on the properties of the microstructure and the architectural arrangement. We can use the above relationship in the weak form of the global potential problem to obtain:
For the right hand side, we note that y0 is only defined on the global level and we can pull this out of the integral. The current source Iv will now be an average current source at the global level. Thus, we can rewrite the above equation as:
Note that the above equation has the exact same form has the single level potential problem. We can therefore use the same finite element approximation as for the single level potential problem to obtain:
If we again identify the following matrix and vector in shorthand:
We have the following matrix vector equations for the global problem:
Finally, there is a question of localization, or how do we compute the potenital and its gradient at the microscopic level once we have calculated F1 and F0. Recall that the total potential and its gradient were approximated as:
We saw however, previously, that we can rewrite the the gradient of F1 as:
By using the Identity matrix times the gradient of F0, we can rewrite the complete potential gradient as:
By examination of the above equation, we can write the total potential as:
The above two equations represent the process of localization. Localization is calculation of quantities at the microscopic level from the macroscopic field equations and the characteristic solution of microscopic equation.
c. Summary of the Finite Element Approximation for the Homogenization Potential Analysis
In summary, there are four steps for hierarchical analysis of tissues using the homogenization method. These are in order of application:
1. Solve the local microscopic problem for as all designated RVE's to give characteristic potential solution:
where the matrices in this equation are defined above.
2. Calculate the macroscopic materials properties: this is homogenization
3. Solve the macroscopic problem using the homogenized macroscopic material properties from step 2.
where the matrix quantities are as defined above.
4. Compute the potential at the microscopic level using characteristic potential gradient and global gradient. This is called Localization.
Potential Localization:
Potential Gradient Localization:
BME 506 Home
