Rayleigh Taylor Instability Challenge

Q1. What are some practical CFD models that have been based on the mathematical analysis of Rayleigh Taylor waves? In your own words, explain how these mathematical models have been adapted for CFD calculations.

Rayleigh waves are a type of surface acoustic waves that travel along the surface of solids. They can be produced in materials in many ways, such as by a localized impact or by piezo-electric transduction, and are frequently used in non-destructive testing for detecting defects. Rayleigh waves are part of the seismic waves that are produced on the earth by earthquakes. When guided in layers they are referred to as Lamb waves, Rayleigh–Lamb waves, or generalized Rayleigh waves.

Q2. Perform the Rayleigh Taylor instability simulation for 3 different mesh sizes with the base mesh being 0.5 mm. Compare the results by showing the animations( Attach them from your google drive). Also, explain why a steady-state approach might not be suitable for these types of simulation.

In this challenge, Rayleigh Taylor Instability has been studied. For this setup geometry of 20cm*20cm block is used one above other. Common layer in-between is shared using Share Prep option in Ansys Workbench module.

In the Meshing module, three mesh sizes are used and a difference in the animations are shown by the end of the report. Mesh sizes are 0.5mm 0.25mm and 0.15mm

Simulation Setup: Transient solver with Pressure based type simulation is used. Reason for not using Steady-state simulation is discussed further. Since there is not much concern about temperature energy equation is not used. The flow model is taken as Laminar flow. This involves two phases (water and air) in Multiphase - Volume fraction is turned on which takes care of two phases interaction. In the material library, Water is loaded from the fluent database and the default value is used for simulation. Under Solving module Patch is done in which the first Volume fraction option is used and for water case, fff-water-surface is provided as 1 and fff-air-surface is provided as 0. Before this step water and air are assigned at proper places. As a sanity check, velocity contours are made and verified as top part shows unit velocity and the bottom part shows zero velocity.

Case 1:

Mesh Size 0.5mm

Number of Elements: 3200

Number of Nodes: 3321

Animation:

Case 2:

Mesh Size 0.25mm:

Number of Elements: 12800

Number of Nodes: 13041

Animation:

Case 3:

Mesh Size 0.15mm:

Number of Elements: 35378

Number of Nodes: 35778

Animation:

RESULTS ANALYSIS:

For comparing results five points were taken in the domain and their air volume fraction, water volume fraction and velocity profile are plotted respectively to compare the variation of result with each mesh size. Five different lines used for the comparison are shown below:

Line (0,0.1) and (0.02,0.01) [Along x-axis]:

Line (0,0.02) and (0.02,0.02) [Along x-axis]:

Line (0,0.03) and (0.02,0.03) [Along x-axis]:

Line (0.005,0) and (0.005,0.04) [Along y-axis]:

Line (0.015,0) and (0.015,0.04) [Along y-axis]:

Thus the plots show that velocity seems to be stabilized somewhat for fine refined mesh than other two mesh sizes. One issue is only 100 sampling points were used for plotting these data, the graph would have been better if there are still a high number of sampling points used. As well the air volume fraction and water volume fraction sums to unity which proves the accuracy of physics of the problem. Along the x-axis, there are a lot of fluctuations throughout three mesh sizes so no evident result can arrive for this case. But on the other hand along the y-axis, we can see the more stable curve in more refined mesh than other tow mesh sizes used.

Throughout the simulation, report files have been generated and some Data parsing has been done using Matlab:

clc
clear all
close all

A=load('report-file-0.txt');
flow_time=A(:,2);
volume_integral_velocity=A(:,5);
air_surface_velocity=A(:,4);
water_surface_velocity=A(:,6);

subplot(2,1,1)
plot(flow_time,volume_integral_velocity);
xlabel('Flow time(s)')
ylabel('Volume Integral velocity (m/s)')
title('Volume Integral Velocity Vs Flow time for Mesh size 0.25mm')
grid on
% figure(2)
% plot(flow_time,air_surface_velocity);
% figure(3)
% plot(flow_time,water_surface_velocity);

B=load('case 2 output files.txt');
time=B(:,1);
volume_average_velocity=B(:,2);
flow_time=B(:,3);
volume_integral_velocity=B(:,5);
% figure(4)
% plot(time,volume_integral_velocity)
% xlabel('Time steps count')
% ylabel('Volume Integral velocity (m/s)')
subplot(2,1,2)
plot(flow_time,volume_integral_velocity)
xlabel('Flow time (s)')
ylabel('Volume Integral velocity (m/s)')
title('Volume Integral Velocity Vs Flow time for Mesh size 0.15mm')
grid on

In this case, two mesh size of 0.25mm and 0.15mm report files of flow time, delta time, volume integral velocity, volume integral average velocity were extracted for every iteration. Then Velocity volume integral vs flow time were plotted for two cases as shown below:

It can be seen that the velocity profile has been gradually starting to converge at a lower time in 0.15mm mesh size than others. In other words, Rayleigh instability are getting converged quickly in much fine mesh size. Even mesh size of 0.15mm yields just 33K nodes there is a lot of room for improvements in mesh size. While using mesh size of 0.1mm an error occurred as:

which gives a major shoot-up in the continuity equation. This might be due to CFL factor which goes beyond a certain number and provides an unstable solution since CFL is inversely proportional to mesh size. In that case, CFL number should be given much lower value and simulation should be made (in all three cases CFL number is 0.2)

The issue using Steady-state solver:

The first case is since it is using a Coupled form of Pressure velocity relation convergence occurs by the end of the second iteration even residual is 1e-6. Even at the worst case if the mesh is made super fine, the steady-state solver will do all calculations based on iterations level, but we need to see the result as function as time evolution. As well no-slip boundary condition will be activated in this steady-state solver which deviates physics of the problem. Converged solution for steady-state solver:

Q3. Define the Atwood Number. Find out the Atwood number for this case and explain how the variation in Atwood number affects the behavior of the instability. How can we use the Atwood number to validate our simulation result?

The Atwood number is a dimensionless number in fluid dynamics used in the study of hydrodynamic instabilities in density stratified flows. It is a dimensionless density ratio related to the density of heavier fluid and the density of a lighter fluid. (The difference in the fluid densities divided by their sum). Atwood number is an important parameter in the study of Rayleigh–Taylor instability. For Atwood number close to 0, RT instability flows take the form of symmetric “fingers” of fluid; for Atwood number close to 1, the much lighter fluid “below” the heavier fluid takes the form of larger bubble-like plumes.

Atwood Number = `(rho1-rho2)/(rho1+rho2)` = (997-1.225)/(997+1.225) = 0.9975.

As defined above Atwood number, in this case, is close to unity, which provides light fluid below and heavy fluids form bubble-like plumes