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.
Convergence problems with non-newtonian power law fluid
Posted Jan 1, 2013, 10:51 p.m. EST Fluid & Heat, Microfluidics, Studies & Solvers Version 4.2, Version 4.3, Version 4.3a 11 Replies
Please login with a confirmed email address before reporting spam
Hi,
I'm modeling fluid flow through a microfluidic bifurcating network. Among other things, I am calculating the volumetric flow rate for a given applied pressure and comparing it to experiment. The problem I am having is that when I am using a Non-newtonian power law fluid the solution will not converge. I have run the same model/mesh/solver with a newtonian fluid (glycerol) and the solution converges. I have tested smaller models (1 inlet; 2 or 4 outlets) with the power law fluid and the solution converges, but when I do the full scale (64 outlets), it will not converge.
Boundary conditions
This model has one inlet. I set it to either Pressure (no viscous stress) or Laminar inflow (pressure). There are 64 outlets. They are all set to Pressure (no viscous stress) of 0. The walls are no slip. I mostly keep both inlet and outlets at Pressure (nvs).
Solver
The default solver is an iterative solver with multigrid. I got a warning that I needed to increase the "factor in error estimate." So I changed it from the default value of 20 to 40, but this had no effect. I also got another warning that there was an "ill conditioned pre-conditioner." I'm vaguely aware of what this means in a linear-algebra context, but I don't know how I can fix this in COMSOL.
Under Study > Solver config > Solver 1 > Stationary Solver > Fully coupled, I changed the "Nonlinear method" from the defuault of "Linear (Netwon)" to "Highly non-linear," but the solution still did not converge.
I noticed that in the model library there is a model that uses a non-newtonian Carreau model fluid.
www.comsol.com/showroom/gallery/171/
In this model, the solver is a Direct solver - PARDISO. I tried this solver on my model, and after ~8 hours the solver appeared to make little progress.
Mesh
The mesh I'm using is a free tetrahedral mesh. I vary the mesh fineness by a factor I call "C", where the maximum element size is (the width of the channel)/C. In this model I uploaded, C=5 (coarse), but I typically have C=7 or 8 (finer). I'm wondering if my meshing is part of the problem. The mesh elements are all about the same size—they are not finer closer to the walls of the channels. Though this may be a problem for the non-newtonian fluid, it seems to be OK when I am using glycerol as the fluid.
I'm also wondering if I should do some sort of swept mesh to create a hexahedral mesh since my geometry is composed entirely of rectangular prisms.
I'm hoping to gain a better intuition of what's going on and how to work this out.
I would be very grateful for any tips/insights into this problem, thanks!
I'm modeling fluid flow through a microfluidic bifurcating network. Among other things, I am calculating the volumetric flow rate for a given applied pressure and comparing it to experiment. The problem I am having is that when I am using a Non-newtonian power law fluid the solution will not converge. I have run the same model/mesh/solver with a newtonian fluid (glycerol) and the solution converges. I have tested smaller models (1 inlet; 2 or 4 outlets) with the power law fluid and the solution converges, but when I do the full scale (64 outlets), it will not converge.
Boundary conditions
This model has one inlet. I set it to either Pressure (no viscous stress) or Laminar inflow (pressure). There are 64 outlets. They are all set to Pressure (no viscous stress) of 0. The walls are no slip. I mostly keep both inlet and outlets at Pressure (nvs).
Solver
The default solver is an iterative solver with multigrid. I got a warning that I needed to increase the "factor in error estimate." So I changed it from the default value of 20 to 40, but this had no effect. I also got another warning that there was an "ill conditioned pre-conditioner." I'm vaguely aware of what this means in a linear-algebra context, but I don't know how I can fix this in COMSOL.
Under Study > Solver config > Solver 1 > Stationary Solver > Fully coupled, I changed the "Nonlinear method" from the defuault of "Linear (Netwon)" to "Highly non-linear," but the solution still did not converge.
I noticed that in the model library there is a model that uses a non-newtonian Carreau model fluid.
www.comsol.com/showroom/gallery/171/
In this model, the solver is a Direct solver - PARDISO. I tried this solver on my model, and after ~8 hours the solver appeared to make little progress.
Mesh
The mesh I'm using is a free tetrahedral mesh. I vary the mesh fineness by a factor I call "C", where the maximum element size is (the width of the channel)/C. In this model I uploaded, C=5 (coarse), but I typically have C=7 or 8 (finer). I'm wondering if my meshing is part of the problem. The mesh elements are all about the same size—they are not finer closer to the walls of the channels. Though this may be a problem for the non-newtonian fluid, it seems to be OK when I am using glycerol as the fluid.
I'm also wondering if I should do some sort of swept mesh to create a hexahedral mesh since my geometry is composed entirely of rectangular prisms.
I'm hoping to gain a better intuition of what's going on and how to work this out.
I would be very grateful for any tips/insights into this problem, thanks!
Attachments:
11 Replies Last Post Jan 4, 2013, 12:59 a.m. EST
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
are you using a boundary mesh along the no-slip walls ? this would help if not used so far. your model is so nicely symmetric that you should be able to cut it in two, perhaps even more but then you would remove any possible variations on the output flow.
if you use laminar input flow pressure driven, be avare that normally (was so for earlier versions) the pressure includes the laminar inflow length, so the true pressure drop must be measred as a differential from your inlet (and not the total pressureat input) to your outlet
--
Good luck
Ivar
are you using a boundary mesh along the no-slip walls ? this would help if not used so far. your model is so nicely symmetric that you should be able to cut it in two, perhaps even more but then you would remove any possible variations on the output flow.
if you use laminar input flow pressure driven, be avare that normally (was so for earlier versions) the pressure includes the laminar inflow length, so the true pressure drop must be measred as a differential from your inlet (and not the total pressureat input) to your outlet
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Steve,
I agree with Ivar that you should exploit the half symmetry for computational efficiency.
Also, it is true that swept meshing of regular domains is a more efficient way of meshing (ie less elements per domain for same accuracy). For your case, to sweep multiple connected sections you will need to do some partitioning. See attached mph file for an example of how to partition and define a Mesh Control face to hide the partition from the physics set-up. Also, see Mesh 2 for an example of how to use Mapped and Swept mesh techniques.
These techniques may or may not help with convergence. I'd try adding a boundary layer mesh to your current set-up as Ivar suggested. But, cut the model in half first to save memory and time!
Regards,
Josh Thomas
AltaSim Technologies
I agree with Ivar that you should exploit the half symmetry for computational efficiency.
Also, it is true that swept meshing of regular domains is a more efficient way of meshing (ie less elements per domain for same accuracy). For your case, to sweep multiple connected sections you will need to do some partitioning. See attached mph file for an example of how to partition and define a Mesh Control face to hide the partition from the physics set-up. Also, see Mesh 2 for an example of how to use Mapped and Swept mesh techniques.
These techniques may or may not help with convergence. I'd try adding a boundary layer mesh to your current set-up as Ivar suggested. But, cut the model in half first to save memory and time!
Regards,
Josh Thomas
AltaSim Technologies
Attachments:
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Thanks for the reply, Ivar.
I am currently trying a boundary layer in the mesh.
I was aware that the laminar inflow measures the specified pressure from the beginning of the defined entrances length, but thanks for pointing that out as it's easy to miss: it threw me off the first time I used laminar inflow.
Do you have any insights into the solvers, i.e. direct vs. iterative? It seems like there's so many options to choose, but maybe the default is best?
-Steve
I am currently trying a boundary layer in the mesh.
I was aware that the laminar inflow measures the specified pressure from the beginning of the defined entrances length, but thanks for pointing that out as it's easy to miss: it threw me off the first time I used laminar inflow.
Do you have any insights into the solvers, i.e. direct vs. iterative? It seems like there's so many options to choose, but maybe the default is best?
-Steve
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Josh,
I agree cutting the model in half is a good idea for saving memory/computation time. The reason I did not do that, is that initially I was intending on creating assymetric defects in the network.
In the real world, my goal is to have perfectly even flow through each output. I was planning on using COMSOL to calculate the effect differences in the geometry would have on the even-ness of flow. However, for the "control" simulation where the network is perfectly symmetric, halfing the model is a good idea.
I agree cutting the model in half is a good idea for saving memory/computation time. The reason I did not do that, is that initially I was intending on creating assymetric defects in the network.
In the real world, my goal is to have perfectly even flow through each output. I was planning on using COMSOL to calculate the effect differences in the geometry would have on the even-ness of flow. However, for the "control" simulation where the network is perfectly symmetric, halfing the model is a good idea.
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
I'm sure a boundary mesh will help, but OK it all depends on how fine youre mesh was before. And for the half symmetry, I'm rater sure 1/2 is OK, bt I'm far from sure that your current model will be fully symmetric (in flow) through the 32 outlets on one half ;) At least not for any laminar flow rate
I'm, not that a CFD specialist to play a lot with the solvers for that physics, so I tend to use the defult choice (I almost alaways start with the default choice, then if it takes too long or does not converge I start to look at the solver tree, and check if it makes sen, this is partiuclarly worth for multiple physics, when solving for one dependent variable before another has some importance, one have often a driving and a "slave" variable, then I tend to solve for the driver first when COMSOL set's up segregated sequence, except is the model is rather liht and I'm sure to enter my RAM size, then sometimes I try a "direct" one
--
Good luck
Ivar
I'm sure a boundary mesh will help, but OK it all depends on how fine youre mesh was before. And for the half symmetry, I'm rater sure 1/2 is OK, bt I'm far from sure that your current model will be fully symmetric (in flow) through the 32 outlets on one half ;) At least not for any laminar flow rate
I'm, not that a CFD specialist to play a lot with the solvers for that physics, so I tend to use the defult choice (I almost alaways start with the default choice, then if it takes too long or does not converge I start to look at the solver tree, and check if it makes sen, this is partiuclarly worth for multiple physics, when solving for one dependent variable before another has some importance, one have often a driving and a "slave" variable, then I tend to solve for the driver first when COMSOL set's up segregated sequence, except is the model is rather liht and I'm sure to enter my RAM size, then sometimes I try a "direct" one
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi Steve,
The convergence problem is most probably due to the power-law viscosity. With the constants you provided the viscosity is infinite at zero shear flow rate. That results in solution matrix ill-conditioning and high viscosity gradients in the channels. A finer mesh and/or a boundary layer mesh will improve the solution resolution around the walls and probably help with convergence, however, it’s better to avoid the infinite viscosity altogether. I recommend that you switch to a Carreau model, or any other user-defined viscosity model, that involves realistic upper and lower bounds on viscosity.
Nagi Elabbasi
Veryst Engineering
The convergence problem is most probably due to the power-law viscosity. With the constants you provided the viscosity is infinite at zero shear flow rate. That results in solution matrix ill-conditioning and high viscosity gradients in the channels. A finer mesh and/or a boundary layer mesh will improve the solution resolution around the walls and probably help with convergence, however, it’s better to avoid the infinite viscosity altogether. I recommend that you switch to a Carreau model, or any other user-defined viscosity model, that involves realistic upper and lower bounds on viscosity.
Nagi Elabbasi
Veryst Engineering
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Josh,
I took your advice and rebuilt the network with the mirror symmetry (see attached). I added a boundary layer as Ivar suggested. Thanks for sending me the example swept mesh. Is the partitioning necessary because the model does not have a constant depth, and beginning and end faces of the mesh sweep must match? When meshing the new model I opted to stick with the free tetrahedral mesh because partitioning and sweeping was too tedious. The tedium was due to an...
Unrelated problem: When I am selecting domains or boundary in the "Graphics" window, the domain/boundary that I click does not correspond to the domain/boundary that is selected. Sometimes, when clicking on the same boundary multiple times, sometimes different selections are made. The clicked-domain to selected-domain relationship appears to be mostly random. The only reliable way for me to choose boundaries/domains is to use the box-selector tool. I have found that this bug occurs in 4.3 and not 4.2. Versions 4.2, 4.3 and 4.3a are all installed in the station I work at. Most of my recent work was done in 4.3a and cannot be opened in 4.2.
Ivar,
I just ran the new model (halved and mirrored) with glycerol as the material using the boundary layer meshing (4 layers instead of the default 8 -- I wanted to save time, but I will increase that number if I need to). I am currently solving the same model with the power-law fluid. The flow through all of the outputs is very even (see the pictures-- glycerol at 92 psi). Fortunately, this agrees very well with experiment. I can show you videos if you don't believe me ;)
Nagi
Using a Carreau model instead of a power-law fluid is an excellent suggestion, thanks! I was using a power-law model because our fluid was experimentally measured and fit a power-law very well. It makes sense that the viscosity at a no-slip boundaries would cause problems. I will try to work out the coefficients for the Carreau model which apply over the expected shear rates.
I took your advice and rebuilt the network with the mirror symmetry (see attached). I added a boundary layer as Ivar suggested. Thanks for sending me the example swept mesh. Is the partitioning necessary because the model does not have a constant depth, and beginning and end faces of the mesh sweep must match? When meshing the new model I opted to stick with the free tetrahedral mesh because partitioning and sweeping was too tedious. The tedium was due to an...
Unrelated problem: When I am selecting domains or boundary in the "Graphics" window, the domain/boundary that I click does not correspond to the domain/boundary that is selected. Sometimes, when clicking on the same boundary multiple times, sometimes different selections are made. The clicked-domain to selected-domain relationship appears to be mostly random. The only reliable way for me to choose boundaries/domains is to use the box-selector tool. I have found that this bug occurs in 4.3 and not 4.2. Versions 4.2, 4.3 and 4.3a are all installed in the station I work at. Most of my recent work was done in 4.3a and cannot be opened in 4.2.
Ivar,
I just ran the new model (halved and mirrored) with glycerol as the material using the boundary layer meshing (4 layers instead of the default 8 -- I wanted to save time, but I will increase that number if I need to). I am currently solving the same model with the power-law fluid. The flow through all of the outputs is very even (see the pictures-- glycerol at 92 psi). Fortunately, this agrees very well with experiment. I can show you videos if you don't believe me ;)
Nagi
Using a Carreau model instead of a power-law fluid is an excellent suggestion, thanks! I was using a power-law model because our fluid was experimentally measured and fit a power-law very well. It makes sense that the viscosity at a no-slip boundaries would cause problems. I will try to work out the coefficients for the Carreau model which apply over the expected shear rates.
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
good to hear you are getting there ;)
For the regular distribution of your flow, it's only my gu's feeling that says me do not expect a symmetric geoemtry to give fully symmetric outlet flows for any case of fluid velocity, vicosity ... but I might be wrong, I tend to believe what COMSOL gives out, provided my physics is set up correctly ;)
--
Good luck
Ivar
good to hear you are getting there ;)
For the regular distribution of your flow, it's only my gu's feeling that says me do not expect a symmetric geoemtry to give fully symmetric outlet flows for any case of fluid velocity, vicosity ... but I might be wrong, I tend to believe what COMSOL gives out, provided my physics is set up correctly ;)
--
Good luck
Ivar
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Steve-
Yes, the partioning would be needed for swept meshing because the source and destination face(s) need to be compatible.
One thing I noticed that may or may not help. I saw that you have generated boundary layers on some artificial unphysical boundaries (see attached png). With these boundaries removed from the boundar layer operation, you will need to check "trimming" instead of "splitting" to avoid errors. Trimming unfortunately doesn't give you as much refinement at the corners. Overall, with this technique you will save about 6% in total element number. See attached mph. The reason I say this may not help is that you are losing some of the refinement at the corners with this technique. The benefit is that you are removing a lot of unnecessary refinement on the unphysical boundaries.
I suppose In the end, if you are getting convergence and good accuracy in a reasonable amount of time then no need to fiddle with perfecting the mesh.
Best regards,
Josh Thomas
AltaSim Technologies
Yes, the partioning would be needed for swept meshing because the source and destination face(s) need to be compatible.
One thing I noticed that may or may not help. I saw that you have generated boundary layers on some artificial unphysical boundaries (see attached png). With these boundaries removed from the boundar layer operation, you will need to check "trimming" instead of "splitting" to avoid errors. Trimming unfortunately doesn't give you as much refinement at the corners. Overall, with this technique you will save about 6% in total element number. See attached mph. The reason I say this may not help is that you are losing some of the refinement at the corners with this technique. The benefit is that you are removing a lot of unnecessary refinement on the unphysical boundaries.
I suppose In the end, if you are getting convergence and good accuracy in a reasonable amount of time then no need to fiddle with perfecting the mesh.
Best regards,
Josh Thomas
AltaSim Technologies
Attachments:
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
I ran a parameter sweep of the power-law fluid (wax) at pressures of 70, 105, 140, 175, 210, 280, 350 and 420 psi. The solutions converged for the first three pressures (70, 105, 140), But did not converge at 175. What surprised was that comsol then solved for pressures of 157.5, 148.75, 144.375, 142.19, and 141.75. At 141.75 (attached) it converged. I guess COMSOL was trying to find the max pressure it could converge at - pretty clever.
I think that at higher pressures, the shear rate difference is large and creates a large velocity gradient that is not resolvable with my current mesh. I will try re-solving the simulations with the boundary layer number set to the default of 8 (as opposed to 4). Regardless of whether this works or not, I will next try Nagi's suggestion of the Carreau model fluid with the boundary layer of 4.
Josh,
I see what you're saying about meshing the unphysical boundary-it's a good point. However, as I mentioned earlier, it is not easy for me to select large numbers of boundaries because of this apparent bug in v. 4.3 (perhaps I should make a separate thread for it). I agree that fixing the boundaries would result in a more efficient mesh, but currently I'm not too concerned with the time to solve (I set it to solve and come back in several hours).
Ivar,
This bifurcating network design yields surprisingly even outputs. Partly because the Reynolds number is so low (~10^-4 to 10^-6), it is closer to creeping flow than laminar. I could probably solve it as creeping flow (neglecting the Stoke's term), but in a few experiments it didn't seem to save much time, and I don't see the harm in solving for laminar flow. What can result in uneven flow is if the vertical section of each generation of bifurcations are very short compared to their widths. In my example, the aspect ratio is 3, but even at 1.5 it still works well. Also, I found (I think) that a network with more generations is less susceptible to unevenness of flow due to inertial effects of a short vertical section immediately after a bend.
I think that at higher pressures, the shear rate difference is large and creates a large velocity gradient that is not resolvable with my current mesh. I will try re-solving the simulations with the boundary layer number set to the default of 8 (as opposed to 4). Regardless of whether this works or not, I will next try Nagi's suggestion of the Carreau model fluid with the boundary layer of 4.
Josh,
I see what you're saying about meshing the unphysical boundary-it's a good point. However, as I mentioned earlier, it is not easy for me to select large numbers of boundaries because of this apparent bug in v. 4.3 (perhaps I should make a separate thread for it). I agree that fixing the boundaries would result in a more efficient mesh, but currently I'm not too concerned with the time to solve (I set it to solve and come back in several hours).
Ivar,
This bifurcating network design yields surprisingly even outputs. Partly because the Reynolds number is so low (~10^-4 to 10^-6), it is closer to creeping flow than laminar. I could probably solve it as creeping flow (neglecting the Stoke's term), but in a few experiments it didn't seem to save much time, and I don't see the harm in solving for laminar flow. What can result in uneven flow is if the vertical section of each generation of bifurcations are very short compared to their widths. In my example, the aspect ratio is 3, but even at 1.5 it still works well. Also, I found (I think) that a network with more generations is less susceptible to unevenness of flow due to inertial effects of a short vertical section immediately after a bend.
Attachments:
Please login with a confirmed email address before reporting spam
Posted:
1 decade ago
Hi
for me the flow path suggest that one rounds off several corner therein, to avoid stagnation regions where you might get clogging.
Indeed a t low "creep" flow rates the pressure driven flow should be rather uniform (from my guts feeling again to be checked as you have done ;),
What is the speed of a pressure wave through such a system ? , probably much shorter than the flow itself
--
Good luck
Ivar
for me the flow path suggest that one rounds off several corner therein, to avoid stagnation regions where you might get clogging.
Indeed a t low "creep" flow rates the pressure driven flow should be rather uniform (from my guts feeling again to be checked as you have done ;),
What is the speed of a pressure wave through such a system ? , probably much shorter than the flow itself
--
Good luck
Ivar