[Floing]

Hello,

i have a question concerning my simulation of a bubble rising in water with an free surface model like in the tutorial " Free surface flow over a bump" (which works fine so far).

The geometry i am using is an 4 degree piece of a cylinder with my bubble in the center of it. What i'm tying to do is to implement an inlet boundary condition above my risng bubble with its rising velocity as inlet velocity. With this inlet condition the bubble should (almost) remain at its place at the center of my geometry and should not rise.

This velocity should vary over time since the bubble will rise faster with time. I am trying to find an volume weighted average velocity for the whole bubble and already used the following expression to calculate this velocity but it did not work as expected:

volumeAve(Air.Velocity v)@Blase

This Expression does not change over time/time step and calculates very small velocitis (which could be because the average velocity is made out of cells where no air is present which leads to an smaller average velocity)

Is it possible to set something like this up? I would be very thankfull if you can help me.

Floing

[ghorrocks]

centre of mass of the bubble could be something like:

volumeInt(y*Air.vf)@Blase/volumeInt(Air.vf)@Blase

This should give the y height of the centre of mass of the bubble at that time step.

[Floing]

Hello Glenn,

thanks for your help. I tried to model my velocity with your function but it resulted in the following error:

The function 'volumeInt' referenced in parameter 'Normal Speed' in object '/FLOW:Flow Analysis 1/DOMAIN:Blase/BOUNDARY:Oeffnung/BOUNDARY CONDITIONS/MASS AND MOMENTUM' has an invalid argument, 'y*Air.Volume Fraction'. Only arguments that consist of a single recognised variable name are supported by the solver.

Floing

[ghorrocks]

To work around that you can define an additional variable with the expression y*Air.vf. Then the function becomes volumeInt(AV)@Blase/volumeInt(Air.vf)@Blase

[Floing]

Hello again,

thank you Glenn, I found this expression working for me as long as I dont use it as inlet velocity:

volumeInt(Air.Velocity v)@Blase/volumeInt(Air.Volume Fraction)@Blase

with this definition from Ansys help:

"23.4.1.4. Volume Integral

A volume integral is computed by summing the product of the cell volume and the selected field variable"


This should model my volume average velocity for my air bubble but everytime I set this as my inlet velocity, the solver crashes with an fatal overflow error. I know that ther is an FAQ to this but I dont know if it can solve my problem (Already tried a finer mesh with the same result and smaller, adaptiv time steps). I assume that my physics are correct since my simulation works without this velocity and double precision is always on.

Does somebody have an idea why this happens?

Floing

[ghorrocks]

The problem with your equation is it will average the entire velocity field, not just the velocity field in the bubble. You probably want to multiply Air.Velocity V by Air.vf so the velocity field in the water is ignored.

If you are setting your inlet velocity to this expression I suspect you will have problems converging. The easiest way to do this will be to set the inlet velocity to a guessed value and see the result. Then adjust the inlet velocity and rerun with a closer value. Continue until you are close enough.

[Floing]

Do you know why this does not converge? Ther is no obvious reason for me why it should not work.

[ghorrocks]

For it to converge it needs to be numerically stable. Google numerical stability and you will find a zillion reasons why it might not converge.

[Floing]

Hi,

thank you Glenn. I tried to make my simulation numerical stable but I failed to do so. The problem is, that I have to implement some sort of controller for my inlet velocity, so I can’t just guess an inlet velocity like you mentioned. Furthermore it doesn’t converge with a fixed value either. Different Mesh sizes were used to solve my Bubble radius with 10, 15, 20 and 25 cells per radius. The 4 degree extend is done with one cell via sweep method.

The added picture shows my results for 20 and 25 Mesh cells. The velocity should rise to about 30 cm/s for a 2mm bubble but it stops at around 16/9 cm/s. Three different Methods were used to initialize my bubble (step-Function/tanh/CEL-Function) but the effect for my overall velocity is small. Currently I’m using the tanh function which I found in the Ansys help section to initialize a smeared bubble or drop (Picture). I’m adding 10% of the Auftrieb 1 expression to the current velocity at the upper opening to hold my bubble. Bigger values than this lead to an even faster instability.
I know that it isn’t somebodies job to solve my problems but a suggestion would mean a lot to me.


My Geometry is still the same like in my first post and this is my setup:

LIBRARY:
CEL:
EXPRESSIONS:
Auftrieb 1 = volumeInt(Air.Produkt)@Blase/volumeInt(Air.Volume \
Fraction)@Blase
Auftrieb 2 = volumeInt(Air.Volume Fraction*Air.velocity \
v)@Blase/volumeInt(Air.Volume Fraction)@Blase
Blasenradius = 1 [mm]
Geschwindigkeit = Air.Volume Fraction*Air.Velocity v
VF Int = Verschmierung(1[mm]-sqrt(x^2+(y-6[mm])^2+z^2))
VFA = step((Blasenradius-radius)/1 [m])
VFAir = 0.5*tanh(dist/delta)+0.5
delta = 0.05[mm]
dist = Blasenradius-radius
radius = sqrt(x^2+(y-6[mm])^2+z^2)
uinlet = abs(areaAve(Water.Velocity)@Oeffnung)
uinletnew 1 = uinlet+0.1*Auftrieb 1
uinletnew 2 = uinlet+0.1*Auftrieb 2
END
FUNCTION: Verschmierung
Argument Units = mm
Option = Interpolation
Result Units = m/m
INTERPOLATION DATA:
Data Pairs = -0.1,0,0.1,1
Extend Max = On
Extend Min = On
Option = One Dimensional
END
END
END
ADDITIONAL VARIABLE: Produkt
Option = Definition
Tensor Type = SCALAR
Units = [m s^-1 ]
Variable Type = Unspecified
END

FLOW: Flow Analysis 1
SOLUTION UNITS:
Angle Units = [rad]
Length Units = [m]
Mass Units = [kg]
Solid Angle Units = [sr]
Temperature Units = [K]
Time Units = [s]
END
ANALYSIS TYPE:
Option = Transient
EXTERNAL SOLVER COUPLING:
Option = None
END
INITIAL TIME:
Option = Automatic with Value
Time = 0 [s]
END
TIME DURATION:
Option = Total Time
Total Time = 0.1 [s]
END
TIME STEPS:
First Update Time = 0.0 [s]
Initial Timestep = 1e-6 [s]
Option = Adaptive
Timestep Update Frequency = 1
TIMESTEP ADAPTION:
Maximum Timestep = 1e-4 [s]
Minimum Timestep = 1e-6 [s]
Option = Number of Coefficient Loops
Target Maximum Coefficient Loops = 5
Target Minimum Coefficient Loops = 2
Timestep Decrease Factor = 0.8
Timestep Increase Factor = 1.2
END
END
END
DOMAIN: Blase
Coord Frame = Coord 0
Domain Type = Fluid
Location = Blase
BOUNDARY: Boden
Boundary Type = OPENING
Location = Boden
BOUNDARY CONDITIONS:
FLOW DIRECTION:
Option = Normal to Boundary Condition
END
FLOW REGIME:
Option = Subsonic
END
MASS AND MOMENTUM:
Option = Opening Pressure and Direction
Relative Pressure = 0 [Pa]
END
END
FLUID: Air
BOUNDARY CONDITIONS:
VOLUME FRACTION:
Option = Value
Volume Fraction = 0
END
END
END
FLUID: Water
BOUNDARY CONDITIONS:
VOLUME FRACTION:
Option = Value
Volume Fraction = 1
END
END
END
END
BOUNDARY: Mantel
Boundary Type = WALL
Location = Mantel
BOUNDARY CONDITIONS:
MASS AND MOMENTUM:
Option = Free Slip Wall
END
END
FLUID PAIR: Air | Water
BOUNDARY CONDITIONS:
WALL ADHESION:
Option = None
END
END
END
END
BOUNDARY: Oeffnung
Boundary Type = OPENING
Location = Oeffnung
BOUNDARY CONDITIONS:
FLOW REGIME:
Option = Subsonic
END
MASS AND MOMENTUM:
Option = Cartesian Velocity Components
U = 0 [m s^-1]
V = -uinletnew 1
W = 0 [m s^-1]
END
END
FLUID: Air
BOUNDARY CONDITIONS:
VOLUME FRACTION:
Option = Value
Volume Fraction = 0
END
END
END
FLUID: Water
BOUNDARY CONDITIONS:
VOLUME FRACTION:
Option = Value
Volume Fraction = 1
END
END
END
END
BOUNDARY: W1 1
Boundary Type = SYMMETRY
Location = W1,W2
END
DOMAIN MODELS:
BUOYANCY MODEL:
Buoyancy Reference Density = 997 [kg m^-3]
Gravity X Component = 0 [m s^-2]
Gravity Y Component = -g
Gravity Z Component = 0 [m s^-2]
Option = Buoyant
BUOYANCY REFERENCE LOCATION:
Option = Automatic
END
END
DOMAIN MOTION:
Option = Stationary
END
MESH DEFORMATION:
Option = None
END
REFERENCE PRESSURE:
Reference Pressure = 1 [atm]
END
END
FLUID DEFINITION: Air
Material = Air at 25 C
Option = Material Library
MORPHOLOGY:
Option = Continuous Fluid
END
END
FLUID DEFINITION: Water
Material = Water
Option = Material Library
MORPHOLOGY:
Option = Continuous Fluid
END
END
FLUID MODELS:
ADDITIONAL VARIABLE: Produkt
Option = Fluid Dependent
END
COMBUSTION MODEL:
Option = None
END
FLUID: Air
ADDITIONAL VARIABLE: Produkt
Additional Variable Value = Geschwindigkeit
Option = Algebraic Equation
END
FLUID BUOYANCY MODEL:
Option = Density Difference
END
END
FLUID: Water
ADDITIONAL VARIABLE: Produkt
Additional Variable Value = Geschwindigkeit
Option = Algebraic Equation
END
FLUID BUOYANCY MODEL:
Option = Density Difference
END
END
HEAT TRANSFER MODEL:
Fluid Temperature = 25 [C]
Homogeneous Model = True
Option = Isothermal
END
THERMAL RADIATION MODEL:
Option = None
END
TURBULENCE MODEL:
Option = Laminar
END
END
FLUID PAIR: Air | Water
Surface Tension Coefficient = 0.072 [N m^-1]
INTERPHASE TRANSFER MODEL:
Option = Free Surface
END
MASS TRANSFER:
Option = None
END
SURFACE TENSION MODEL:
Curvature Under Relaxation Factor = 0.5
Option = Continuum Surface Force
Primary Fluid = Water
Volume Fraction Smoothing Type = None
END
END
INITIALISATION:
Option = Automatic
FLUID: Air
INITIAL CONDITIONS:
VOLUME FRACTION:
Option = Automatic with Value
Volume Fraction = VFAir
END
END
END
FLUID: Water
INITIAL CONDITIONS:
VOLUME FRACTION:
Option = Automatic
END
END
END
INITIAL CONDITIONS:
Velocity Type = Cartesian
CARTESIAN VELOCITY COMPONENTS:
Option = Automatic with Value
U = 0 [m s^-1]
V = 0 [m s^-1]
W = 0 [m s^-1]
END
STATIC PRESSURE:
Option = Automatic with Value
Relative Pressure = 0 [Pa]
END
END
END
MULTIPHASE MODELS:
Homogeneous Model = On
FREE SURFACE MODEL:
Interface Compression Level = 2
Option = Standard
END
END
END
OUTPUT CONTROL:
MONITOR OBJECTS:
MONITOR BALANCES:
Option = Full
END
MONITOR FORCES:
Option = Full
END
MONITOR PARTICLES:
Option = Full
END
MONITOR POINT: Auftrieb
Coord Frame = Coord 0
Expression Value = Auftrieb 1
Option = Expression
END
MONITOR POINT: Durchschnitt
Coord Frame = Coord 0
Expression Value = uinlet
Option = Expression
END
MONITOR POINT: Summe
Coord Frame = Coord 0
Expression Value = uinletnew 1
Option = Expression
END
MONITOR RESIDUALS:
Option = Full
END
MONITOR TOTALS:
Option = Full
END
END
RESULTS:
File Compression Level = Default
Option = Standard
END
TRANSIENT RESULTS: Results
File Compression Level = Default
Include Mesh = No
Option = Selected Variables
Output Variables List = Absolute Pressure,Air.Velocity,Water.Velocity \
w,Water.Velocity v,Water.Velocity u,Velocity u,Velocity w,Velocity \
v,Air.Velocity u,Air.Velocity v,Air.Velocity \
w,Water.Velocity,Water.Superficial Velocity,Air.Superficial \
Velocity,Pressure,Air.Volume Fraction,Water.Volume Fraction
OUTPUT FREQUENCY:
Option = Time Interval
Time Interval = 0.0002 [s]
END
END
END
SOLVER CONTROL:
ADVECTION SCHEME:
Option = High Resolution
END
CONVERGENCE CONTROL:
Maximum Number of Coefficient Loops = 30
Minimum Number of Coefficient Loops = 1
Timescale Control = Coefficient Loops
END
CONVERGENCE CRITERIA:
Residual Target = 0.000001
Residual Type = RMS
END
MULTIPHASE CONTROL:
Initial Volume Fraction Smoothing = Volume-Weighted
END
TRANSIENT SCHEME:
Option = Second Order Backward Euler
TIMESTEP INITIALISATION:
Option = Automatic
END
END
END
END

Floing

[ghorrocks]

Have you done a mesh sensitivity study? Also time step size? Convergence tolerance?

Also for 2D axisymmetric wedges you should check the sensitivity to the wedge angle.

Finally: How big are your bubbles? Does the skin of your bubble act as a slip or a no-slip surface to the liquid? Some surfactants can act as a no-slip surface to the liquid as they form a strong skin around the air bubble which does not move relative to the bubble. For liquids with low surfactancy or large bubbles they act more like slip conditions. You might be somewhere in the middle to make things complex

[Floing]

Hello Glenn,

thank you for your ongoing help. My bubble diameters reach from 0.5-4mm (the bubble in the postet picture above has a diameter of 2mm) for my air bubble in water. I never asked myself if this would be free or no slip, but as you mentioned I guess it would be more like free slip.

About the timesteps, I thought this could be done with adaptive timesteping and set Coefficient Loops to 2-5 (I will set the lower border to a value my timestep will never reach). For convergence tolerance I used 10e-4, 10e-5 and 10e-6 but I could not find a point were my simulation runs stable or the results were near my wanted values. In a paper with a similar problem (bubble rise without the controlle part) 10e-6 was used as tolerance. I get this small OK’s from my solver during the simulation which should mean that it doesn’t converge in the desired Area. Should I use even lower values? My problem seems to be somewhere else.

So, in my opinion, the Mesh would be the point to look for after. This are the Mesh statistics from the CFX-Solver (although I know that there are more points to look after then just those three mentioned here).

Domain Name | Orthog. Angle | Exp. Factor | Aspect Ratio |
+----------------------+---------------+--------------+--------------+
| | Minimum [deg] | Maximum | Maximum |
+----------------------+---------------+--------------+--------------+
| Blase | 11.2 ! | 2 OK | 14 OK |
+----------------------+---------------+--------------+--------------+
| | %! %ok %OK | %! %ok %OK | %! %ok %OK |
+----------------------+---------------+--------------+--------------+
| Blase | <1 0 100 | 0 0 100 | 0 0 100 |
+----------------------+---------------+--------------+--------------+
Maximum Face Angle for Case Fluiddynamik CFX

Min: 90 [degree]
Max: 92.0001 [degree]
Minimum Face Angle for Case Fluiddynamik CFX

Min: 3.99993 [degree]
Max: 88.0001 [degree]

(Added pictures of my Mesh)

The FAQ says that the Aspect ratio has to be around 1-2 for surface tension modeling. Since my bubble is in the center of my geometry the ratio there is close to this value and gets bigger with my wedge radius but I don’t want to model any surface tension at the edge of my wedge so this should be ok. I don’t get how I can increase the orthogonal Angle of my Mesh but I did not pay much attention to it since the number of cells which are below the desired value is <1% (which could be one of my countless mistakes). Like I mentioned in my previous post I used 4 different mesh cell sizes to model my bubble but nothing came close to the results I want for a bubble with a diameter of 2mm (30 cm/s rising velocity). In my understanding this would be some sort of mesh sensitive study (if I am mistaken here I’m very sorry). About the wedge angle, I was told that using a 4-degree Wedge would be common for problems with rotational periodicity. This would be a point to look after.

My overall physics are modeled after similar bubble cases in this forum and from papers.



Floing

PS: I just recogniced that the aspect ratio for the first 4 cells is not between 1-2 it goes up to 12. Can this be the answer to my problem or is there still something else I forgot? I feel a bit like an idiot , since it seems my problem is an very obious one. In my oppinion I have to improve my mesh but I dont really know in which direction (This would be a question for another forum part)

There is always the possibility that I made a major mistake somewhere since this are my first steps in CFD-simulation. Would someshing like that be easier to model as quasi 2d with an rectangle geometry extruded for one mesh cell? It should be easier to match the mesh criteria.

[ghorrocks]

I guess it would be more like free slip.

That does not sound like a very good reason to me. I think you should do some research into this and determine the answer properly.

Your second paragraph suggests the problem is indeed somewhere else. Fix the problem so it runs stably before doing mesh, convergence tolerance and time step size sensitivity studies.

A 4 degree wedge is typical for 2D axisymmetric models. However, surface tension modelling is highly sensitive to mesh quality - far more sensitive than any just about other model. So your wedge angle of 4 degrees could be a problem here.

I would suggest trying modelling a 3D rectangular block with perfectly 1:1:1 aspect ration hexes. You can put symmetry planes in two sides, but it does mean you are modelling 25% of the 3D geometry rather than a single slice. But your model is small so it should run quickly.

I think you will find the requirement for accurate surface tension modelling is much tighter than 1-2 (where did you get that from?) In my experience, if you want the Laplacian pressure in a bubble to be correct within about 10% you need to keep aspect ratio <1.2.