Add damping and stiffness coefficients for fluid flow

[JuliaMota]

Currently in fluid flow, it is possible to find the pressure field in the journal bearings. From the pressure field, the stiffness and damping coefficients can be calculated. For this, it is first necessary to calculate the oil film forces in the directions N and T, i.e., the direction opposite the eccentricity and tangential direction.

[flaviorangel]

We are still working on this issue. Even though the pressure matrix is being calculated correctly based on two references we are using, and even though the error in both N and T seems to be small enough, we haven't been able to obtain the stiffness and damping matrices correctly. One of our concerns was if the attitude angle, eccentricity, and load are correct, so we tested how the first two behave as we increase the spin speed.

Using the fluid_flow_short_load example from test_fluid_flow.py, starting with speed 0, then 1, then increasing exponentially, we plotted the first eccentricity and then only the center of the rotor. The sommerfeld seems to be behaving correctly. Here is the result:

[flaviorangel]

We found out that the program is not being able to calculate a correct eccentricity based on the load. If the eccentricity is given, the returned load makes sense with the rest of the code, but if the load is given instead, the eccentricity error, even though small is enough to cause significant error throughout the rest of the code. Here is an example of when the eccentricity is given:

nz = 8 ntheta = 132 nradius = 11 omega = 100.2np.pi/60 p_in = 0 p_out = 0 radius_rotor = 0.1999996 radius_stator = 0.1999996 + 0.000194564 length = (1/10) (2 radius_stator) eccentricity = 0.0001 viscosity = 0.015 density = 860

The returned load is 42.88929075582885.

The resulting forces are:

força x: -0.7383085341492546 ;force y: 43.57995114226335 força x ishida: 0.0 ;force y ishida: 42.84760766277614 força x java: 3.552713678800501e-15 ;force y java: 42.84760766277615 força x friswell: 7.105427357601002e-15 ;force y friswell: 42.889290755828846

The resulting force along the x-axis should be zero and the one along the y-axis should be equal to the load. Notice that the results are reasonably close.

Now, if we instantiate a FluidFlow with the same parameters as above, but giving the load instead, this is the eccentricity the program returns: 0.0001470009050421523. Even though the error seems small for a small number, it is enough to mess with the resulting forces:

força x: -2.51858175322009 ;force y: 204.66216491157854 força x ishida: 1.4210854715202004e-14 ;force y ishida: 199.746819623002 força x java: 0.0 ;force y java: 199.74681962300207 força x friswell: 1.4210854715202004e-14 ;force y friswell: 199.9411377127039

The resulting force along x-axis makes sense, as it is reasonably close to zero (ours still require calibration), but the resulting force along y-axis does not make sense.

As we could not find an error in our code, we concluded that the Friswell attempt to calculate the eccentricity based on the load is a bold one, not attempt by his pairs, as it involves a quartic equation in eccentricity_ratio**2. Therefore, the user ideally should not give the load and let the program calculate it.

[JuliaMota]

With the current code of PR #357 we can numerically calculate the stiffness matrix. However, we have not yet found a way to correctly calculate the damping matrix. I am studying more about it and soon we will generate more results.

[flaviorangel]

We are still stuck on how to calculate the damping matrix. The authors suggest a small perturbation in the journal center velocities to do so. As the system works as a spring–mass one, a change in the force is expected. This way, we should be able to calculate the damping matrix, given that: F = -kx -cv F¹ = -kx -cv¹ c = (F - F¹)/(v¹ - v)

The problem is we are not sure how to implement this perturbation. Because our rotor works in steady state, these velocities are zero (it's not moving) and we cannot add them given our modeling. We have tried a few attempts on how to calculate the damping matrix and have been comparing with Friswell's non numerical way. The following results were obtained using the FluidFlow object from the fluid_flow_short_eccentricity function in test section.

These are the values expected for the C matrix based on Friswell: cxx = 61544.53796395076 cxy = -46951.22012390657 cyx = -46951.22012390657 cyy = 138845.4565778727

Attempt 1- Changing angular frequency

A small change in angular frequency moves the equilibrium position of the rotor (if every related attributes are recalculated). The difference between previous and later x and y should stand for velocity in above equation. Omega was changed in plus 10%:

cxx = -99051.9665987171, rel_error = -2.609435538483302 cxy = -8952.526726291537, rel_error = 0.8093228098723445 cyx = -26686194.25908057, rel_error = 567.3812729180285 cyy = -2411954.810501671, rel_error = -18.371506925391575

Omega plus 1%:

cxx = -184410.69544226676, rel_error = -3.99637793284408 cxy = -10492.997309889843, rel_error = 0.7765127874803188 cyx = -43988236.10346737, rel_error = 935.8922892180492 cyy = -2502937.489571938, rel_error = -19.02678712909949

[flaviorangel]

Attempt 2- Changing angular frequency, but replacing the rotor velocity with the fluid velocity

The challenge of this approach is calculating a radial velocity for the rotor using the one of the fluid. If we take the fluid that is touching the wall of the rotor, this value is equal to the angular frequency. We can then decompose this frequency in its x and y components, but these values will vary with the angle. A solution would be the sum of the x and y components of every point along theta. This result is zero for both axis if we take the fluid touching the wall of the rotor, but a value different than zero should occur for another position along the gap. Omega plus 10%:

cxx = -0.6291838124530058, rel_error = -1.0000102232274912 cxy = -0.35401891044994216, rel_error = 0.9999924598570704 cyx = -169.5122471602426, rel_error = 0.9963896093283009 cyy = -95.37839318151977, rel_error = -1.0006869392454913

Omega plus 1%: cxx = -0.5360541802297808, rel_error = -1.0000087100203847 cxy = -0.5853065906849947, rel_error = 0.9999875337299192 cyx = -127.86719223441372, rel_error = 0.9972765949021779 cyy = -139.61557079007676, rel_error = -1.0010055465568062

(Notice that the erro seems smaller just because the actual values are way smaller than they were supposed to be; that way, the relative error fraction is approaching one)

[flaviorangel]

Attempt 3- Changing angular frequency, keeping the rotor in same position, replacing the rotor velocity with the fluid velocity

This time we don't recalculate the rotor position when changing the angular frequency. Even though the rotor will certainly be out its equilibrium position, we would hope the C matrix will have a greater impact than the K matrix in the force change.

Omega plus 1%: cxx = 0.004670117762131378, rel_error = 0.999999924118079 cxy = 0.00612167775716015, rel_error = -1.0000001303837842 cyx = -10.58282618085178, rel_error = 0.9997745995492147 cyy = -13.872166600279822, rel_error = -1.0000999108429054

Omega plus 10%: cxx = 0.004670117762216926, rel_error = 0.999999924118079 cxy = 0.006121677757348445, rel_error = -1.0000001303837842 cyx = -10.58282618082552, rel_error = 0.9997745995492151 cyy = -13.872166600417975, rel_error = -1.0000999108429065

[ross-bott]

Hi there! I have marked this issue as stale because it has not had activity for 45 days. Consider the following options:

• If the issue refers to a large task, break it in smaller issues that can be solved in less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.

[JuliaMota]

As we were not successful in the attempts described above, I think that to find the damping matrices it is necessary to consider the time variable in the equations. We currently have only the stationary case. However, introducing this variable can be very laborious and cause significant changes in the model.

On the other hand, I have the suggestion that the movement of the rotor in relation to time may be directly linked to Sommerfeld's number. This has already been confirmed in the case of approaches to the short bearing. We would need to find a way to generalize this concept to all bearing sizes, considering the numerical results.

[JuliaMota]

Since some problems of this issue have already been resolved, it will be closed. The discussion on the damping coefficients will be recorded in the issue #435 .