Note: This discussion is about an older version of the COMSOL Multiphysics® software. The information provided may be out of date.

Discussion Closed This discussion was created more than 6 months ago and has been closed. To start a new discussion with a link back to this one, click here.

Rotation of a solid body under fluidic vortex

Please login with a confirmed email address before reporting spam

Dear COMSOL Community,

I am trying to simulate the rotation of a biological cell in a microfluidic vortex using Fluid Surface Interface (FSI). I have modeled my simulation based off the FSI example for comsol, however the model will not converge. The initial values can be computed, however the model reaches a convergence error before the first time step is realized.

I have attempted to constrain the input to have a smooth transition and to be parabolic (see functions and parameters) with little success.

Ideally, the circle (cell) would be contained to only move rotationally, without translation (even though I cant even get to this step).

Any suggestions would be appreciated and the model is attached.


14 Replies Last Post Jan 25, 2012, 2:22 a.m. EST
Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 22, 2012, 8:38 a.m. EST
Hi

in V4 you can make a "step" easier by using the built-in step function, then you can plot it too to see if it's as you want, as I do not really understand your Heavyside function there

If you want to have a closer to laminar flow, on the input side use a x-velocity of the type

"Vx_max*4*s(1-s)" or "Vx_ave*6*s(1-s)"

this will give you a parabolic shape, if you multiply with the step function you can get it to grow nicely too, but this is perhaps more appropriate for higher fluid velocities

In your case you could try to remove the segregated steps, (push all into the first step and delete the two others), if you have enough RAM it will probably solve better

In a true example you should probably (once later) add gravity loads, then do not forget to add output gravity load pressure depending on hight "Y" to simulate a large container, else you will have excessive return flow in the output region

To see the radial pressure acting on your cell, plot a polar plot of "p" with atan2(y,x) as Theta angle expression. You might need to remove the average pressure value to get a more representative view

To see how the mesh distorts, use "wireframe" in your surface plot

you can see your particle doing a 1/4 circle in 4 seconds with a 200 um/s flow velocity, you might need to consider solver remeshing and "hyperelastic" ALE mesh smoothing as you get many "inverted elements"

--
Good luck
Ivar
Hi in V4 you can make a "step" easier by using the built-in step function, then you can plot it too to see if it's as you want, as I do not really understand your Heavyside function there If you want to have a closer to laminar flow, on the input side use a x-velocity of the type "Vx_max*4*s(1-s)" or "Vx_ave*6*s(1-s)" this will give you a parabolic shape, if you multiply with the step function you can get it to grow nicely too, but this is perhaps more appropriate for higher fluid velocities In your case you could try to remove the segregated steps, (push all into the first step and delete the two others), if you have enough RAM it will probably solve better In a true example you should probably (once later) add gravity loads, then do not forget to add output gravity load pressure depending on hight "Y" to simulate a large container, else you will have excessive return flow in the output region To see the radial pressure acting on your cell, plot a polar plot of "p" with atan2(y,x) as Theta angle expression. You might need to remove the average pressure value to get a more representative view To see how the mesh distorts, use "wireframe" in your surface plot you can see your particle doing a 1/4 circle in 4 seconds with a 200 um/s flow velocity, you might need to consider solver remeshing and "hyperelastic" ALE mesh smoothing as you get many "inverted elements" -- Good luck Ivar

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 1:20 a.m. EST

Hi

in V4 you can make a "step" easier by using the built-in step function, then you can plot it too to see if it's as you want, as I do not really understand your Heavyside function there

If you want to have a closer to laminar flow, on the input side use a x-velocity of the type

"Vx_max*4*s(1-s)" or "Vx_ave*6*s(1-s)"

this will give you a parabolic shape, if you multiply with the step function you can get it to grow nicely too, but this is perhaps more appropriate for higher fluid velocities

In your case you could try to remove the segregated steps, (push all into the first step and delete the two others), if you have enough RAM it will probably solve better

In a true example you should probably (once later) add gravity loads, then do not forget to add output gravity load pressure depending on hight "Y" to simulate a large container, else you will have excessive return flow in the output region

To see the radial pressure acting on your cell, plot a polar plot of "p" with atan2(y,x) as Theta angle expression. You might need to remove the average pressure value to get a more representative view

To see how the mesh distorts, use "wireframe" in your surface plot

you can see your particle doing a 1/4 circle in 4 seconds with a 200 um/s flow velocity, you might need to consider solver remeshing and "hyperelastic" ALE mesh smoothing as you get many "inverted elements"

--
Good luck
Ivar


Ivar,

Thank you so much, my model is now converging properly! A few things:

1) Is there a way of plotting "rotation rate" of the cell vs. time? I may want to modulate the input velocity over time to see how it effects the rotation rate and this would be helpful.

2) Can you explain how (and what) hyperelastic ALE meshing smoothing works? And how to do it of course.

3) Is there a way to constrain the cell about a point so that it can not translate in x/y but only rotate? I run into convergence problems at higher velocities and or time locations because the cell has started translating.

Thanks in advance,
Jakrey
[QUOTE] Hi in V4 you can make a "step" easier by using the built-in step function, then you can plot it too to see if it's as you want, as I do not really understand your Heavyside function there If you want to have a closer to laminar flow, on the input side use a x-velocity of the type "Vx_max*4*s(1-s)" or "Vx_ave*6*s(1-s)" this will give you a parabolic shape, if you multiply with the step function you can get it to grow nicely too, but this is perhaps more appropriate for higher fluid velocities In your case you could try to remove the segregated steps, (push all into the first step and delete the two others), if you have enough RAM it will probably solve better In a true example you should probably (once later) add gravity loads, then do not forget to add output gravity load pressure depending on hight "Y" to simulate a large container, else you will have excessive return flow in the output region To see the radial pressure acting on your cell, plot a polar plot of "p" with atan2(y,x) as Theta angle expression. You might need to remove the average pressure value to get a more representative view To see how the mesh distorts, use "wireframe" in your surface plot you can see your particle doing a 1/4 circle in 4 seconds with a 200 um/s flow velocity, you might need to consider solver remeshing and "hyperelastic" ALE mesh smoothing as you get many "inverted elements" -- Good luck Ivar [/QUOTE] Ivar, Thank you so much, my model is now converging properly! A few things: 1) Is there a way of plotting "rotation rate" of the cell vs. time? I may want to modulate the input velocity over time to see how it effects the rotation rate and this would be helpful. 2) Can you explain how (and what) hyperelastic ALE meshing smoothing works? And how to do it of course. 3) Is there a way to constrain the cell about a point so that it can not translate in x/y but only rotate? I run into convergence problems at higher velocities and or time locations because the cell has started translating. Thanks in advance, Jakrey

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 2:09 a.m. EST
Hi

2) there are 3 ALE smoothing approaches, Lagrange, Winslow and hyperelastic, check the doc, you select the one in one of the node drop-down fields, but thee are at different places depending on which module you call

3) to constrain the cell, you would need to add a weak constraints that forces the average position of the cell border to stay at a given location (X0,Y0) (two weak equations one for X (u) and one for Y (v)

1) Plotting rotations are slightly trickier, as you need to extract the ange. For small rotations you have the curlU variables or via an atan2(V,U) where V,U are some average value. But as your cell is "rigid", you could take a closer look at the rigid body BC there you have as output a rotation you can follow, also valid for large angles, as its based on quaternions development, that removes any singularity

--
Good luck
Ivar
Hi 2) there are 3 ALE smoothing approaches, Lagrange, Winslow and hyperelastic, check the doc, you select the one in one of the node drop-down fields, but thee are at different places depending on which module you call 3) to constrain the cell, you would need to add a weak constraints that forces the average position of the cell border to stay at a given location (X0,Y0) (two weak equations one for X (u) and one for Y (v) 1) Plotting rotations are slightly trickier, as you need to extract the ange. For small rotations you have the curlU variables or via an atan2(V,U) where V,U are some average value. But as your cell is "rigid", you could take a closer look at the rigid body BC there you have as output a rotation you can follow, also valid for large angles, as its based on quaternions development, that removes any singularity -- Good luck Ivar

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 4:57 a.m. EST
I am trying to become more versed in using weak constraints. Do you have any recommendations for resources/models regarding this?

In any case, what would the expression(s) be for the two weak constraints?

I also see that there is a pointwise constraint, would this be of any use?

Thank you,
Jakrey
I am trying to become more versed in using weak constraints. Do you have any recommendations for resources/models regarding this? In any case, what would the expression(s) be for the two weak constraints? I also see that there is a pointwise constraint, would this be of any use? Thank you, Jakrey

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 10:28 a.m. EST
Hi

here you have 2 ways that seem to work ;) its 4.2a (166)

enable or disable the two global weak constraints, and or the rigid connector , or disable all 3, solve and update the derived values (note an error will be generated if the rigid connector are not defined because a variable is floating

--
Good luck
Ivar
Hi here you have 2 ways that seem to work ;) its 4.2a (166) enable or disable the two global weak constraints, and or the rigid connector , or disable all 3, solve and update the derived values (note an error will be generated if the rigid connector are not defined because a variable is floating -- Good luck Ivar


Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 2:40 p.m. EST
Ivar,

Thanks once again for the prompt replies. I am running 4.2 and cant seem to open the 4.2a file? Can you please resave and attach the file so I can access it?

(I should be receiving a trial of 4.2a here any day but currently I am without).

Thanks again,
Jakrey
Ivar, Thanks once again for the prompt replies. I am running 4.2 and cant seem to open the 4.2a file? Can you please resave and attach the file so I can access it? (I should be receiving a trial of 4.2a here any day but currently I am without). Thanks again, Jakrey

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 4:21 p.m. EST
Hi

sorry I have only the latest version, I'll try to add a few snapshots when I'm back by my COMSOl WS tomorrow

--
Good luck
Ivar
Hi sorry I have only the latest version, I'll try to add a few snapshots when I'm back by my COMSOl WS tomorrow -- Good luck Ivar

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 23, 2012, 5:47 p.m. EST
That would be greatly appreciated.

Awaiting your reply,
Jakrey
That would be greatly appreciated. Awaiting your reply, Jakrey

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 24, 2012, 1:54 a.m. EST
Hi

start 2D Solid + ALE physics + time solver
draw 2 concentric circles of a 1:2 diameter ratio (use other default settings not very important)
To easily see the mesh deformation, make 2 materials, one for each domain, i.e. E=2E10], nu=0.3, rho = 7500 for the inner circle, and E=1, nu=0, rho = 1 for the outer annulus
In solid (of both domains) fix the largest outer diameter, apply a point load on one of the 4 "corners points" of the inner circle, give it a tangential Point force to force a rotation and displacement motion of the inner circle.
then define a Rigid Connector BC on the inner circle perimeter and select Prescribed Displacement (of the rotation point, this is missing in the text ;) and select both X=0, Y=0
Define the ALE on the external annulus, Define a Free Deformation domain condition on this annulus
Select the inner circle periphery and apply a Prescribed Mesh Displacement BC of dx=u and dy=v the imposed displacement from the solid physics

Then define an AVERAGE operator aveop1() on the spatial frame, select the inner circle periphery
Define two variables U0 = aveop1(u) and V0 = aveop1(v), these two values represent the CoG of the circle and by plotting these over time you see the displacement of the centre of the inner circle, they will also serve later for the Global weak constraint

Solve with default mesh (note the Include geometry non-linearities is checked by default because of the ALE

Select the Result default "Stress"2D plot Plot, replace mises by disp for the displacement,. and select WIREFRAME color and style
Now that we have solved once and we have plots and results, one can rapidly return to the solver node and check the Plot While Solving to observe the phenomena live for next solving sequence
Add a 1D line plot, with a Global Plot and define U0, V0 the two CoG coordinates, this gives the displacement, which is "0" with the rigid connector
Add a new 1D plot with a Global plot and type: solid.thz_rig1 and 0.5*aveop1(solid.curlUZ)
the former is the Rigid connector angle (valid also for large displacements), the latter is the small angle rotation from linear theory, (the antisymmetric displacement tensor components)

Resolve to see the mesh curl

Now disable the Rigid connector and resolve to see the part move and rotate freely (the plot of the RC BC lags an error just ignore it with an OK

Ready for the weak constraints:

Define in the solid physics a Global Constrain and type U0, then duplicate the constraint and type V0
Solve and observe that you get the same result as for the Rigid connector (RC-BC).
Finally you can disable one of the two and see how you constrain the centre of the inner circle along X or Y
Observe the coG displacement on the 1D plot

Hope I got everything here ;)

--
Good luck
Ivar
Hi start 2D Solid + ALE physics + time solver draw 2 concentric circles of a 1:2 diameter ratio (use other default settings not very important) To easily see the mesh deformation, make 2 materials, one for each domain, i.e. E=2E10], nu=0.3, rho = 7500 for the inner circle, and E=1, nu=0, rho = 1 for the outer annulus In solid (of both domains) fix the largest outer diameter, apply a point load on one of the 4 "corners points" of the inner circle, give it a tangential Point force to force a rotation and displacement motion of the inner circle. then define a Rigid Connector BC on the inner circle perimeter and select Prescribed Displacement (of the rotation point, this is missing in the text ;) and select both X=0, Y=0 Define the ALE on the external annulus, Define a Free Deformation domain condition on this annulus Select the inner circle periphery and apply a Prescribed Mesh Displacement BC of dx=u and dy=v the imposed displacement from the solid physics Then define an AVERAGE operator aveop1() on the spatial frame, select the inner circle periphery Define two variables U0 = aveop1(u) and V0 = aveop1(v), these two values represent the CoG of the circle and by plotting these over time you see the displacement of the centre of the inner circle, they will also serve later for the Global weak constraint Solve with default mesh (note the Include geometry non-linearities is checked by default because of the ALE Select the Result default "Stress"2D plot Plot, replace mises by disp for the displacement,. and select WIREFRAME color and style Now that we have solved once and we have plots and results, one can rapidly return to the solver node and check the Plot While Solving to observe the phenomena live for next solving sequence Add a 1D line plot, with a Global Plot and define U0, V0 the two CoG coordinates, this gives the displacement, which is "0" with the rigid connector Add a new 1D plot with a Global plot and type: solid.thz_rig1 and 0.5*aveop1(solid.curlUZ) the former is the Rigid connector angle (valid also for large displacements), the latter is the small angle rotation from linear theory, (the antisymmetric displacement tensor components) Resolve to see the mesh curl Now disable the Rigid connector and resolve to see the part move and rotate freely (the plot of the RC BC lags an error just ignore it with an OK Ready for the weak constraints: Define in the solid physics a Global Constrain and type U0, then duplicate the constraint and type V0 Solve and observe that you get the same result as for the Rigid connector (RC-BC). Finally you can disable one of the two and see how you constrain the centre of the inner circle along X or Y Observe the coG displacement on the 1D plot Hope I got everything here ;) -- Good luck Ivar

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 24, 2012, 5:10 a.m. EST
Ivar,

Thanks again for this detailed walkthrough! I believe I have done everything according to your explanation, however I am not getting the same results from disabling the Rigid Connector vs. the two global constraints. It might just be the late night, and the need for rest.

Anyways, I have attached my model if you can have a look.

Thanks again,
Jakrey
Ivar, Thanks again for this detailed walkthrough! I believe I have done everything according to your explanation, however I am not getting the same results from disabling the Rigid Connector vs. the two global constraints. It might just be the late night, and the need for rest. Anyways, I have attached my model if you can have a look. Thanks again, Jakrey


Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 24, 2012, 6:27 a.m. EST
Hi

first of all invert your two material selections: stiff material small circle, soft material external

Then for some reason your global Constraints do not work !? just as if U0 and V0 are not calculated ?, strange

now it could also come from a version issue, I do not like to run older versions, I remake them

Try to write out the equation "aveop1(u)" fully in the Global Constraint and do not refer to U0 (normally I would expect it to be the same never knows)

--
Good luck
Ivar
Hi first of all invert your two material selections: stiff material small circle, soft material external Then for some reason your global Constraints do not work !? just as if U0 and V0 are not calculated ?, strange now it could also come from a version issue, I do not like to run older versions, I remake them Try to write out the equation "aveop1(u)" fully in the Global Constraint and do not refer to U0 (normally I would expect it to be the same never knows) -- Good luck Ivar

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 24, 2012, 2:24 p.m. EST
That seemed to work... weird?

Now I will try to implement this idea into my original model :) Thanks for the example model. Ill report back when I have worked on it more.
That seemed to work... weird? Now I will try to implement this idea into my original model :) Thanks for the example model. Ill report back when I have worked on it more.

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 24, 2012, 5:59 p.m. EST
Ivar,

I noticed that the fsi coordinate variables were defined as up(x), up (y) . If you change the expression to mean(x), mean (y) it seems to constrain the cell to only rotation however the fluid streamlines aren't correct.

What is the effect of changing this expression on the physics?
Ivar, I noticed that the fsi coordinate variables were defined as up(x), up (y) . If you change the expression to mean(x), mean (y) it seems to constrain the cell to only rotation however the fluid streamlines aren't correct. What is the effect of changing this expression on the physics?

Ivar KJELBERG COMSOL Multiphysics(r) fan, retired, former "Senior Expert" at CSEM SA (CH)

Please login with a confirmed email address before reporting spam

Posted: 1 decade ago Jan 25, 2012, 2:22 a.m. EST
Hi

My guess haven't checked (one would have to dig into the "equation views" and there perhaps) is that the "up()" are selecting the fluid side (which for me is correct) while the "mean()" takes the average of both Boundary sides, and then probably the solid side is "0" at best hence you will get something like 1/2 off expected response.

In FSI several boundaries are de-doubled so you do not have continuity across them, therefore it's important to distinguish on which side you consider the variables. i.e. U the fluid velocity is normally only defined on the fluid side

--
Good luck
Ivar
Hi My guess haven't checked (one would have to dig into the "equation views" and there perhaps) is that the "up()" are selecting the fluid side (which for me is correct) while the "mean()" takes the average of both Boundary sides, and then probably the solid side is "0" at best hence you will get something like 1/2 off expected response. In FSI several boundaries are de-doubled so you do not have continuity across them, therefore it's important to distinguish on which side you consider the variables. i.e. U the fluid velocity is normally only defined on the fluid side -- Good luck Ivar

Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.