Understanding Diffusion Vs Convection in a Multiphasic Model

[kk1101]

Hi Dr. Ateshian,

I am sorry for another question. I created this model of an arterial wall in a fluid bath and I am looking at the movement of a neutral solute across the vessel wall.

With a concentration of 1.25mM of solute in the inner wall of the bath and an actual pressure of 13.33kPa, I prescribed a normal traction of -13.33kPa on the inner wall. I also prescribed an effective fluid pressure of 10.11kPa on the inner wall from the effective pressure equation pe = p-RT*osm*deltaC.

I set the effective solubility to 0.8 just as a mock value less than 1. However, I am seeing 0 differences when I vary the osmotic coefficient, and this was surprising to me. I expect to see that the osmotic coefficient regulates the diffusion through the wall and the net solute flux follows the form:

Jnet_solute = Jdiffusion + Jadvection. Is my model correctly showing that advection is the sole contributor to movement of solute across the material or have I set up a boundary condition incorrectly?

Thank you.

[ateshian]

Hi,

The governing equations for a multiphasic material can be found in the Theory Manual. To relate the governing equations to the finite element model, please keep in mind the following: The nodal degrees of freedom in the finite element code are the solid displacement u, the effective fluid pressure (p-tilde, as shown in eq.(9), which depends on the osmotic coefficient and actual solute concentrations c), and the effective solute concentrations (c-tilde as shown in eq.(10), which depends on c and the partition coefficient). When the solute is electrically neutral, its partition coefficient is equal to its solubility. The mixture stress is given in eq.(1), it uses the actual fluid pressure p.

When you prescribe boundary conditions on pressure and concentration on a multiphasic material, you are prescribing p-tilde and c-tilde. When you prescribe a mixture traction, it corresponds to the traction of the mixture stress in eq.(1), which uses the actual pressure p and calculates the effective traction from the stress produced by the state of strain in the solid matrix (evaluated from the gradient of the displacement u).

The molar flux of the solute is given by eq.(13). This formula uses the gradient of the effective solute concentration (c-tilde) for the diffusion component, and the gradient of the effective fluid pressure (p-tilde) for the advection component.

I encourage you to look at these formulas and substitute into them the assumptions of your model to see if you observe the effect(s) that you expect.

To help you figure out the boundary conditions, please see the attached document ( kk1101.pdf ). In summary, the fluid in the external bath may also be treated as a multiphasic mixture: There is no solid matrix in that bath, so the normal component of the mixture traction is -p*. The solubility of a solute in a bath is always equal to 1, and if the solute is neutral it follows that its concentration c* is equal to its effective concentration c-tilde* (since one assumes that the electric potential in the external bath is zero).

Best,

Gerard

[kk1101]

Hi Dr. Ateshian,

Thank you for the kind response. Forgive me for following up with another question. I went through the theory and the attached PDF. Given these equations, is my understanding of these three boundary conditions correct in my attached PDF (FEBio_BC_Math)?

I am not completely understanding equation 2 in your PDF. Shouldn't the matrix stress (sigma^e) balance the fluid pressure in the tissue [sigma+pI = sigma^e]? This is because we are saying the fluid pressure in the bath is 0. Therefore p-tilda != p^*-tilda?

Thank you,
Keshav

[kk1101]

Hi,

I am sorry for yet another question. When testing the bounds of my model, I am noticing that if I set deff = d0, I get this weird behavior where the effective diffusivity increases for a period and then decreases. Is this because of a problem in the way I have meshed the model? I have tried a few other meshes and biased it closer to the direction of the prescribed concentration in both the bath and aortic wall, but I am struggling to get the model to converge in any case when deff = d0. I have attached my model.

The problem appears independent of the loads I have prescribed in both the lumen and the bath, and it seems independent of the axial boundary condition on the whole structure. I purposely have set minimal boundary conditions on the bath - I am perfectly fine accepting how much ever fluid flows into the bath. I have attached a photo of the "weird" bath behavior. I am so sorry; I genuinely cannot understand why it is behaving this way.

Thank you so much!

[ateshian]

Hi,

When you say "the effective diffusivity increases for a period and then decreases" how are you assessing the diffusivity? The diffusivity models you are using are constant diffusivity models, so in a strict sense diffusivity cannot change in your model. Do you mean that the solute flux increases then decreases? I don't see this behavior in the file you attached (I had to remove the PLOT_MUST_POINTS request otherwise I couldn't see the progress of the analysis). Also, I beg you to please use a biased mesh in the domain of Material2, it is impossible to get good diffusion results with such a coarse and uniform mesh.

Best,

Gerard

[kk1101]

Hi Gerard,

Thank you as always for the kind response. I apologize for the lack of clarity in my question. Attached is a pdf of a powerpoint which describes what I am seeing in the previously attached model and the changes I made (including the biased mesh). Are these changes OK or am I artificially solving the problem?

Thank you!

[kk1101]

Hi Dr. Ateshian,

I know you are incredibly busy, and I do apologize. I was wondering if you had any more thoughts on this.

I am sorry.

Thank you.

[ateshian]

Hi,

I looked at the PDF document you attached but I am not able to figure out the question without having to read our prior exchanges and look at your models to understand what you did. I wish I could respond quickly, but it would be helpful if you could outline (a) what are you trying to do, (b) how are you setting up your model (geometry, initial conditions, boundary conditions, loads), (c) what is the exact issue you are encountering that you would like me to address. I apologize for making you do all this work, but I am not able to respond quickly to forum questions unless the problem is presented clearly to me.

Best,

Gerard

[kk1101]

Hi Gerard,

I am so sorry; it was not my intention at all to be disrespectful of your time. I updated the powerpoint with all of the requisite information. I attached both the powerpoint and the model to this post. We are almost done with the model, but I want to be sure that the results we are presenting are correct. Thank you for taking the time over the last several months to answer my numerous questions; I am very grateful. I am sorry if not all of my questions were valid, or if I should have worked harder before asking.

Please let me know if there is any more information I should provide.

The Problem:

When the bath modulus is set to E = 0.001kPa, we see an increase and then a decrease in effective concentration in the bath as Deff approaches D0. I am solving this behavior by increasing the bath elastic modulus to E= 0.01kPa. Is this OK?

[ateshian]

Hi,

Thanks for sharing these details. Since I was still a little concerned about your mesh, I decided to create my own model based on the information provided in your PDF and FEBio files. Since I wanted to include a lot more elements in the radial direction, I decided to create a wedge model and use axisymmetry to constrain the mesh size ( Model_kk1101.fs2 ). I also modified your boundary conditions as follows: Instead of prescribing zero effective fluid pressure at the interface between the vessel and bath, I placed the zero effective fluid pressure and zero effective solute concentration at the edge of the bath (the outer boundary). So I did not prescribe zero effective fluid pressure on the top and bottom surfaces of the bath as you did. I also prescribed the fluid pressure and solute concentrations on the inner lumen of the blood vessel as step increases, instead of ramping them up. The results I got were entirely what one would expect: At steady state the final solute concentration decreased nearly linearly from the inner lumen to the outer boundary of the blood vessel, and remained essentially equal to zero in the entire bath.

Please take a look at this model and see if the variations that caused some concerns in your earlier posts persist after using this model.

Best,

Gerard

[kk1101]

Hi Dr. Ateshian,

Thank you so much for taking the time on this. Why and how did you choose to set the solver to "Full Newton"? I thought Broyden would be a good choice for this model.

Thank you.

[ateshian]

Hi,

In my experience, multiphasic problems don't converge well with quasi-Newton solvers. So I always use the full Newton solver by default. Now I expect that there may be cases where Broyden might work just fine, but I almost never start with it when solving multiphasic problems.

Best,

Gerard