Advection/Diffusion+Reversible Surface Binding

[john quinn]

All,
I have been working through the examples and am impressed with the performance of Elmer but it looks like I cannot run the numerical model I need without writing a new solver and its been 20years since I looked at Fortran or maybe I am missing something obvious that I could use!!
The model describes advection/diffusion of a species in solution through a rectangular flow cell. The majority of the species exits the flow cell but a small amount is bound reversibly to one wall of the flow cell due to affinity binding that can be implemented as a boundary condition at this wall. The partial differential equation for binding of a species to a surface is well know. Below I outline details of the model as it might help spring an idea to mind. Also where can I find all DLL files in the Tutorial regarding Electroosmotic Injection? Thoughts,
John.

Model Specification. Mass transport of the species, A, from the inlet through the flow can be modeled by the Advection/Diffusion equation and I noted that that solver is available. The flow cell is 0.03mm (height, ymax)*3.0mm(Length ,xmax) *0.5mm (Width,zmax). The width ,z, can be neglected in the simulation because of the aspect ratio. . A 32x32 elements work fine. Binding occurs only at the elements in contact with the floor of the flow cell. Diffusion in x direction is insignificant in this binding region. A simplified laminar flow model for advection with reduced dimensions will also work fine and is given as follows. I use d in the notation for the partial derivative for convenience.

dA/dt = D ( d^2A/dx^2 + dA/dt) - 4.v.y(1-y) dA/ dx
where A is the concentration of the species A at the binding surface (where AB is the concentration of the species bound at the binding surface at a given time.
D is the diffusivity of the species A
v is the average fluid velocity inside the flow cell
The concentration of A entering the flow cell is described by function G with respect to time and is the inlet boundary condition.
A(t,0,y) = G(t)
Only a small amount of A is bound relative to that leaving the flow cell by flow thus the boundary condition is.
Exit : dA(t,xmax,y)/dx = 0 (convective flow)
Roof: dA(t,x,ymax)/dy = 0 (no binding, no flux)
The analyte species, A, will be transported in the bulk via convection and will diffuse across the non-stirred boundary layer where a small fraction forms a complex with immobilized binding sites, B. It is a simple 1:1 reaction where A+B ↔AB. Thus the free analyte concentration just above the binding surface ( i.e floor of the flow cell inside non-stirred boundary layer ) is given by;
D . dA(t,x,y)/dy = dAB(t,x)/dt = kon . A(t,x,0) . (Bmax-AB(t,x)) - koff.AB(t,x)
kon is the association rate constant
koff is the dissociation rate constant
Bmax is the maximum concentration of binding sites at the binding surface
AB is the concentration of complex AB at the surface at any given time.
I need to be able to simulate the flux of A and AB as a function of time by summing over a region (Xmax-x) . Ideally the z dimension could be included to sum over a specified area relative to the inlet.

[raback]

Hi

Quickly reading through your problem setup my guess is that you might be able to do your simulation with the existing AdvectionDiffusion solver (coupled with the flow solver). It depends on the fact whether you can know [AB] a priori, or if there should be a balance for it too.

The solution would involve expressing the BC in two parts

D . dA(t,x,y)/dy = dAB(t,x)/dt = kon . A(t,x,0) . (Bmax-AB(t,x)) - koff.AB(t,x)
kon is the association rate constant
koff is the dissociation rate constant
Bmax is the maximum concentration of binding sites at the binding surface
AB is the concentration of complex AB at the surface at any given time.

Assuming you would solve for "Concentration" the following BC definitions might be quite close

Code: Select all

Mass Transfer Coefficient = Variable Time, Coordinate 1
  Real MATC "-kon * (Bmax-AB(tx(0),tx(1)))"
Concentration Flux = Variable Time, Coordinate 1
  Real MATC "koff * AB(tx(0),tx(1))"

It's advicable to introduce the complex BCs gradually starting from constant expressions to check the signs etc.

-Peter

[john quinn]

Racback,
This is a great starting point.
Many Thanks,
I will let you know how it goes,
Happy new year!

John.

[john quinn]

Using the boundary conditions as suggested should work for a simplifed case where I could use an analytical expression to determine AB as a function of time using the concentration of the advected species as an input parameter. However, to obtain the most useful information I need to solve the ODE for the formation of complex as a function of time. Is there a way to couple an ODE solver with the advection/diffusion solver ?
Are there any links where I can find open source fortran solvers ?

[raback]

Hi

It would probably be easiest to write the ODY solver using the standard PDE machinery of Elmer. Ideally the matrix equation would of course not be needed but since [AB] is present only at boundaries the additional cost is not significant. And perhaps there is some minor diffusive transport on the surface as well which would also regularize the solution.

So you could choose some time-dependent solver as a starting point and build up it to an equation for [AB] balance. A minimalistic 2nd order solver is, for example, the wave equation in tests/WaveEqu. To apply a solver at an boundary only give the corresponding boundary "Body Id" and refer to the new body id in all the sif definitions.

These minimalistic instructions don't certainly help all the way but may perhaps lead you to one possible track.

-Peter

[john quinn]

I noticed that there AdvectionReaction Solver but I cannot find the F90 for this. Is it available?

[john quinn]

It would appear that my project is almost identical to the electorchemistry problem decribed
Navier-Stokes equations with diffusion and electromigration (AdvectionDiffusion model) for convective mass transfer

http://www.elmerfem.org/wiki/index.php/Electrochemistry

Simplisticly if I could subsitiute my surface mass reaction for the surface electrochemical component then this seems a good start.
Thoughts?

[hazelsct]

Hello John,

I don't think you need Nernst-Planck. That is for the specific case where the electric field moves ions around. For your case, I think you could just set Convection = Computed, and specify the x-velocity to be 4*v*y*(1-y).

That doesn't help matters with your boundary conditions, but I think it gets the internal equation right.

-Adam

[john quinn]

In your directions below you gave some boundary conditions.
I took your advice and set the boundary analyte concnetration to constant values. I began with zero and noted that when I did this the analyte dispears at the wall and so the wall behaves as an infinite sink. I get a distinct boundary layer with an anlyte gradient developing. This in itself is useful when investigating transport limitation.
However I am have problems with my sif file accepting the boundary conditions as you defined below. Are there any instructions on the syntax you are using to declare a variable as a function of time.


Mass Transfer Coefficient = Variable Time, Coordinate 1
Real MATC "-kon * (Bmax-AB(tx(0),tx(1)))"
Concentration Flux = Variable Time, Coordinate 1
Real MATC "koff * AB(tx(0),tx(1))"

[john quinn]

Raback,

I was able to determine the concentration gradient of my analyte close to a sensing surface in my flow cell where I assume the surface is an infinite sink so I just set a BC where the analyte = 0 at this surface. I am still trying to use an analytical expression for the the accumulation at the wall as I gave you before but I cannot figure out the structure of the BC you provided. Does coordinate 1 contain the analyte concentration as you gave below?
Given that tx(0) is Time and tx(1) is coordinate 1. Is coordinate 1 the x-axis of my flow cell?
In your BC Bmax-AB has component (tx(0),tx(1)). Where can I find the definition of this structure. What is it doing?

Mass Transfer Coefficient = Variable Time, Coordinate 1
Real MATC "-kon * (Bmax-AB(tx(0),tx(1)))"

I have an analytical equation for the mass balance of complex at the wall but need it to read in the changing concentration of analyte at a small distance from the boundary. I have used the SaveLine routine to output the analyte concentration at a distance of 0.5um from the sensing surface and this is a good approximation for reality. However I am having trouble figuring out how to do this within the program so that I can enter that analyte concentration into the analytical expression with each time step!
Many thanks for the help,
John.