After launching FENSAP-ICE, create a new project directory by clicking on the icon:

Enter the name of the new project directory in the Project name box, and browse to position it within your home directory.

A message window will ask about the unit settings. Accept the defaults to keep SI units for this project.

Create a new run with FENSAP as the flow solver, by clicking on the icon:

Name this run: FENSAP_clean_4deg. A new run corresponds to a grid file icon, followed by the configuration and solution icons as shown below:

Download the 2_In-Flight_Icing.zip file here .

Unzip 2_In-Flight_Icing.zip to your working directory.

Assign the grid file naca0012.grid provided in the tutorials subdirectory ../workshop_input_files/Input_Grid/Naca0012/ by double-clicking on the grid icon and by browsing to this directory.

Double-click the config icon (gear) to open the graphical and input parameter windows of FENSAP.

The graphical window displays the airfoil and boundary edges in the far field. To zoom in on the airfoil, use the mouse wheel or Ctrl + left-click to draw a bounding box around it. To center the airfoil, press Backspace or right-click the axis marker and select Fit to view. The option Reset view returns to the original view.

Solve the Navier-Stokes equations (viscous flow) with the energy equation (Full PDE). For turbulence, select the Spalart-Allmaras model with low free stream turbulence (use a very low Eddy/laminar viscosity ratio of 1.e-5). Turbulence is then only generated by the airfoil. Surface roughness will not be considered in this calculation as this is a clean airfoil.

Set the following reference flow values in the Conditions panel:

The Reynolds number, Mach number and Adiabatic stagnation temperature are automatically updated by the graphical interface.

The initial solution (pressure, temperature, density) is set from the reference flow conditions of step 8. Impose the initial velocity using Velocity angles and an Angle of attack (AoA) of 4 degrees and hit Enter. You can check the initial velocity vector displayed in the graphical window. Click the cube icon   to toggle the display.

Go to the Boundaries panel to set the boundary conditions.

FENSAP-ICE uses a 4-digit boundary condition designation scheme where the first digit determines the main category: 1000 boundary conditions are Inlets, 2000 boundary conditions are Walls, 3000 boundary conditions are Exits, and 4000 boundary conditions are Symmetry planes. Specifically, 4100 is for X-symmetry, 4200 is for Y-symmetry, and 4300 is for Z-symmetry. See The Boundary Face Table within the FENSAP-ICE User Manual for more details.

Click BC_1000. Set the Type as Supersonic or far-field, for which all primitive variables are required, and click Import reference conditions to automatically set the boundary conditions from the reference and initial values set in steps 8 and 9.

Click BC_2001. Surface type shall be No-slip, with a specified temperature on the wall. Specifying a surface temperature produces heat fluxes from the airfoil surface to the air which will be used by ICE3D to calculate heat transfer from the water and ice surface. The final surface temperature is calculated by ICE3D, and the temperature set at this step is discarded. The value of the surface temperature should be several degrees above the adiabatic stagnation temperature in order to compute heat fluxes with the correct sign on the entire aircraft surface. For convenience, right-click in the Temperature box, choose Copy from… → Adiabatic stagnation temperature + 10 to assign the surface temperature. The BC_2001 settings should look like the following:

Repeat this step for the boundary conditions 2002, 2003, and 2004.

The BC_4300 boundary condition is the Z-symmetry plane that does not require any boundary condition specifications.

Go to the Solver panel. Select the Steady option in the Time integration box. Set the value of the CFL number to 50. With this strategy, the local time step is computed from the characteristic velocity, speed of sound and length of each element.

Set the Maximum number of time steps to 300 to achieve steady-state.

Uncheck the Use variable relaxation option. This option is sometimes helpful to improve robustness when the solution diverges after a few iterations due to numerical instabilities.

Select the Streamline upwind (SU) artificial viscosity scheme with Cross-wind dissipation set at 1e-9 and the order at 100% Second order. For more details on artificial dissipation, see Artificial Dissipation in the FENSAP-ICE User Manual.

Go to the Out panel and save the flow solution every 20 iterations by overwriting the solution file.

Compute the forces acting on the airfoil by selecting the option Drag direction based on inlet BC under Forces. In this case, the drag direction matches the angle of attack imposed on the inlet BC_1000. Select the positive Y as the lift direction.

Set the Reference area of the airfoil to 0.05334 m² to compute the lift and drag coefficients. This is the planform area of the airfoil as it appears in the grid. For correct lift and drag coefficient calculations, the planform area should be accurately specified. The lift and drag forces will be updated on Graphs at every 20 iterations as set in step 12.

Click the Run button to open the execution menu and set the Number of CPUs to 4 or more, if possible.

Launch the calculation by clicking on the Start menu button.

Go to the Graphs panel to monitor the convergence of the run.

The first graph is the average residual which shows the convergence of the flow equations. This value should reduce by at least three orders of magnitude from its initial value to obtain a suitable solution. Ideally it should reach 1e-15 (machine zero) for perfect convergence. Many things affect the level of convergence, primarily the quality and resolution of the mesh. For example, a coarse boundary layer mesh will not be able to produce perfect convergence since there is not enough resolution to properly capture the velocity profiles. In 3D unstructured grids, transitioning too soon to isotropic tetras above a fine boundary layer can also have the same effect.

The lift and drag coefficients, as global integral quantities, are also shown in Graphs panel. At convergence, when the simulation reached steady state, these values should not change. Mass inflow, outflow, and deficit graphs show the evolution of the total mass flow rates through inlets (negative) and exits (positive), and their difference which should be small at convergence (mass conservation). These quantities are computed every time the solution is saved, as specified in the Out panel. If flow calculations diverge (residuals keep going up and reaching very high values), it is the result of improper specification of inlet and exit boundary conditions in most cases and the mass flow graphs will show amplifying oscillations.

GMRES - Navier-Stokes – Residual change graph shows the performance of linear system convergence. The residual change should ideally be lower than 0.1. Small values (<0.5) mean that the linear system is solved successfully; while a value of 1.0 means that the solution process is not updating the results any longer. This can happen if the CFL is too high, if the mesh quality is bad, or if the boundary conditions are imposed incorrectly.

At convergence, the lift and drag coefficients should read 0.46259 and 0.00958 respectively.

Figure 2.2: Average Residual and Lift Coefficient Convergence

Click the View with… (or View) button at the bottom of the Graphs panel to view the solution. If a default post processor is not set, the application prompts you to choose one (Viewmerical, Fieldview, etc.). Choose Viewmerical to view the solution in that case. It is recommended to set Viewmerical as the default post processor (See ICE3D Ice Accretion on the NACA0012) since no file format conversion is needed.

In Viewmerical, you can choose which data field to display in the Data panel, and choose a color range. The figures below are created using the color scheme Spectrum 2 – 16.

Figure 2.3: Pressure and Mach Number Fields of the Clean NACA0012 at an AoA of 4 Degrees

Create a new FENSAP run within the same project directory and name it: FENSAP_rough_4deg.

To ease data entry, drag & drop the config icon (blue gear) of the initial run (Flow Solution on the Clean NACA0012 Airfoil) onto this new run. The reference parameters of the first run are then automatically copied into the new run.

Double-click the config icon to open the configuration window. In the Surface roughness box choose Specified sand-grain roughness and set the value to 0.0005 m (default). This value has been determined to be a good setting for most icing cases based on the validation data available in the literature. Later in the tutorial, the proprietary ice roughness computation model (beading) of FENSAP-ICE will be used to compute the ice roughness distribution over the airfoil as this ice roughness naturally changes with time. For now, the classical approach of specifying a uniform roughness on the airfoil is followed.

All the other settings regarding the Inlet and Wall boundary conditions, and solver settings are imported from the previous tutorial. Launch the calculation by clicking Run and Start menu.

At convergence, the lift and drag coefficients read 0.42669 and 0.01968 respectively. That is a 7.8% loss in lift and a 105.4% increase in drag. The increase in drag due to roughness is quite high in this case, partly because the roughness height is significant for this size of airfoil (0.5334 m chord) and that the whole surface is set as rough. In reality, only about the first 10% of the chord gets iced in average. For icing calculations, the flow solution should be computed with roughness set everywhere since there is no knowledge of droplet impingement and icing limits a priori.

Figure 2.4: Average Residual and the Total Heat Flux Convergence

Go back to the main project window and set Viewmerical as the default post processor by going to the Settings menu, opening Preferences, switching to Postprocessing tab, and selecting VIEWMERICAL. Click OK.

Back in the project window, right-click the solution icon soln of the first run labeled FENSAP_clean_4deg and choose View with VIEWMERICAL. The program will launch and show an isometric display of the entire grid showing the first solution field which is Density.

Go back to the project window again and repeat the above step for the second run, FENSAP_rough_4deg. A window will open asking if you would like to append the solution in the previously launched Viewmerical window. Click Append.

You should now have both solutions loaded in Viewmerical.

The data sets will be listed on the right, as data-soln and data-soln-1, with the boundaries listed under each item. You can rename the data sets to help with their identification by double-clicking on their names. Change data-soln to clean and data-soln-1 to rough.

To display the data side by side, click the first data set that you just renamed to clean, and choose Horizontal-Left under the Split screen menu.

At this point, the Viewmerical window should look like this:

Go to the Data tab, and click the lock icon next to Shared to apply any changes here to both data sets.

You can choose any data field to display in the data list, which currently shows Density. Switch to Velocity Magnitude instead. Then change the color range to Spectrum 2 – 32.

Align the view angle with the Z-symmetry plane by right-clicking on the 3D axes on the lower left, and choosing Top (Z). Alternatively, you can left-click the Z axis itself, or press 5 on the keyboard.

Zoom in on the front part of the airfoil. You can use Ctrl + left-click to draw a zoom box, or scroll the mouse wheel to zoom in and middle-click to pan. Using the camera icon on the upper left corner, you can take a snapshot of the solution window to capture the following image:

Figure 2.5: Velocity Magnitude Field of the Clean (Left) and Rough (Right) Airfoil at an AoA of 4 Degrees

To use bold fonts for the legend, click   on the top left corner of the window and select Command window; then type BIGFONTS in the command line of 3dview console and hit Enter; now the legend fonts become bold.

For a more in-depth quantitative comparison, you will use 2D data plots. Click the Query tab and enable 2D Plot.

Change the Cutting plane to Z and the horizontal axis to X.

First, you will compare the pressure coefficient on both solutions. On the lower right, there is access to data sets and solution fields. Click the solution field section and choose Pressure Coefficient.

The color and thickness of the data curves displayed in the graph can be changed via the cube menu on the top right and by choosing Curve Settings. Set the first curve red and the second blue, and both curve widths to 2 and click OK.

The Cp graph of both runs should look like the following:

Figure 2.6: Distribution of Pressure Coefficient on the Clean and Rough Airfoil at an AoA of 4 Degrees

You can draw a zoom box by Shift + left-click, and zoom out by middle-clicking your mouse. Left-clicking on the curves will display the value and location. Zooming in before left-clicking on the curves will help. You will see that the pressure coefficient of the rough solution is slightly higher on the suction side (negative Cp), resulting in the aforementioned 7.8% loss in lift.

Next, you will compare the heat fluxes in both cases. Keeping everything the same, switch to Classical Heat Flux in the lower right corner of the window. Also, change the horizontal axis to Y, which better highlights the leading-edge area. Applying roughness of 0.5 mm increases the heat fluxes. This in turn will increase the ice accretion rate in ICE3D. It is crucial that the flow solution for icing is computed with roughness, otherwise the computed ice thickness will be lower and a lot of runback will take place.

Figure 2.7: Distribution of Classical Heat Flux on the Clean and Rough Airfoil at an AoA of 4 Degrees

Create a new run in the same project directory as the previous tutorials by clicking on the new run icon. Select the DROP3D solver, and name it DROP3D_MVD.

Drag & drop the blue config icon of the FENSAP_rough_4deg run (Flow Solution on the Rough NACA0012 Airfoil) onto the config icon of DROP3D_MVD. The input parameters of FENSAP will be automatically copied into DROP3D.

Open the configuration window by double-clicking the DROP3D config icon.

In the Model panel, the default configuration should appear which sets Physical model and Particle type to Droplets, and Droplet drag model to Water – default.

In the Conditions panel, the reference flow conditions imposed for the air solution (FENSAP) should have been copied automatically into DROP3D. In addition to the flow conditions the droplet reference conditions need to be set.

Set the Liquid Water Content (LWC) to 0.55 g/m³ and the Droplet diameter to 20 microns. The Droplet distribution box should read Monodisperse, which means that the diameter that is set represents the MVD of the cloud. For certification purposes, the Appendix C is available in FENSAP-ICE to pick the LWC and droplet size based on the free stream temperature and cloud types. You can temporarily enable Appendix C and click Configure to see the charts and experiment which MVD to get the matching LWC. Once you are done, click Cancel and disable the appendix to return to the original settings.

The Droplet initial solution is based on constant reference values and an inlet velocity AoA of 4 degrees, automatically imported after conversion. Select Velocity angles under Droplet initial solution and set the Angle of attack (X-Y) to 4 deg. The LWC field will be initialized with the reference LWC value throughout the domain. The Dry initialization option sets the initial condition to zero LWC everywhere except the inlets. This option may be needed depending on the type of simulations, especially in the case of internal flows. Leave it unchecked for now.

Go to the Boundaries panel. Select the Inlet boundary and click the Import reference conditions button to automatically set the Inlet conditions. The boundary type should be set to Supersonic or far-field.

Go to the Solver panel. The local time step is computed from the local velocity, drag and the length of each element. Set the CFL number to 20 and Maximum number of time steps to 300. DROP3D uses different time steps than FENSAP since the governing PDEs are different.

Go to the Out panel. Save the solution at every 40 iterations by overwriting the solution file.

Click the Run button at the bottom of the panel to switch to the run window. Run the computation using 4 or more CPUs if possible.

The calculations stop when the convergence level reaches the convergence limits set on the residual and on the total collection efficiency. Otherwise, the simulation continues until DROP3D reaches 300 iterations. In the Graphs panel, look at Average Residual, Total Beta (Collection Efficiency), and Change in total Beta curves. Convergence of water catch is in general achieved when these curves level off at low values of residual. Often the solution in the wake of the droplet flow is still converging while the impingement at the surfaces is fully converged. If you wish to converge the wake and the shadow zones further, Convergence level in the Advanced solver settings of the Solver panel should be reduced. The droplet wake usually is not of interest and it is sufficient to achieve convergence of the total beta alone. However, in some cases like turbomachinery computations, the wake of one stage sets the inlet conditions of the next stage, therefore it is crucial to converge the wake there as well.

Figure 2.8: Average Residual, Total Beta, and Change in Total Beta Curves

View the solution in Viewmerical by clicking the View button. LWC will be available for the whole domain while collection efficiency will only be displayed on the walls. Examine the LWC distribution in the area close to the airfoil, as indicated in the figure below. The blue region is called the shadow zone, where no droplets exist. In between the shadow zone and the free stream, there are bands of high LWC concentrations which are the enrichment zones forming due to the constriction of stream tubes in the continuum domain. These features can be of special interest for aircraft components downstream.

Figure 2.9: LWC at an Angle of Attack of 4 Degrees, Showing the Shadow Zone (Blue Region)

Switch to collection efficiency in the data box, and plot the 2D curve against Y axis in the Query panel. The cutting plane should be Z. The maximum beta occurs at the stagnation point, just below the leading edge in this case. The points on the upper and lower surfaces where beta becomes zero are the impingement limits. In rime icing cases, all the water that impinges is frozen instantly, therefore icing limits are the same as impingement limits. In glaze icing, water can runback and freeze past the impingement limits. Maximum beta is usually no more than 1.0, and reduces as the droplet flow becomes tangent to the surface.

Figure 2.10: Collection Efficiency on the Surface of the Airfoil at an Angle of Attack of 4 Degrees

Create a new DROP3D run within the same project and name it DROP3D_Lang_D.

Drag & drop the blue config icon of the calculation performed in Monodispersed Calculation onto the config icon of this new run.

Double-click the config icon and go to the Conditions panel. Select the Langmuir-D distribution in the Droplet Distribution box of the Droplets reference conditions section. Click the View distribution button.

The droplet diameters are on the horizontal axis, and the weights (the percentage of droplets of a given diameter contained in the cloud) are on the vertical axis. The individual weights are shown with the blue curve, and the overall sum, cumulative weight, is shown with the red curve. On the red curve, the data points are plotted at the mid-range of their cumulative weight intervals. For example, the 20 microns droplet, which happens to be the MVD, covers the cumulative weight range of 35% to 65% and it is therefore plotted at 50% cumulative weight on the red curve.

Keep all the other settings the same, and run the calculation. The individual runs will be executed one after the other, and the results will be combined.

Display all 7 droplet size solutions simultaneously in Viewmerical. First click the View button and choose the first solution Distribution.01/droplet, this will open Viewmerical. Go back to the DROP3D run and click the View button again and choose Append to load the second solution Distribution.02/droplet in the same Viewmerical window. Repeat this operation for all available solutions.

In Viewmerical, go to the Data panel and click Shared. Switch the data field to Collection Efficiency-Droplet. Go to the Query tab, enable 2D plot, and switch the Cutting plane to Z. The graph should display the 7 individual beta distributions. You can draw a zoom box by Shift + left-click, and also you can rename the curves by renaming the original data set names in the Objects panel if you wish.

Figure 2.11: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Langmuir D

The curve with the lowest beta corresponds to the smallest droplet size, and the one with the largest beta corresponds to the largest droplet size. Small droplets are less ballistic, tend to follow the airflow and avoid the aircraft therefore reducing their collection efficiency and impingement limits. Larger droplets are more ballistic and they do not tend to follow the airflow therefore their collection efficiency and impingement are usually higher than the smallest droplets. In general, this information is crucial to properly design the IPS power requirements and coverage.

Using the View button again, choose New Window, this time first load the composite droplet solution that is listed at the end of the drop-down menu, then Append with distribution 4 which is the MVD. Go to Data, click Shared, choose Collection efficiency-Droplet as the data field, and in Query, display the two data sets using the Z Cutting plane. You can name the curves by renaming the data sets in the Objects panel.

Figure 2.12: Comparison of Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Monodisperse vs. Langmuir D

The composite solution is fairly close to that of the MVD. The impingement limits of the composite solution will always be further back due to the inclusion of larger droplets in the distribution. The maximum beta of the composite is lower than the MVD here. This is not always the case. Based on the size and shape of the impingement surface, the composite solution can have a maximum beta that is several times higher than the MVD. In this case, however, the results of the MVD and the distribution are close.

Load the largest and smallest droplet distribution solutions on a new Viewmerical window using New Window and Append, and display them side by side showing the LWC using the procedure learned in Post-Processing Two Solutions with Viewmerical. Observe the difference in the shadow zones. The smallest droplets follow the airfoil very closely but avoiding it while the largest droplets barely change their path and hit almost straight on, leaving a larger shadow zone.

Figure 2.13: LWC Distribution and Shadow Zones for 44.4 Micron Droplets (Left) and 6.2 Micron Droplets (Right)

If not already done, in the FENSAP-ICE main menu, go to Settings → Preferences → Postprocessing and set the Default post-processor software to CFD-Post. Select Write CFD-Post launch files. Click OK to close the window.

Right-click on the MULTI_FENSAPDROPICE run’s config icon and select View previous log/graph. Click View Ice at the bottom of the Execution panel to view the results in CFD-Post.

A View with CFD-Post window will appear. Select -All files- from the drop-down menu and click OK to close the window. CFD-Post will automatically load the ICE3D solutions of every shot.

After opening CFD-Post, a Domain Selector window will request confirmation to load the following domains: ice swimsol, map grid, and map swimsol. Click OK to proceed.

Go to the Calculators tab and double-click on Macro Calculator. The Macro Calculator’s interface panel will be activated and displayed.

Select the Ice Cover – 3D-View macro script from the Macro drop-down list. This will bring up the user interface which contains all input parameters required to view ICE3D output solutions in the CFD-Post 3D Viewer.

The default settings inside the Macro Calculator panel will allow you to automatically output the ice shape of the first shot of the multishot simulation. Output the ice shape at the end of the multishot simulation of Multishot Ice Accretion with Automatic Mesh Displacement, this corresponds to the ice shape of shot 3, by specifying 3 besides the MultiShot Num and then by clicking Calculate. Figure 2.24: Ice View in CFD-Post, Final Ice Shape shows the output the final ice shape.

Figure 2.24: Ice View in CFD-Post, Final Ice Shape

You can also generate and save animations that highlight the ice shape evolution of your multishot simulation. Follow these steps to create and save a custom animation.

Select Ice Cover – 2D-Plot from the Macro drop-down list to create 2D-plots of the multishot simulation. You will create a 2D-Plot that contains all the ice shapes generated by the multishot simulation.

Make sure that Multi-shot # is set to 1.

Change Plot’s Title from default, ICE SHAPE PLOT, to Multishot Ice Shape at -7.5 C (3 shots).

Select Multi-Shots in 2D-Plot (with). The macro will generate a series of 2D plot curves, starting from the assigned shot number in Multi-shot # to the last shot of the simulation.

Inside 2D-Plot (with),

To center the 2D-Plot around the leading edge of the NACA0012, in 2D-Plot (with),

Leave the other default settings unchanged and click Calculate to create a 2D-Plot of the multiple ice shapes in a floating ChartViewer of CFD-Post. Adjust the output window’s size. Figure 2.25: 2D-Plot in CFD-Post, Ice Shapes of the Multishot Simulation shows the output of the macro.

Figure 2.25: 2D-Plot in CFD-Post, Ice Shapes of the Multishot Simulation

The 2D-Plot macro can also export all plotted curves to a CSV format file and simultaneously save the plot as a figure. Keep all input parameters above unchanged and follow these steps.

Create a new run and select FENSAP as flow solver. Name it SWEEP_clean.

Drag & drop the config icon of the run FENSAP_clean_4deg onto the config icon of this new run.

These runs will not be used for icing computations; therefore, you do not need to compute heat fluxes on the walls. Double-click the config icon and proceed to the Boundaries panel. Switch the wall boundary conditions parameters from Temperature to Heat flux, and apply 0 (adiabatic walls). Change all four wall boundary conditions like so. In the Solver panel, increase the Maximum number of time steps to 1000.

Click the Run button and go to the Sweep panel. Enable the sweep mode by choosing BC Inlet – Angle of attack (X-Y). Enter the Minimum and Maximum angles as 0 and 20, and the Number of steps to 11. This will make the step size 2 degrees between the runs.

Start the calculations. A total of 11 flow solutions will be done. This may take some time depending on the Number of CPUs available for the job.

Each run’s graphs and log files will be collected under a different name specifying the angle of attack. In the Graphs panel, you can switch to the Lift Coefficient graph and cycle the solution steps to see the lift for each angle of attack. There is no automatic way of reporting the converged lift values for each step and produce a graph; this step should be done manually. Record the converged value of lift and drag for each angle of attack. The data should look like this:

Create a new run and select FENSAP as the solver. Name it SWEEP_rough.

Drag and drop the config icon of the previous run, SWEEP_clean, onto the config icon of this run.

Double-click the config icon. In the Model panel, set Surface roughness to Sand-grain roughness – BC type. This will enable wall by wall specification of the roughness amount.

Go to the Boundaries panel. Choose BC_2002 which covers the first 20% of the airfoil, and set the Sand-grain roughness in the BC wall parameters section to 0.003 meters. This value is commonly used when performing analyses for icing certification. Visit the other wall boundary conditions and set this value to 0.

Click the Run button, and go to the Sweep panel. Set the Sweep variable to BC Inlet - Angle of attack (X-Y), Minimum and Maximum angles to 0 and 20, and set 11 steps. Start the computations.

Like before, get the converged lift and drag coefficient readings from the lift and drag convergence graphs of each run executed in the sweep.

Create a new run and select FENSAP as the solver. Name it SWEEP_iced.

Drag and drop the config icon of the previous run, SWEEP_clean, onto the config icon of this run.

Switch the grid from the clean naca0012 to the last displaced grid in the multishot run completed in Multishot Ice Accretion with Automatic Mesh Displacement. To do this, double-click the grid icon which currently holds naca0012.grid, navigate to the run directory labeled run_MULTI_FENSAPDROPICE in the project directory, and choose the file grid.ice.000004.

Double-click the config icon. On the Model panel, enable surface roughness and choose the Sand-grain roughness – file option. Click the folder icon to the right and browse to the same multishot run directory. Choose the file named roughness.dat.ice.000003.This is the beading roughness distribution obtained at the last shot by ICE3D.

This grid is too coarse at the leading edge and reducing the CFL or using variable relaxation in the Solver panel is helpful to improve robustness. Activate the Use variable relaxation option and leave the other settings as they are.

Figure 2.26: Leading Edge Detail on the Displaced Grid of the Multishot Run

Click the Run button, and go to the Sweep panel. Set the Sweep variable to BC Inlet - Angle of attack (X-Y), Minimum and Maximum angles to 0 and 20, and set 11 steps. Start the computations.

On the Graphs panel, cycle the solution steps (or angle of attacks) and look at the Residual – Average and Forces - Lift Coefficient graphs to ensure convergence is achieved for each angle of attack. Since the grid is very coarse at the leading edge of the horns, the average residual decreases by only two orders of magnitude at some solution steps (angle of attack of 6 and 8 degrees). For this tutorial, you will consider these results as accurate and converged. However, manual remeshing around the ice is strongly suggested to obtain accurate aerodynamic results.

Like before, get the converged lift and drag coefficient readings from the lift and drag convergence graphs of each run executed in the sweep.

When the data is assembled in the above graphs using Excel, the overall trend appears. Roughness alone is responsible for lowering the lift slope and decreasing the stall angle considerably. The increase in drag due to roughness alone is massive at high angles. With the ice shape included, lift slope and stall angle further decreases, and drag increases even more. Although this example is only showing the effects on a 2D NACA0012 airfoil, similar behavior is expected on 3D wings, helicopter blades, propeller blades, turbo-fan and compressor blades, etc.

Figure 2.27: Lift and Drag Coefficients on a Clean, a Rough, and an Iced Airfoil

After launching FENSAP-ICE, create a new project directory by clicking on the icon below:

Enter the name of the new project directory, FLUENT_ICING, in the Project name box, and browse to position it within your home directory.

A message window will ask about the unit settings. Accept the defaults to keep SI units for this project.

Launch Fluent from your computer. Click Show More Options in the Fluent Launcher window. Under General Options, set your Working Directory to the FLUENT_ICING directory.

Select Dimension as 3D, pick Double Precision under Options, and assign 2 to 4 CPUs under Solver Processes. Press Start menu to close the Fluent Launcher.

Read the case file by going to File → Read → Case. Browse to and select the file ../workshop_input_files/Input_ Grid/Naca0012/naca0012.cas.h5.

From the top bar navigation menu, select Physics → Solver → Operating Conditions.... Set the Operating Pressure (pascal) to 101325 Pa. Press OK.

From the side menu, select General. Ensure the Solver is set to Type: Pressure-Based, Velocity-Formulation: Absolute, and Time: Steady.

From the side menu, select Models → Energy and ensure it is turned on. Then double-click Viscous to open the Viscous Model menu. There are different turbulence models that can be selected. For icing applications using FENSAP-ICE with Fluent, it is strongly recommended to use the popular k-ω SST model. Therefore, change the Model to k-omega (2 eqn) and SST. In the Options section, enable Viscous-Heating and Production Limiter to be consistent with FENSAP. In the Model Constants section, change the Energy Prandtl Number and Wall Prandtl Number to 0.9 and the Production Limiter Clip Factor to 10. Press OK.

From the side menu, click Materials → Fluid and double-click air to open the air properties. Set the Density to ideal-gas. Set the Cp (Specific Heat) to 1004.6882 j/kg.K. This value is equal to 7/2 R air when air is treated as an ideal gas. In FENSAP, the gas constant R is always 287.05376 j/kg.K. Set the Thermal Conductivity to 0.023439363 W/m.K and Viscosity to 1.6801754e-05 Kg/m.s. These values match the previous FENSAP tutorial, In-Flight Icing Using FENSAP Within FENSAP-ICE, and have been computed using the equations presented in the ANSYS FENSAP-ICE User Manual. Click Change/Create to save the air properties, then press Close.

---

Note:  For simplicity, thermal conductivity and viscosity equations presented in the User’s manual are shown below.

where refers to the ambient air static temperature, and , and are equal to 0.00216176 W/m/K^3/2, 288 K and 17.9*10^-6 Pa.s, respectively.

---

From the side menu, right-click pressure-inlet-4 boundary under Boundary Conditions and change the type to pressure-far-field. Double-click it to open and set the far field boundary conditions.

Double-click the wall-5 boundary. In the Momentum panel, set the Shear Condition to No Slip. In the Thermal panel, set the Thermal Conditions to Temperature, and set the Temperature to 280.929174208 K, and press OK. This temperature corresponds to the Adiabatic stagnation temperature + 10 set in FENSAP. Repeat this process for the wall-6, 7, 8 boundaries.

Ensure that symmetry-9 and symmetry-10 boundaries are set to symmetry type.

Select Reference Values from the side menu. Under Compute from, select pressure-far-field-4. Set the Area to 0.05334 m² and the Length to 0.5334 m. These reference values will be used only for postprocessing purposes. For instance, Force coefficients use the reference area, density, and velocity.

Click Solution → Methods on the side menu. Set the Pressure-Velocity Coupling scheme to Coupled. Under Spatial Discretization, set the Gradient to Green-Gauss Node Based and the remaining options to Second or Second Order Upwind.

Double-click Solution → Controls on the side menu. Set the Flow Courant Number to 50. Set Under-relaxation Factors for Turbulent Viscosity and Energy to 0.9.

From the side menu, double-click Monitors → Residuals and modify the Absolute Criteria for convergence to 1e-12 for all parameters. Make sure that the Print to Console and the Plot are enabled and ensure that Monitor Check and Convergence are selected for all parameters.

Monitor the drag and lift coefficients during the simulation.

From the side menu, double-click Solution → Initialization and initialize with Hybrid Initialization. Click Initialize.

Double-click Run Calculation from the side menu. Set the Number of iterations to 500. Click Calculate to start this simulation.

Monitor the convergence of this calculation in the Graphics and the Console windows:

Figure 2.34: Scaled Residuals

Once the simulation is completed, go to Files → Write → Case & Data and save this calculation in the project directory FLUENT_ICING. Name this simulation naca0012_clean.

Launch Fluent from your computer. Click Show More Options. Under General Options, set your Working Directory to the FLUENT_ICING directory.

In the Fluent Launcher window, select Dimension as 3D, pick Double Precision under Options, and assign 2 to 4 CPUs under Solver Processes. Press Start menu to close the Fluent Launcher.

Read the case file by going to File → Read → Case. Browse to and select the file naca0012_clean.cas.h5 created in the previous section.

Double-click the wall-5 boundary. In the Wall Roughness section of the Momentum panel, select High Roughness (Icing) under Roughness Models, choose Specified Roughness under Sand-Grain Roughness, and set Roughness Height (m) to 0.0005 m as a constant value. In the Thermal panel, set the Temperature to 280.929174208 K, and press OK. Repeat this process for wall boundaries wall-6, wall-7, and wall-8.

In this tutorial, reports of lift and drag coefficients convergence are not required. Only plots are sufficient to monitor lift and drag coefficients convergence. To cancel these reports, go to Solution → Monitors → Report Files from the side menu, right-click each report, and Delete them.

Double-click Solution → Initialization from the side menu. Initialize with Hybrid Initialization. Click Initialize.

Double-click Run Calculation from the side menu. Increase the Number of iterations to 1000. Click Calculate.

Once the simulation is completed, go to Files → Write → Case & Data and save this calculation in the project directory FLUENT_ICING. Name this simulation naca0012_rough.

Take a look at the convergence history of this simulation in the Graphics window located at the right of your screen. The following two figures show the convergence of residuals and lift and drag coefficients. You can enlarge and move the legend box in the Graphics windows by dragging one side of or whole the box. Also, examine the convergence of lift and drag coefficients.

Figure 2.37: Scaled Residuals

Figure 2.38: Convergence of Lift and Drag Coefficients of the Rough Airfoil

Create a new run in the project directory FLUENT_ICING by clicking the new run icon. Select the DROP3D solver, and name it DROP3D_MVD.

Right-click the grid icon and select the Define option. Browse and select naca0012_rough.cas.h5. Click Open. When the grid is loaded, the Grid converter dialog window will open to prompt for the selection/verification of the boundary condition types. The data fields will be converted to FENSAP standards.

Use the default settings and click Next. Another window will appear after to double-check the Reference condition.

Press OK.

Use the default settings and click Next. The next window will show details of the process for converting the Fluent case and data files to FENSAP-ICE grid and solution files (.soln, .hflux, .surface, and .params).

Once complete, click Finish. The DROP3D run now has the FENSAP format grid and solution files displayed with Fluent overlaid.

Open the configuration window by double-clicking the DROP3D config icon.

In the Model panel, the default configuration should appear which sets Physical model and Particle type to Droplets, and Droplet drag model to Water – default.

In the Conditions panel, the reference flow conditions imposed for the air solution (Fluent) should have been copied automatically into DROP3D. Verify that the following settings have been copied correctly:

In addition to the flow conditions, the droplet reference conditions should also be set. Therefore, set the Liquid Water Content (LWC) to 0.55 g/m³ and the Droplet diameter to 20 microns. The Droplet distribution box should read Monodisperse, which means that the diameter that is set represents the MVD of the cloud.

For certification purposes, the Appendix C is available in FENSAP-ICE to pick the LWC and droplet size based on the free stream temperature and cloud types. You can temporarily enable Appendix C and click Configure to see the charts and experiment which MVD to get the matching LWC. Once you are done, click Cancel and disable the appendix to return to the original settings.

The Droplet initial solution is based on constant reference values and an inlet velocity AoA of 4 degrees, automatically imported after conversion. Select Velocity angles under Droplet initial solution and set the Angle of attack (X-Y) to 4 deg. The LWC field will be initialized with the reference LWC value throughout the domain. The Dry initialization option sets the initial condition to zero LWC everywhere except the inlets. This option may be needed depending on the type of simulations, especially internal flows. Leave it unchecked now.

Go to the Boundaries panel. Select the Inlet boundary and click the Import reference conditions button to automatically set the Inlet conditions. The boundary type should be Supersonic or far-field.

Go to the Solver panel. The local time step is computed from the local velocity, drag and the length of each element. Set the CFL number to 20 and Maximum number of time steps to 300.

Go to the Out panel. Save the solution at every 40 iterations by overwriting the solution file.

Click the Run button at the bottom of the panel to switch to the run window. Run the computation using 4 or more CPUs if possible.

The calculations stop when the convergence level reaches the convergence limits set on the residual and on the total collection efficiency. Otherwise, the simulation continues until DROP3D reaches 300 iterations. In the Graphs panel, you can look at Average Residual, Total Beta (Collection Efficiency), and Change in total Beta curves.

Figure 2.39: Average Residual, Total Beta, and Change in Total Beta Curves

Convergence of water catch is in general achieved when these curves level off at low values of residual. Often the solution in the wake of the droplet flow is still converging while the impingement at the surfaces is fully converged. If you wish to converge the wake and the shadow zones further, Convergence level in the Advanced solver settings of the Solver panel should be reduced. The droplet wake usually is not of interest and it is sufficient to achieve convergence of the total beta alone. However, in some cases like turbomachinery computations, the wake of one stage sets the inlet conditions of the next stage, therefore it is crucial to converge the wake there as well.

To post process the results of this section, you can skip the next section and go directly to Post-Processing Using Viewmerical.

Create a new DROP3D run within the same project and name it DROP3D_Lang_D.

Drag & drop the blue config icon of the calculation performed in Monodispersed Calculation onto the config icon of this new run.

Double-click the config icon and go to the Conditions panel. Select the Langmuir-D distribution in the Droplet Distribution box of the Droplets reference conditions section. Click the View distribution button.

The droplet diameters are on the horizontal axis, and the weights (the percentage of droplets of a given diameter contained in the cloud) are on the vertical axis. The individual weights are shown with the blue curve, and the overall sum, cumulative weight, is shown with the red curve. On the red curve, the data points are plotted at the mid-range of their cumulative weight intervals. For example, the 20 microns droplet, which happens to be the MVD, covers the cumulative weight range of 35% to 65% and it is therefore plotted at 50% cumulative weight on the red curve.

Keep all the other settings the same, and run the calculation. The individual runs will be executed one after the other, and the results will be combined.

Go to the Graphs panel to monitor the convergence of each run. Look at Average Residual, Total Beta (Collection Efficiency), and Change in total Beta curves to ensure the simulation is converged sufficiently for each run. You can also look at the log file of each run through the Log panel.

Go back to the main project window and set Viewmerical as the default post processor by going to Settings → Preferences → Postprocessing → VIEWMERICAL. Click OK.

Back in the project window, right-click the solution icon droplet of the DROP3D_MVD run and choose View with VIEWMERICAL. The program will launch and show an isometric display of the entire grid showing the first solution field which is Droplet LWC.

You can also see the solution after launching the simulation, directly, by clicking the View button at the bottom of the Graphs panel.

You can choose any field to display in the data list. Go to the Data tab and choose Droplet LWC (kg/m3) to display LWC. Then change the color range to Spectrum 2 – 16.

Align the view angle with the Z-symmetry plane by right-clicking on the 3D axes on the lower left, and choosing Top (Z). Alternatively, you can left-click the Z axis itself.

Zoom in on the airfoil. You can use Ctrl + left-click to draw a zoom box, or scroll the mouse wheel to zoom in and middle-click to pan.

To use bold fonts for the legend, click   on the top left corner of the window and select Command window; then type BIGFONTS in the command line of the 3dview console and hit Enter. The legend fonts now become bold.

Using the camera icon on the upper left corner, you can take a snapshot of the solution window to capture the following image.

Figure 2.40: LWC over NACA 0012 at an AoA of 4 Degrees, Showing the Shadow Zone (Blue Region)

Examine the LWC distribution in the area close to the airfoil. The blue region is called the shadow zone, where no droplets exist. In between the shadow zone and the free stream, there are bands of high LWC concentrations which are the enrichment zones forming due to the constriction of stream tubes in the continuum domain. These features can be of special interest for aircraft components downstream.

Go to the Data tab and choose Collection efficiency-Droplet. Collection efficiency is only displayed on the walls of your geometry. Go to Objects tab and uncheck BC_1000 and BC_4300 to display the collection efficiency distribution only on the walls (BC_2000, BC_2001, BC_2002, and BC_2003).

Use left mouse buttons to rotate, middle mouse button to pan, and right-mouse button to zoom in the airfoil surface to obtain following figure.

Figure 2.41: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees

For a more in-depth quantitative view, it would be possible to create 2D data plots using Viewmerical. Click the Query tab and enable 2D Plot.

Change the Cutting plane to Z and the horizontal axis to Y.

On the lower right corner of Viewmerical, you can directly modify data sets and solution fields. Leave them as they are now.

The color and thickness of the data curve displayed in the graph can be changed by left clicking on the cube menu located on the top right and by choosing Curve Settings. Set the curve color to red and the curve widths to 2 and press OK.

Finally, the following 2D plot is generated.

Figure 2.42: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Monodisperse

The maximum beta occurs at the stagnation point, just below the leading edge in this case. The points on the upper and lower surfaces where beta becomes zero are the impingement limits. In rime icing cases, all the water that impinges is frozen instantly, therefore icing limits are the same as impingement limits. In glaze icing, water can runback and freeze past the impingement limits. Maximum beta is usually no more than 1.0, and reduces as the droplet flow becomes tangent to the surface.

To save data points of this collection efficiency distribution, go to the cube menu on the top right and choose Save one file. A new window pups up to browse and name the file that should contain these data points.

Close the Viewmerical window.

You can also open and compare several solution files using Viewmerical. Let’s display simultaneously all 7 droplet size solutions obtained in DROP3D_Lang_D run.

Back in the project window, right-click the solution icon droplet of the DROP3D_Lang_D run, choose View with VIEWMERICAL, and select the first solution Distribution.01/droplet; this will open Viewmerical. Go back in the project window, right-click the droplet icon of the same run, choose View with VIEWMERICAL, and choose Append to load the second solution Distribution.02/droplet in the same Viewmerical window. Repeat this operation for all the available solutions.

In Viewmerical, go to the Data panel and click Shared. Switch the data field to Collection Efficiency-Droplet.

Go to Query tab, enable 2D plot, and switch the Cutting plane to Z. The graph should display 7 individual beta distributions. You can draw a zoom box by Shift + left-click and you can rename the curves by double-clicking on the original name of each dataset in the Objects panel and entering the new name in the window Rename dataset. In addition, you can change the color and thickness of the data curve displayed in the graph via the cube menu on the top right and by choosing Curve Settings.

Figure 2.43: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Langmuir D

The curve with the lowest beta corresponds to the smallest droplet size, and the one with the largest beta corresponds to the largest droplet size. Smallest droplets are less ballistic, tend to follow the air flow and avoid the aircraft therefore reducing their collection efficiency and impingement limits. Larger droplets are more ballistic and they do not tend to follow the airflow. Therefore, their collection efficiency and impingement are usually higher than the smallest droplets. In general, this information is crucial to properly design the IPS power requirements and coverage.

To compare the composite result to that of the MVD, right-click the solution icon droplet of the DROP3D_Lang_D run and choose View with VIEWMERICAL. Choose New Window and select the composite droplet solution that is listed at the end of the drop-down menu. Go back to the project window and right-click the solution icon droplet of the DROP3D_MVD run, and Append.

Go to Data, click Shared, choose Collection efficiency-Droplet as the data field. Go to Query tab, enable 2D plot, and switch the Cutting plane to Z to display the two data sets. You can name the curves by renaming the data sets in the Objects panel, and change the color and thickness of the data curves using Curve Settings. The following figure shows resultant 2D-plots.

Figure 2.44: Comparison of Collection Efficiency on the Surface, Langmuir D vs. Monodisperse

The composite solution is fairly close to that of the MVD. The impingement limits of the composite solution will always be further back due to the inclusion of larger droplets in the distribution. The maximum beta of the composite is lower than the MVD here. This is not always the case. Based on the size and shape of the impingement surface, the composite solution can have a maximum beta that is several times higher than the MVD. In this case, however, the results of the MVD and the distribution are close.

Load the largest and smallest droplet distribution solutions on a new Viewmerical window using New Window and Append, and display their LWC side by side.

Click the first data set, largest droplet size, and choose Horizontal-Left under Split screen menu.

Go to the Data tab, and click the lock icon next to Shared to apply any changes here to both data sets. Keep LWC data field to display in the data list and change the color range to Spectrum 2 – 16. Align the view angle with the Z-symmetry plane and zoom in to capture the following image:

Figure 2.45: LWC Distribution and Shadow Zones for 44.4 Micron Droplets (Left) and 6.2 Micron Droplets (Right)

Observe the difference in the shadow zones. The smallest droplets follow the airfoil very closely but avoiding it while the largest droplets barely change their path and hit almost straight on, leaving a larger shadow zone.

Now, let’s compare the Fluent airflow results obtained on the clean airfoil, Flow Solution on the Clean NACA0012 Airfoil and on the rough airfoil, Flow Solution on the Rough NACA0012 Airfoil. Viewmerical only reads FENSAP format files, therefore comparisons will be done using the Fluent solution files converted in FENSAP format.

Right-click the airsol icon naca0012_rough.soln of the DROP3D_MVD run, select View with Viewmerical, and choose New Window to display the airflow solution.

In the case of the clean airfoil solution, create a new DROP3D run and name it DROP3D_Clean. Right-click the grid icon and select the Define option to browse and select naca0012_clean.cas.h5, and then follow the same procedure used in step 2 of DROP3D Droplet Impingement on the NACA0012 (Starting from Fluent Airflow), to convert Fluent case and data files into FENSAP-ICE grid and soln files. Once this is done, right-click the airsol icon naca0012_clean.soln of this run and Append this solution in Viewmerical.

In the Objects panel of Viewmerical, rename the first dataset to Fluent-rough and the second dataset to Fluent-clean. Go to Data panel, click Shared and choose Pressure as the data field. Go to Query tab, enable 2D plot, and switch the Cutting plane to Z. The graph should display 2 individual pressure distributions. Change the color and thickness of the data curve using Curve Settings. Set the first curve red and the second blue, and both curve widths to 2 and click OK. You will see that the pressure of the rough solution is slightly higher on the suction side, resulting in reduction of lift (10%). The lift coefficient of the clean airfoil is 0.4480 while the lift coefficient of the rough airfoil is 0.4023.

Figure 2.46: Comparison of Pressure Distributions on the Surface of the Airfoil at an AoA of 4 Degrees, Rough vs. Clean surface

On the lower right corner of your Viewmerical screen, you can modify to solution field. Change Pressure to Classical heat flux to display the following image.

Figure 2.47: Comparison of Classical Heat Flux on the Surface of the Airfoil at an AoA of 4 Degrees, Rough vs. Clean Surface

As shown in the figure above, applying roughness of 0.5 mm increases the heat fluxes. This in turn will increase the ice accretion rate in ICE3D as cooling effects have increased. It is crucial that the flow solution for icing is computed with roughness, otherwise the computed ice thickness will be much smaller and a lot of runback will take place.

After launching FENSAP-ICE, create a new project directory by clicking on the icon below:

Enter the name of the new project directory, CFX_ICING, in the Project name box, and browse to position it within your home directory.

A message window will ask about the unit settings. Accept the defaults to keep SI units for this project.

Launch CFX from your computer to open the CFX launcher window. Set the Working Directory to the CFX_ICING directory.

Select CFX → CFX-Pre from the top menu to launch CFX-Pre.

Import a previously saved .cfx file. Select File → Open Case. Browse to and select the file ../workshop_input_files/Input_Grid/Naca0012/naca0012.cfx.

In the top menu bar, go to Edit, and select Options. In CFX-Pre → General → Beta Options → Physics Beta Features, click to Enable Beta Features.

Double-click Domain 1 under Simulation → Flow Analysis 1 to edit the domain settings.

Set-up boundary conditions inside Simulation → Flow Analysis 1 → Domain 1.

Under Simulation → Flow Analysis 1, double-click Initialization to edit the Global Initialization settings. Set the Velocity Type to Cartesian, with Option set to Automatic with Value with [U, V, W] components of [102.5495844, 7.170965501, 0] m/s. Set the Static Pressure Option to Automatic with Value and the Relative Pressure to 0 Pa. Set the Temperature Option to Automatic with Value and the Temperature to 265.67 K. Set the Turbulence Option to Intensity and Eddy Viscosity Ratio, the Fractional Intensity Option to Automatic with Value with a Value of 0.0008, and the Eddy Viscosity Ratio Option to Automatic with Value with a Value of 1e-5. Click OK.

Under Simulation → Flow Analysis 1, go to Solver and double-click Solver Control. In the Basic Settings tab, set the Advection Scheme and Turbulence Numerics to High Resolution. Set Max. Iterations to 2000. Set the Timescale Control to Local Timescale Factor and the Timescale Factor to 2.0. Set the Convergence Criteria Residual Type and Target to RMS and 1e-20, respectively.

On the side menu bar panel, edit the default Air Ideal Gas properties located under Simulation → Materials.

Right-click the Simulation Control and go to Insert → Execution Control. Inside Execution Control, select Double Precision, and set Number of Processors to 4, if available. Set the Solver Input File to naca0012_clean.def and change the path such that it lies within the CFX_ICING project directory.

Right-click Expressions on the side menu panel, and select Insert → Expression. Name the expression Surface Heat Flux. In the expression Definition panel, enter the equation from the figure below, and click Apply. This equation computes the averaged heat flux, which is important for monitoring the convergence of the simulation.

Save the CFX case by going to Files → Save Case As and save this calculation in the project directory CFX_ICING. Name this simulation naca0012_clean.cfx.

Select Run Solver and Monitor from the top menu to start this simulation. A window may pop up requesting to choose the location of your .def file. Click Rename and save the .def file to the project directory CFX_ICING. Name the file naca0012_clean.def.

The CFX-Solver Manager will open and the simulation will start. Monitor the convergence of this calculation using the graph and the out file monitor windows:

Once the simulation is completed, CFX will save the solution file automatically in the project directory CFX_ICING and its name is naca0012_clean_001.res.

Launch CFX from your computer the open the CFX launcher window. Set the Working Directory to the CFX_ICING directory.

Select CFX → CFX-Pre from the top menu to launch CFX-Pre.

Open the case file by going to File → Open Case. Browse to and select the file naca0012_clean.cfx created in the previous section.

Double-click the Wall_2000 boundary. In the Wall Roughness section of the Boundary Details panel, select High Roughness (Icing) under Option and set Sand-Grain Roughness to 0.0005 m. Click OK. Repeat this process for wall boundaries Wall_2001, and Wall_2002.

On the top menu bar panel, go to Execution Control. Inside Execution Control, select Double Precision, and set Number of Processors to 4, if available. Set the Solver Input File to naca0012_rough.def and change the path such that it lies within the CFX_ICING project directory.

Save the CFX case by going to Files → Save Case As and save this calculation in the project directory CFX_ICING. Name this simulation naca0012_rough.cfx.

Select Run Solver and Monitor from the top menu to start this simulation. A warning window may pop up requesting to confirm the location and name of your .def file. Check to see that it is correctly named naca0012_rough.def and located in the project directory CFX_ICING. Click Continue.

The CFX-Solver Manager will open and the simulation will start. Monitor the convergence of this calculation using the graph and out file monitor windows.

Figure 2.68: Momentum and Mass Residuals

Figure 2.69: Heat Transfer Residual

Figure 2.70: User Points: Wall Heat Flux Monitor

Create a new run in the project directory CFX_ICING by clicking the new run icon. Select the DROP3D solver, and name it DROP3D_MVD.

Right-click the grid icon and select the Define option. Browse and select naca0012_rough_001.res. Click Open. When the grid is loaded, the Grid converter dialog window will open to prompt for the selection/verification of the boundary condition types. CFX data fields will be converted to FENSAP standards.

Leave the rest as default and click Next. Another window will appear to double check the conversion of datafields.

Use the default settings and click Next. Another window will appear to double-check the Reference conditions. The FENSAP-ICE converter will attempt to import reference conditions by reading grid dimensions and boundary conditions from the CFX run.

Press OK. Another window will appear where the reference conditions can be modified if necessary. In this case, the boundary conditions may not have been read precisely. Set the Reference length to 0.5334, the Reference static temperature to 265.67 and the Reference velocity to 102.8. During conversion, all reference conditions use the SI unit system.

Click Next. The next window will show details of the process for (soln, .hlux, and .surface files). A .params file is created that converting the CFX .res file to FENSAP-ICE grid and solution files (.soln, .hflux, and .surface files). A .params file is created that contains the conversion settings.

Once complete, click Finish. The DROP3D run now has the FENSAP format grid and solution files displayed with CFX overlay.

Open the configuration window by double-clicking the DROP3D config icon.

In the Model panel, the default configuration should appear which sets Physical model and Particle type to Droplets, and Droplet drag model to Water – default.

In the Conditions panel, the reference flow conditions imposed to the air solution (CFX) should have been copied automatically into DROP3D. Verify that the following settings have been copied correctly:

In addition to the flow conditions, the droplet reference conditions should also be set. Therefore, set the Liquid Water Content (LWC) to 0.55 g/m³ and the Droplet diameter to 20 microns. The Droplet distribution box should read Monodisperse, which means that the diameter that is set represents the MVD of the cloud.

For certification purposes, the Appendix C is available in FENSAP-ICE to pick the LWC and droplet size based on the free stream temperature and cloud types. You can temporarily enable Appendix C and click Configure to see the charts and experiment which MVD to get the matching LWC. Once you are done, click Cancel and disable the appendix to return to the original settings.

The Droplet initial solution is based on constant reference values and an inlet velocity AoA of 4 degrees, automatically imported after conversion. Select Velocity angles under Droplet initial solution and set the Angle of attack (X-Y) to 4 deg. The LWC field will be initialized with the reference LWC value throughout the domain. The Dry initialization option sets the initial condition to zero LWC everywhere except the inlets. This option may be needed depending on the type of simulations, useful for internal flows or for strong secondary flows. Leave it unchecked now.

Go to the Boundaries panel. Select Inlet_1000 and click the Import reference conditions button to automatically set the Inlet conditions by using the droplets reference conditions and the initial droplet velocity vector.

Go to the Solver panel. The local time step is computed from the local velocity, drag and the length of each element. Set the CFL number to 20 and Maximum number of time steps to 300.

Go to the Out panel. Save the solution at every 40 iterations by overwriting the solution file.

Click the Run button at the bottom of the panel to switch to the run window. Run the computation using 4 CPUs, if possible.

The calculations stop when the convergence level reaches the convergence limits set on the residual and on the total collection efficiency. Otherwise, the simulation continues until DROP3D reaches 300 iterations. In the Graphs panel, you can look at Average Residual, Total Beta (Collection Efficiency), and Change in total Beta curves to monitor convergence.

Figure 2.72: Average Residual, Total Beta, and Change in Total Beta Curves

Convergence of water catch is in general achieved when these curves level off at low values of residual. Often the solution in the wake of the droplet flow is still converging while the impingement at the surfaces is fully converged. If you wish to converge the wake and the shadow zones further, the Convergence level in the Advanced solver settings of the Solver panel should be reduced. The droplet wake usually is not of interest and it is sufficient to achieve convergence of the total beta alone. However, in some cases like turbomachinery computations, the wake of one stage sets the inlet conditions of the next stage, therefore it is crucial to converge the wake there as well.

To post process the results of this section, you can skip the next section and go directly to Post-Processing Using Viewmerical.

Create a new DROP3D run within the same project and name it DROP3D_Lang_D.

Drag & drop the blue config icon of DROP3D_ Lang_D onto the config icon of this new run.

Double-click the config icon and go to the Conditions panel. Select the Langmuir-D distribution in the Droplet Distribution box of the Droplets reference conditions section. Click the View distribution button.

The droplet diameters are on the horizontal axis, and the weights (the percentage of droplets of a given diameter contained in the cloud) are on the vertical axis. The individual weights are shown with the blue curve, and the overall sum, cumulative weight, is shown with the red curve. On the red curve, the data points are plotted at the mid-range of their cumulative weight intervals. For example, the 20 microns droplet, which happens to be the MVD, covers the cumulative weight range of 35% to 65% and it is therefore plotted at 50% cumulative weight on the red curve.

Keep all the other settings the same, and run the calculation. The individual runs will be executed one after the other, and the results will be combined.

Go to the Graphs panel to monitor the convergence of each run. Look at Average Residual, Total Beta (Collection Efficiency), and Change in total Beta curves to ensure the simulation is converged sufficiently for each run. You can also look at the log file of each run through the Log panel.

Go back to the main project window and set Viewmerical as the default post processor by going to Settings → Preferences → Postprocessing → VIEWMERICAL. Click OK.

Back in the project window, right-click the solution icon droplet of the DROP3D_MVD run and choose View with VIEWMERICAL. The program will launch and show an isometric display of the entire grid showing the first solution field which is Droplet LWC.

You can also see the solution after launching the simulation, directly, by clicking the View button at the bottom of the Graphs panel.

You can choose any field to display in the data list. Go to the Data tab and choose Droplet LWC (kg/m3) to display LWC. Then change the color range to Spectrum 2 – 16.

Align the view angle with the Z-symmetry plane by right-clicking on the 3D axes on the lower left, and choosing Top (Z). Alternatively, you can left-click the Z axis itself.

Zoom in on the airfoil. You can use Ctrl + left-click to draw a zoom box, or scroll the mouse wheel to zoom in and middle-click to pan.

To use bold fonts for the legend, click   on the top left corner of the window and select Command window; then type BIGFONTS in the command line of the 3dview console and hit Enter. The legend fonts now become bold.

Using the camera icon on the upper left corner, you can take a snapshot of the solution window to capture the following image.

Figure 2.73: LWC over NACA0012 at an AoA of 4 Degrees, Showing the Shadow Zone (Blue Region)

Examine the LWC distribution in the area close to the airfoil. The blue region is called the shadow zone, where no droplets exist. In between the shadow zone and the free stream, there are bands of high LWC concentrations which are the enrichment zones forming due to the constriction of stream tubes in the continuum domain. These features can be of special interest for aircraft components downstream.

Go to the Data tab and choose Collection efficiency-Droplet. Collection efficiency is only displayed on the walls of your geometry. Go to the Objects tab and uncheck BC_1000 and BC_4300 to display the collection efficiency distribution only on the walls (BC_2000, BC_2001, BC_2002, and BC_2003).

Use the left mouse button to rotate, middle mouse button to pan, and right mouse button to zoom into the airfoil surface to obtain following figure.

Figure 2.74: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees

For a more in-depth quantitative view, it would be possible to create 2D data plots using Viewmerical. Click the Query tab and enable 2D Plot.

Change the Cutting plane to Z and the horizontal axis to Y.

On the lower right corner of Viewmerical, you can directly modify data sets and solution fields. Leave them as they are for now.

The color and thickness of the data curve displayed in the graph can be changed by left clicking on the cube menu located on the top right and by choosing Curve Settings. Set the curve color to red and the curve widths to 2 and press OK.

Finally, the following 2D plot is generated.

Figure 2.75: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Monodisperse

The maximum beta occurs at the stagnation point, just below the leading edge in this case. The points on the upper and lower surfaces where beta becomes zero are the impingement limits. In rime icing cases, all the water that impinges is frozen instantly, therefore icing limits are the same as impingement limits. In glaze icing, water can runback and freeze past the impingement limits. Maximum beta is usually no more than 1.0, and reduces as the droplet flow becomes tangent to the surface.

To save data points of this collection efficiency distribution, go to the cube menu on the top right and choose Save one file. A new window appears to browse and name the file that should contain these data points.

Close the Viewmerical window.

You can also open and compare several solution files using Viewmerical. Let’s display simultaneously all 7 droplet size solutions obtained in the DROP3D_Lang_D run.

Back in the project window, right-click the solution icon droplet of the DROP3D_Lang_D run, choose View with VIEWMERICAL, and select the first solution Distribution.01/droplet; this will open Viewmerical. Go back in the project window, right-click the droplet icon of the same run, choose View with VIEWMERICAL, and choose Append to load the second solution Distribution.02/droplet in the same Viewmerical window. Repeat this operation for all the available solutions.

In Viewmerical, go to the Data panel and click Shared. Switch the data field to Collection Efficiency-Droplet.

Go to Query tab, enable 2D plot, and switch the Cutting plane to Z. The graph should display 7 individual beta distributions. You can draw a zoom box by Shift + left-clicking. You can rename the curves by double-clicking on the original name of each dataset in the Objects panel and then by entering the new name in the window Rename dataset. In addition, you can change the color and thickness of the data curve displayed in the graph via the cube menu on the top right and by choosing Curve Settings.

Figure 2.76: Collection Efficiency on the Surface of the Airfoil at an AoA of 4 Degrees, Langmuir D

The curve with the lowest beta corresponds to the smallest droplet size, and the one with the largest beta corresponds to the largest droplet size. Smallest droplets are less ballistic, tend to follow the air flow and avoid the aircraft therefore reducing their collection efficiency and impingement limits. Larger droplets are more ballistic and they do not tend to follow the airflow. Therefore, their collection efficiency and impingement limits are usually higher than the smallest droplets. In general, this information is crucial to properly design the IPS power requirements and coverage.

To compare the composite result to that of the MVD, right-click the solution icon droplet of the DROP3D_Lang_D run and choose View with VIEWMERICAL. Choose New Window and select the composite droplet solution that is listed at the end of the drop-down menu. Go back to the project window and right-click the solution icon droplet of the DROP3D_MVD run, and Append.

Go to Data, click Shared, choose Collection efficiency-Droplet as the data field. Go to Query tab, enable 2D plot, and switch the Cutting plane to Z to display the two data sets. You can name the curves by renaming the data sets in the Objects panel, and change the color and thickness of the data curves using Curve Settings. The following figure shows the resultant 2D-plots.

Figure 2.77: Comparison of Collection Efficiency on the Surface, Langmuir D vs. Monodisperse

The composite solution is fairly close to that of the MVD. The impingement limits of the composite solution will always be further back due to the inclusion of larger droplets in the distribution. The maximum beta of the composite is lower than the MVD here. This is not always the case. Based on the size and shape of the impingement surface, the composite solution can have a maximum beta that is several times higher than the MVD. In this case, however, the results of the MVD and the distribution are close.

Load the largest and smallest droplet distribution solutions on a new Viewmerical window using New Window and Append, and display their LWC side by side.

Click the first data set, largest droplet size, and choose Horizontal-Left under Split screen menu.

Go to the Data tab, and click the lock icon next to Shared to apply any changes here to both data sets. Keep LWC data field to display in the data list and change the color range to Spectrum 2 – 16. Align the view angle with the Z-symmetry plane and zoom in to capture the following image.

Figure 2.78: LWC Distribution and Shadow Zones for 44.4 Micron Droplets (Left) and 6.2 Micron Droplets (Right)

Observe the difference in the shadow zones. The smallest droplets follow the airfoil very closely but avoiding it while the largest droplets barely change their path and hit almost straight on, leaving a larger shadow zone.

Now, let’s compare the CFX airflow results obtained on the clean airfoil, Flow Solution on the Clean NACA0012 Airfoil and on the rough airfoil, Flow Solution on the Rough NACA0012 Airfoil. Viewmerical only reads FENSAP format files, therefore comparisons will be done using the CFX solution files converted in FENSAP format.

Right-click the airsol icon naca0012_rough.soln of the DROP3D_MVD run, select View with Viewmerical, and choose New Window to display the airflow solution.

In the case of the clean airfoil solution, create a new DROP3D run and name it DROP3D_Clean.

Right-click the grid icon and select the Define option to browse and select naca0012_clean_001.res, and then follow the same procedure used in step 2 of the DROP3D Droplet Impingement on the NACA0012 (Starting from CFX Airflow), to convert the CFX .res files into FENSAP-ICE grid and soln files. Once this is done, right-click the airsol icon naca0012_clean_001.soln of this run and Append this solution in Viewmerical.

In the Objects panel of Viewmerical, rename the first dataset to CFX-rough and the second dataset to CFX-clean. Go to Data panel, click Shared and choose Pressure as the data field. Go to Query tab, enable 2D plot, and switch the Cutting plane to Z. The graph should display 2 individual pressure distributions. Change the color and thickness of the data curve using Curve Settings. Set the first curve red and the second blue, and both curve widths to 2 and click OK. You will see that the pressure of the rough solution is slightly higher on the suction side.

Figure 2.79: Comparison of Pressure Distributions on the Surface of the Airfoil at an AoA of 4 Degrees, Rough vs. Clean Surface

On the lower right corner of your Viewmerical screen, you can modify to solution field. Change Pressure to Classical heat flux to display the following image.

Figure 2.80: Comparison of Classical Heat Flux on the Surface of the Airfoil at an AoA of 4 Degrees, Rough vs. Clean Surface

As shown in the figure above, applying roughness of 0.5 mm increases the heat fluxes. This in turn will increase the ice accretion rate in ICE3D as cooling effects have increased. It is crucial that the flow solution for icing is computed with roughness, otherwise the computed ice thickness will be much smaller and a lot of runback will take place.