Stiff spring boundary conditions?

[howtophysics]

Let me see if I can line out my question a little better to hopefully get some sort of input. I am trying to understand where a type of boundary condition approx. called stiff spring BCs. I have, among a couple of other examples, an example comsol "dialysis" model that uses it. I have been searching references for quite some time on this topic and I haven't been able to find anything that gives an adequate explanation (a full one) of where this approx. comes from and "how" it works. Even people I know who work with comsol all the time have not heard of this before:

how the "dialysis" model is set up: you define three separate concentrations for 3 regions, which you solve independently - c1 (region 1) c2 (region 2) and c3 (region 3). Region 1 (convection and diffusion of a gas moving through a fluid) is linked to region 2 (diffusion of a gas moving through a membrane) which is in turn linked to region 3 (convection and diffusion of a gas moving through a fluid). They are linked through stiff spring BCs:

between region one and two (it is similar between regions 2 and 3):
(-D∇c1+c1u)⋅n= M(c2 - Kc1) --> (inward flux approx. on the c1 side)
(-Dm∇c2)⋅n= M(Kc1- c2) --> (inward flux approx. on the c2 side)
where:
*** Dm is the diffusivity of the gas in the membrane
*** D is the likewise in the fluid
*** K=c2/c1 is a partition coefficient derived from Henry's law (which I assume means that they are referring to the saturation concentrations under the conditions in question)
*** M is an arbitrary large (?) stiff spring velocity

I am trying to understand this approximation better and determine whether or not it is good to use in my case (I have two regions, a diffusive layer and a channel with convection and diffusion). ANY input that you could give at this point would be greatly appreciated (even if it's a guess or just pointing me in a certain direction given something similar that you've seen).

Questions:

1. why does this work (I can guess intuitively from playing around with it, but I would rather not - I'd like a more complete explination)?
2. where does it come from?
3. how is it derived?

Notes:

• the documentation for the model in question if you would like to check it out: http://www.engineering.uiowa.edu/~cbe_217/Misc/dialysis.pdf
• this guy seems to have used stiff spring boundary conditions in the way that I would like to: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2678100/pdf/1742-4682-6-5.pdf
• another example of someone who used them: cds.comsol.com/access/dl/papers/1530/Clark.pdf
• I did some searching, and I've found other people with similar hesitations/problems:
  http://www.comsol.se/community/forums/general/thread/6473
  http://www.comsol.nl/community/forums/general/thread/3314

 

[sshzp4]

The following might help, a little. You need the basics of ordinary diff. equations (ODE), numerical solutions to ordinary differential equations and a little experience with numerical computations to understand the following.

'Stiff boundary condition' is a term borrowed from numerical solutions to differential equations. Many (nearly most with in engineering) ODE's have smooth, continuous phase spaces (the hypothetical space inhabited by solutions) where numerical solvers don't have to make very tiny increments to their independent variable step size to come up with a high accuracy answer.

However, some problems are 'stiff', where the nature of solution changes quickly. Imagine the trajectory of a bouncing ball. The ball impacts the earth, and the impact event is a 'stiff' boundary condition. So in calculating the numerical solution to such a problem you have to ensure you can detect when the boundary (contact with earth) event takes place. To do this you use special ODE solvers called 'stiff solvers'.

I am not familiar with the physical model you describe, but I assume (until corrected), that the membrane is the stiff event in your current case. The nature of your solution changes rapidly once you exit the c1 side to enter c2 and then exit to c3.

The reason why it's called stiff spring boundary condition is purely historical as most linear ODEs were imagined using mass-spring-damper models and you would impose the boundary condition using a nonlinear spring model.

Till corrected, I will describe your diffusion equations as being similar to chemical kinetics. You have backwards and forwards diffusion which implements the changes in concentrations, that's what the RHS of both equations mean (+=forwards, and -=backwards). The LHS is simply the projection of the diffusion current in space (concentration/space) to model the number of molecules passing through the membrane (which is oriented in $\mathbf{n}$). This is just a mathematical representation of what someone thinks is going on inside the system.

 

[howtophysics]

Thanks for the reply! That does help - and I really appreciate some summary background.

What you said makes a fair bit of sense (but let me try) -- so the solution changes quickly at a boundary. I guess that this would come, in part, from the different solubilities of the gas between the membrane material and the fluid - i.e. if one can hold more gas than the other and you are trying to show how the gas moves through them, then you would need a BC and can deal with the "jumps" in your FE solution.

So I guess then that this technique would be valid (in theory) for any similar system with a "jump" (non smooth phase space) in boundary values... that is, the idea of stiff spring conditions is centered around the boundary difference. Therefore, in theory, this approximation should work/be valid (if we assume it worked here) for my very similar system that simply involves two regions with a boundary jump... (all other approx. [etc.] being as valid as they may be).

 

[sshzp4]

Sure. Let me point this out once, it is also possible to have 'stiffness' within the definition of the system itself, not merely within the boundary values of the R^3 space used to define your system (in your case it is the x, y, z coordinates introduced by the del operator); for example, your value of K could also be modeled using some other equation, which might be very stiff. But that adds a new state variable to your system as well (c1, c2, c3 and K, your phase or state space is now four dimensional). So even if your original system were non-stiff, this new stiffness in the value of K would still make the phase space stiff.

And depending upon how you chose to define your problem, sometime you can get away without considering the stiffness at all. For example, if you consider the position of the dropping ball from 10m, till, say 0.1 cm, above where it hits the ground, you can get away without considering the impact event (dd x = -g).

Finally,"I guess then that this technique...centered around the boundary difference". Yes. But I would rephrase "stiff spring conditions is centered around the boundary difference" to "stiff spring condition is centered around a rapidly changing solution" or "the boundary is modeled using a rapidly changing solution". Another rephrasing would be of "simply involves two regions with a boundary jump" to "involves two regions with a boundary (defined using the state variables)". The boundary is the jump. The boundary exists in physical space (sometimes we call this 'the dimensions occupied by the system'), the jump models what happens at the boundary in the solution space (values of c_i).

Alternatively, the opposite way of thinking about this is: if you notice a jump in the solution space, then that is where the boundary lies.

An example: Plot this function: y=x-x/(1+exp(-B*(-8+2x))) for domain x=-10 to +10; Use a large number for B (say 1e5). The boundary can be observed to be at x=4; if you had constructed an ODE using this function (it would look like dd y +p d y +q y + r =0; p,q,r are f(x)), and tried to simulate it using a non-stiff solver, you will see results different from the actual solution.


(Misc: The function y(x), is essentially the product of a straight line, y1=x, and a 'sigmoid', y2=1-1/(1+exp(..)). So y=y1*y2. The sigmoid y2 is essentially a step function, very discontinuous, which goes from 1 to 0 at x=4. y1 is a smooth line. If you look at the plot and compare it with a plot for coefficient of friction (COF), then you can see that this mathematical model reflects the physics behavior of the COF very well. So with the right choice of numbers, I could use y(x) to be mu(F) where mu is COF and F is applied tangential load. I could go on.. but I should stop)

Sid

 

[howtophysics]

Wow - again - thanks Sid! That really did help me frame the concepts quite a bit better. Are there any books on the subject you could recommend that have helped you in your understanding of these kinds of concepts (relating to FEA)?

 

[sshzp4]

I am not sure if any texts helped me out completely. Most of this just comes from mathematical modeling experience. If you are really interested and have the time to devote, try understanding linear algebra (Gilbert Strang's MIT lectures videos on youtube+his or Jim Hefferson's text on the subject) from the computational point of view. Then there's the topic of numerical simulation of ODE's, my favorite text is Kreyszig/Engineering math, and Hildebrand-Advanced calculus for engineers. Try writing your code for runge-kutta solvers/all flavors.

My most important suggestion: Play with MATLAB or your favorite scientific computing software. ALWAYS write your own code. Understand how programs handle numerical integration. Once you start visualizing the solutions in your head, you will be able to frame your own questions. And that's when life becomes good. Best of luck and enjoy the math!