Create a new project directory by clicking on the New project icon. Name it FENSAP-TURBO. Select the metric unit system.

Create a new FENSAP-TURBO run by clicking on the new run icon and name it TURBO-FLOW.

Set the Number of rows to 5 in the configuration window. This identifies the number of grid files that will be used in this simulation.

Download the 8_FENSAP_ICE_Turbo_Advanced.zip file here .

Unzip 8_FENSAP_ICE_Turbo_Advanced.zip to your working directory.

To assign the grid files, right-click each grid icon, selecting the option Define.

Alternatively, grid files can be assigned within the main input parameter window. Grid files axisymmetric-nacelle, nosecone_ext, fan, bypass and igv are provided in the tutorials subdirectory /workshop_input_files/Input_Grid/Turbo.

Double-click the main config icon at the top of the drop-down list to open the input parameter window.

Go to the Turbo panel. If grid files have not already been assigned to each row, they may be specified here by clicking on the blue browse icon on the right of each row number. Navigate to the default file directory containing the turbofan grid files and arrange them as follows:

Activate rotation for Row 2 and Row 3 by clicking in the check-mark box of each Rotating component, in this case nosecone_ext and fan. Specify a rotation rate of 1380 revolutions per minute (rpm).

The rotation axis must be selected in the Rotation menu at the bottom of the window. In this case, set it to X-axis. Periodic repetitions can be viewed by activating the Display box.

Each component is linked to its adjacent one through the interfaces listed in the Interfaces section. An interface consists of two communicating boundaries; an Exit and an Inlet.

Figure 8.2: Interfaces Are Required to Transfer Boundary Conditions Blue: Nacelle-Nosecone, Purple: Nosecone-Fan, Orange: Fan-Bypass, Green: Fan-Igv

The   button adds an interface. The   button removes an interface. The four interfaces shown in the figure above must be set up. When an interface is added, the pair of rows to be interfaced with their corresponding inlet and exit boundary labels is displayed.

Identify the interface boundaries by opening the grid file of each row (for example, right-click the grids and select View with VIEWMERICAL). Once you know which Exit is connected to which Inlet, fill out the Interfaces table.

Figure 8.3: Interfaces Boundaries Connecting Exits of Each Row to Inlets of Communicating Rows

The relaxation factor for pressure, Press relax., is set to 1 for all interfaces. This value can be decreased if convergence problems occur in complex cases.

Go to the Model panel. The physical model is set to Air. Select the Navier-Stokes option for the Momentum equations and Full PDE for the Energy equation.

Select the Spalart-Allmaras turbulence model. In turbomachines, the eddy to laminar viscosity ratio can be increased to 10. In this tutorial, since the main Inlet is at Far-field, set the Eddy/Laminar viscosity ratio as 1e-5. Keep the Relaxation factor and Number of iterations at 1.

Enable the Surface roughness by using a specified sand-grain roughness of 0.00025 m.

Go to the Conditions panel. In the Reference conditions section, change the following values:

In the Initial Solution section, use Velocity components and set the following values for the three velocity components:

Go to the Boundaries panel. The boundary labels for each component are listed under their corresponding rows. Each boundary that has been already linked to an interface has been disabled. Wall boundary labels may change color according to the type of simulation being run, or based on the selected options.

Set the boundary conditions for each of the following rows:

Row #1

For the inlet boundary BC_1000, choose the inlet type as Subsonic, and click Import reference conditions to assign 110 m/s for Velocity X and 256 K for static Temperature.

Select BC_3000 and choose Subsonic under the Type menu. Click Import reference conditions to assign an exit pressure of 64463.341 Pa. Radial equilibrium must be disabled.

Select BC_5001 and choose Periodic1 in its Turbo part menu. Select BC_5002 and choose Periodic2 in its Turbo part menu.

Row #2

This row contains a wall boundary that rotates (BC_2003) and another that remains stationary (BC_2002). In a rotating frame of reference, stationary walls require a counter-rotating velocity.

To enable the counter-rotating velocity, select the wall boundary BC_2002 and choose the Counter-rotating option in the Rotation dialog box. The boundary label BC_2002 will turn yellow.

Select BC_5003 and choose Periodic1 in its Turbo part menu. Select BC_5004 and choose Periodic2 in its Turbo part menu.

Row #3

Apply the Counter-rotating option to the fan shroud (BC_2005) and the splitter (BC_2007). Use the following table to define the Turbo part of each boundary surface.

Row #4

A radial equilibrium condition accounts for the pressure gradient caused by the tangential velocity.

Select BC_3004, choose Subsonic under the Type menu, and enable Radial Equilibrium option from the drop-down dialog box. Set the following parameters:

Use the following table to define the Turbo part of each boundary surface.

Row #5

Select BC_3005, choose Subsonic under the Type menu, and enable the Radial Equilibrium option from the drop-down dialog box. Set the following parameters:

Use the following table to define the Turbo part of each boundary surface.

All Rows

Adiabatic walls should be imposed on all walls. The heat transfer coefficients required for icing calculations will be calculated as part of the EID computation at the ICE3D step. Navigate through each wall boundary of the components. In the BC wall parameters section, select Heat flux and specify a value of 0 W/m².

Go to the Solver panel. Set the CFL number to 100. Set the Maximum number of time steps to 500. Activate the Use variable relaxation option and keep the default number of Time steps, 300. Keep the default Artificial viscosity settings.

Double-click the Advanced solver settings menu and verify that the convergence level is set to 1e-23. This low value of convergence ensures that a single row does not end prematurely after reaching its convergence limit and stops sending or receiving information from adjoining interfaces.

Go to the Out panel and save the Solution every 40 iterations. Choose to Overwrite the solution file.

Click Run. In the Settings panel, set the Number of CPUs to 5 or higher in the Execution Settings section.

Click the Start button to start the calculation.

Create a new project directory by clicking on the New project icon. Name it FENSAP-TURBO-ICING. Select the Metric unit system from Settings → Units.

Create a new DROP3D-TURBO run by clicking on the new run icon and name it TURBO-MIXED-PHASE-VAPOR. Set the Number of rows to 6.

To import an existing Ansys CFX *.res file, right-click the row 1 grid icon and select the Define option.

The .res solution file solved earlier in Flow Setup in Ansys CFX is located in the /workshop_input_files/Input_Grid/Turbo/CFX sub-directory. This will launch the grid and solution converter for Ansys CFX. For multiple stage turbo applications, click the Multiple grids option.

In general, the boundary conditions conversion options do not need to be modified as this is automatically detected. Verify that the boundary definitions between Ansys CFX and FENSAP are equivalent and click Next.

The Datafields verifies the correspondence between Ansys CFX and FENSAP equivalent data fields. Click Next.

It is necessary to provide all the information requested regarding the reference freestream conditions. These reference values are eventually used in ICE3D to define icing conditions. Currently, some of the reference parameters are not automatically set correctly and you would need to enter these manually as follows:

The Execution and Log panel provide details of the conversion process. Any errors in conversion can be identified in these panels. Click on Finish to complete the CFX to FENSAP-ICE conversion process.

Click Yes when prompted to use the solver settings for reference conditions. The run layout should be populated as shown in the next figure:

Double-click the main config icon of the run to open the DROP3D-TURBO input parameter window.

Under the Turbo panel, you can set up the individual Components and the Interfaces. The graphical user interface automatically sets up each component with its rotational speed and interface connections. The row number corresponds to the order in which Ansys CFX simulation was setup.

---

Note:  To easily identify and view the rows in the graphical window you can click   next to the Components menu.

This will highlight the row that you select and the icon ( ) should turn blue as shown below:

---

Rotation for Row 3 and Row 4 are automatically activated by the check-mark next to the row number. The rotation speeds for Row 3 and Row 4 are 8,000 and 3,800 RPM respectively.

The rotation axis must be selected in the Rotation menu at the bottom of the window. In this case, set it to Z-axis. Periodic repetitions maybe viewed by activating the Display box.

Next, the interfaces should be checked and setup. The figure below shows the interfaces of this current case.

Figure 8.7: Interfaces: Blue: Nacelle-Fan, Orange: Fan-Bypass, Green: Fan-IGV, Cyan: IGV-Rotor, Purple: Rotor-Stator

Each component is linked to its adjacent one through the interfaces listed in the Interfaces section of the Turbo panel. An interface consists of two communicating boundaries; an Exit and an Inlet. The   button adds an interface. The   button removes an interface.

Figure 8.8: Interfaces Boundaries Connecting Exits of Each Row to Inlets of Communicating Rows

---

Note:  To easily identify and view the interfaces in the graphical window you can click   next to the Interfaces menu.

On selecting the row and/or BC the domain and/or interface gets highlighted in the graphical window and the icon ( ) should turn blue as shown below:

---

Proceed to the Model panel.

Set the following under the Particle parameters section:

Set the Vapor model to Enabled and keep the Turbulent Schmidt number to its default value of 0.7. Vapor that exists in the cloud together with the droplets that evaporate and crystals that melt/sublimate will be tracked using this model.

Next, go to the Conditions panel and, under the Reference conditions, verify the following values:

Under the Droplets reference conditions section, set the Droplet diameter to 20 microns. Water density remains at 1,000 kg/m³.

Under the Ice crystal reference conditions section, select Appendix D in the Choose Appendix list and click the Configure button. This will allow you to set either LWC or ICC based on the current Appendix D guidelines for total water content (TWC).

The dark blue curve indicates the temperature isoline based on the reference conditions of this run. The Altitude and Air Temperature are used to determine the allowable values of Total Water Content from the y-axis of the graph. In this case, the total water content is 4.184 g/m³. The Liquid Water Content or the Ice Crystal Content can be changed to add up to the total water content. Set the Liquid Water Content to 0.5 g/m³. Click the OK button to confirm these settings.

Set the following ice crystal properties:

Go to the Particle initial conditions section and choose Velocity components. Set the following:

Set the Vapor initial conditions to Relative humidity and set a value of 100%.

The only boundary condition that is required to be set is the inlet of the first stage. In this case, apply the following boundary conditions for Row05:

For the Inlet boundary BC_1000, select Type Supersonic or far-field, and click Import reference conditions to assign values for Liquid water content, Ice Crystal content, Temperature, Velocity and Relative humidity. Repeat this step for BC_1001.

To post-process the droplet and crystal solutions of this simulation with CFD-Post Turbo, use the following table to set the Turbo part of each boundary surface if this has not been done.

Go to the Solver panel. Set the CFL number to 10 and the Maximum number of time steps to 300. In the Advanced solver settings section set Mass deficit cutoff to 0.1%:

Go to the Out panel, and select the option to save the Solution every 40 iterations. Choose to Overwrite the solution file.

Click the Run button to open the Execution environment. In the Settings panel, set the total number of CPUs to the maximum available. All CPU’s will be used to run each row sequentially. Click the Start menu button to start the calculations.

The ICC enrichment and shadow zones in the figure below show the importance of including the nacelle and spinner components. The resulting shadow zone and the enrichment region surrounding it propagate downstream towards the fan and further down into the bypass. Results below show that only a small amount of ICC/LWC enters the compressor core. In Impingement, Icing and Shedding on Rotating and Stationary Blades, it will be seen that reinjection must be considered to have a more accurate representation of the ICC and LWC concentration and subsequently icing.

In the following two sections, the turbofan air flow, droplet, crystal and vapor solutions will be used to calculate the accretion of ice. Different models within the ICE3D-TURBO framework will be activated, and their corresponding impact on the ice growth will be compared.

Create a new ICE3D-TURBO run by clicking on the new run icon and name it MIXED-PHASE-VAPOR-ICING. Set the number of rows to 6.

Drag & drop the config icon of the previous TURBO-MIXED-PHASE-VAPOR run onto the config icon of this MIXED-PHASE-VAPOR-ICING run to copy the reference and boundary conditions from the airflow, droplet, ice crystal and vapor solutions.

Double-click the config icon to open the ICE3D-TURBO input parameters window. Go to the Model panel and under the Icing model section, select the Glaze-Advanced model, select the Classical heat flux option, and click the Concavity Fix check box. The value should be set to 70 deg.

Enable EID calculation:

For icing simulations where air flow is calculated using adiabatic walls, heat transfer coefficients cannot be readily calculated as a function of the surface convective heat flux and the temperature difference between reference and user-set wall temperatures. EID step uses a propriety technology to post-process the provided adiabatic airflow solution to generate the heat transfer coefficients that are inputs to ICE3D energy equation.

Verify that the Ice crystals dialog box is set to Droplet+Crystals. Activate crystal bouncing by choosing one of the bouncing Modes from the pull-down menu. For this run, select the NTI Bouncing Model. Check the Crystal erosion option to allow crystals to erode the ice layer.

Next, go to the Conditions panel. The Reference conditions are already set and taken from the TURBO-MIXED-PHASE-VAPOR simulation.

---

Note:  The Relative humidity is a factor that determines the surface film evaporation rates. In this setup, the Relative humidity is greyed out, since the vapor transport solution on the surface obtained from the previous section provides the necessary information for ICE3D to calculate accurate film evaporation rates.

---

Go to the Boundaries panel. Each wall contains additional options for icing that include: Disabled, Enabled, Disabled-Sliding, Enabled-Sliding and Sink.

To post-process the icing solutions of this simulation with CFD-Post, use the following table to set the Turbo part of each boundary surface if this has not been done.

Go to the Solver panel. Deselect the Automatic time step option and set the Time step to 1e-4. Set the Total time of ice accretion to 5 seconds. Go to the Out panel and set the Time between solution output to 5 seconds.

Click the Run button at the bottom of the panel to go to the Execution environment. In the Settings panel, set the total Number of CPUs to the maximum available. Click the Start menu button to start the calculation.

The figure below was generated with CFD-Post using the View Ice button in the run panel and by loading all the solution components (all rows). This image shows the ice accretion locations on the 6-stage Turbofan. See the next section for more details of how to post-process turbo icing solutions with CFD-Post.

Figure 8.9: Ice Accretion on All Turbofan Components. Miniature Ice Growth inside Core (Rotor and Stator)

Create a new ICE3D-TURBO run by clicking on the new run icon and name it MIXED-PHASE-ICING. Set the number of rows to 6.

Drag & drop the config icon of the previous MIXED-PHASE -VAPOR-ICING run onto the config icon of this MIXED-PHASE-ICING run to copy the reference and boundary conditions from the airflow, droplets, ice crystal and vapor solutions.

Right-click the config icon of the run and go to the Options menu. Deselect the Use Vapor Solution option to remove the vapor pressure solution from the ICE3D calculation.

Double-click the config icon to open the ICE3D-TURBO input parameters window. First go the Model panel and enable EID:

Next, go to the Conditions panel. In the Model parameters section, set the Relative humidity value to 100%.

---

Note:  The Relative humidity is a factor that determines the surface film evaporation rates. In this setup, since the vapor solution is not being used, the Relative humidity needs to be set. The default value for Relative humidity is 100%.

---

Click the Run button at the bottom of the panel to go to the Execution environment. In the Settings panel, set the total Number of CPUs to the maximum available. Click the Start menu button to start the calculation.

The impact of using constant relative humidity can be seen by comparing the ICE3D solution to the previous run as seen in the figure below. Relative humidity affects the evaporation rates, and consequently the film thickness and ice accretion rate.

Figure 8.15: Film Thickness at 100 % Relative Humidity (Left) and Calculated Relative Humidity (Right)

Create a new ICE3D-TURBO run by clicking the new run icon and name it MIXED-PHASE-ICE-SHEDDING. Set the number of rows to 6.

Drag & drop the config icon of the previous MIXED-PHASE-VAPOR-ICING run onto the config icon of this MIXED-PHASE-ICE-SHEDDING run to copy the reference and boundary conditions from the airflow, droplets, ice crystal and vapor solutions.

Double-click the config icon to open the input parameters window. Go to the Model panel and find the Ice shedding model section. The Delamination and Cracking option uses knowledge of both adhesive strength and cohesive strength in the ice to evaluate if centrifugal forces are sufficient to both delaminate, crack and shed the ice. This model couples the icing calculation to a crack propagation calculation and therefore requires knowledge of the material properties of ice.

In the ice shedding section, set the following options for Ice-Surface interface and Crack detection criteria:

The ice surface interface option determines the adhesive strength between the material and the ice. Ice-Aluminum is chosen for this simulation. The crack detection criteria decide the mechanism used to propagate the crack in the ice. The principal stress criteria are used because the ice is a brittle material and therefore appropriate compared to fracture toughness, which is typically used to materials which can undergo plastic deformation.

The ice materials properties needed for the stress calculation are Youngs modulus, Poisson’s ratio, material density, Principle tensile and shear cohesive strength and fracture toughness. The values listed are all based on average values found in the literature:

Proceed to the Boundaries panel and enable the Icing option for all wall boundaries. Only rotating walls will experience shedding during the accretion time.

The settings for each boundary wall are listed below:

Go to the Solver panel. Deselect the Automatic time step option and set the Time step to 1e-4. Set the Total time of ice accretion to 10 seconds.

Go to the Out panel and set the Time between solution output to 1 seconds. The shedding evaluation is triggered by default at every printout. Change the Shedding every printouts to 4 to trigger shedding at 4 and 8 seconds.

Set Yes for Numbered output files to allow side by side comparisons of the ice before and after shedding.

Click the Run button at the bottom of the panel to go to the Execution environment. In the Settings panel, set the total Number of CPUs to the maximum available. Click the Start menu button to start the calculation.

The shedding calculation produces the following additional files that provide information on ice that sheds:

Load the following blade(row04) solutions in a Viewmerical post-processing window:

To view a side-by-side comparison highlight the ice-shed solution object and choose the Split screen option Horizontal-Right in the drop down menu.

Go to the Data tab. Select Ice thickness as the data field. Set the Color range to Grayscale 64. Click the Shared lock icon to share data between solution objects. Set the global range from -1e-5 to 1e-5.

The image you obtain should look like the one below. The comparison shows that a significant portion of the ice on the blade pressure side has shed.

A closer look at what the delaminated ice pieces are can be seen by observing the shed-ice fragment ID field in the shed-ice solution object. Change the displayed data from Ice Thickness to Shed Ice Fragment ID. Switch the color range to Spectrum 2-32 and enable IsoValues in Enabled – Surfaces mode with Number set to 32. You will get an image like the one below. Enabling iso surfaces clearly outlines the individual fragments. The left image shows the fragments from the 1st analysis at 4 seconds, and the right image shows new fragments that shed at 8 seconds. In the first shedding analysis, 91 fragments were identified, and in the second analysis, 61 additional fragments were identified.

A better understanding of the mass associated with these fragments can be seen in the iceshed.log.row04 file.

A review of the log shows that the largest total amount of ice shed occurred at 8 seconds and the largest piece of ice that shed weighed 0.009 kg and occurred at a radial coordinate of 0.37 m.

The largest piece shed can be viewed in Viewmerical. To do this, click View then select max_ice_shed_piece.grd to load max_ice_shed_piece.grd.row04 alongside the surface mesh-map.grid.row04. Once both grid objects are loaded, change the Cell color of the shed piece to light blue and change the Object type of both grid objects to Smooth shaded:

The following two figures show the comparison of P1 principal stress field before (Left) and after (Right) shedding has taken place at t=4s over the fan (row04). The negative stress indicates regions of ice that are in compression along the blade surface. These are located near the extremity of the ice layer. More continuous principal stress distributions could have been achieved with a finer mesh.

Create a new DROP3D-TURBO run by clicking on the new run icon and name it CRYSTAL-COMPLETE-REINJECTION. Set the number of rows to 6.

Drag and drop the config icon of the previous TURBO-MIXED-PHASE-VAPOR run configured in Droplet and Ice Crystal Impingement to the config icon of this run. This operation copies previous settings and Reference conditions on to this run.

Double-click the config icon to open the DROP3D-TURBO input parameters window and go to the Model panel.

In the Particle parameters section, set the particle type to Crystals.

Make sure the Particle thermal equation and Vapor model options are enabled.

Under the Particle parameters section, enable Particle reinjection by choosing the Complete mode.

The Complete mode uses ICE3D-TURBO to compute the mass of water film detaching from the surface, as well as the mass of crystals that bounced from the surfaces. This mode uses the information from ICE3D-TURBO to perform a simulation of the reinjected particles, using the wall boundaries as inlets, with DROP3D-TURBO.

Each enabled wall family must be subdivided into sections from which particles are reinjected. The walls are divided by axial cuts and are spaced uniformly. To limit the computational time, it is lowered to 2.

In the Particle reinjection section, assign the Number of subdivisions to 2 and their Spacing to Uniform.

In the Icing model section, verify that the following settings are assigned:

Enable EID calculation to post process the adiabatic airflow solution for surface heat transfer coefficient computations:

In the Ice crystals section, double-click Ice crystals to open additional settings. Choose the NTI Bouncing model in Modes and activate Crystals Erosion by clicking on the check box.

Next, go to the Conditions panel. The Reference conditions are already set and taken from the TURBO_MIXED_PHASE simulation.

Under the Ice crystal reference conditions section, select Appendix D in the Choose Appendix list and click the Configure button. This will allow you to set the total water content (TWC) based on the current Appendix D guidelines. Set the remaining ice crystal conditions as follows:

Under the Model parameters section for icing, change the Relative humidity to 100%.

Under the Crystal initial solution section, disable Dry initialization.

Go to the Boundaries panel. Check to see if the walls for each row are as follows:

Additionally, define the inlet conditions of the first stage (nacelle intake). In this case, apply the following boundary conditions for Row05:

For the Inlet boundary BC_1000, click Import reference conditions to assign values for Ice Crystal Content, Temperature and Velocity. Repeat this step for BC_1001.

Go to the Solver panel. The Time integration – Droplets settings are already defined from the previous configuration settings. Ensure the CFL number is set to 10 and the Maximum number of time steps is set to 600. For the Icing simulation settings, open the Time settings dialog window by double-clicking it. Disable the Automatic time step option. Set the Time step to 1e-5 and the Total time of ice accretion to 5 seconds.

Go to the Out panel and set Solution every to 40 iterations. Choose to Overwrite the solution file. For icing file outputs, enter 5 seconds under Time between solution output and select No under Numbered output files.

Click the Run button at the bottom of the panel to go to the Execution environment. In the Settings panel, set the total Number of CPUs to the maximum available. Click the Start menu button to start the calculation.

A comparison of the ICC obtained using 4.18453937 g/m³ ice crystals with and without reinjection (not set up in this tutorial is shown in both figures below.

Figure 8.16: ICC Comparison Without (Left) and With (Right) Reinjection Due to Bouncing

Figure 8.17: ICC Comparison in Bypass, Core IGV and Stator Without (Left) and With (Right) Reinjection

Create a new DROP3D-TURBO run by clicking on the new run icon and name it FILM-COMPLETE-REINJECTION.

Set the number of rows to 6.

Drag and drop the config icon of the TURBO-MIXED-PHASE-VAPOR run configured in Droplet and Ice Crystal Impingement to the config icon of this run. This will copy the previous settings and Reference conditions to this run.

Double-click the config icon to open the DROP3D-TURBO input parameters window and go to the Model panel. Under the Particle parameters section, change the Particle type to Droplets.

Set the Particle reinjection mode to Complete.

The Complete mode for droplet particles only uses ICE3D-TURBO to compute the mass of water film detaching from the surface. This mode uses the information from ICE3D-TURBO to perform a simulation of the reinjected film, using the wall boundaries as inlets, with DROP3D-TURBO.

Set the Vapor model to Disabled.

Each enabled wall family must be subdivided into sections from which particles are reinjected. The walls are divided by axial cuts and are spaced uniformly. To limit the computational time, the number of sub-divisions is lowered to 1.

Under the Particle reinjection section, assign the Number of subdivisions to 1 and their Spacing to Uniform.

Enable EID calculation to post process the adiabatic airflow solution for surface heat transfer coefficient computations:

Next, go to the Conditions panel and under Droplet reference conditions enter the following:

Under the Model parameters section for icing, change the Relative humidity to 100%.

Go to the Boundaries panel and set the following for each row:

Define the inlet conditions of the first stage (nacelle intake). In this case, apply the following boundary conditions for Row05:

For the Inlet boundary BC_1000, click Import reference conditions to assign values for Liquid Water Content, Temperature and Velocity. Repeat this step for BC_1001.

Go to the Solver panel. The Time integration – Droplets settings are already defined from the previous run. Ensure the CFL number is set to 10 and the Maximum number of time steps is set to 600. For the Icing simulation settings, under Time settings unselect the Automatic time step option and set the Time step to 1e-5 and the Total time of ice accretion to 10 seconds.

In the Advanced solver settings section set Mass deficit cutoff to 0.1%:

Click the Run button at the bottom of the panel to go to the Execution environment. In the Settings panel, set the total Number of CPUs to the maximum available. Click the Start menu button to start the calculation.

Figure 8.18: Comparison at the Fan Exit of LWC Before (Left) and After (Right) Accounting for Film Reinjection of Fan and Spinner Walls.

Create a new project and name it nosecone_CHT. Make sure that the metric units system has been selected for this project.

Create a FENSAP run in this project and name it FENSAP_IPS_On.

Double-click the grid icon and select the grid file nosecone_external.grd provided in the tutorials subdirectory ../workshop_input_files/Input_Grid/nosecone.

Figure 8.19: Nose Cone Grid

Double-click the config icon of FENSAP_IPS_On to proceed to the input parameters.

In the Model panel, keep the Navier-Stokes option for the Momentum equations and Full PDE for the Energy equation.

Keep the Spalart-Allmaras turbulence model and the Eddy/Laminar viscosity ratio at the low value of 1e-5.

The following calculation will be carried out in the rotational frame of reference. Select Rotational velocity for the Body forces and set a rotational speed of 1380 rpm along the X axis.

The conditions for this tutorial are set to a typical landing scenario, with low forward velocity and engine rotation rate, at sea level. Typical commercial jets approach at landing speeds in the range of 160 – 220 knots. Here you use 100 m/s which is about 190 knots. Set the Reference conditions as:

In the Initial solution section, select the Velocity components option and set Velocity X equal to 100 m/s while keeping the other two components to zero.

In order to reduce the amount of time required to run this CHT tutorial, a converged flow solution is provided. Change the Initial solution option from Velocity components to Solution restart, and select the solution file nosecone_flow_restart from tutorials subdirectory ../workshop_input_files/Input_Grid/nosecone. This step can be skipped if a full initial flow is desired.

Continue to the Domains tab. There are two domains listed in the panel. Domains in FENSAP-ICE are volume parts marked by different material IDs in the grid file. For more details on creating such grids, see FENSAP-ICE File Formats within the FENSAP-ICE User Manual.

To identify the material domains, right-click the grid icon in the project window and display it with Viewmerical. Click MAT_0 to highlight this domain, which happens to be the nose cone cavity.

Now that Domain_0 is identified as the internal cavity, set its initialization as Rotating and make sure that Rotor (unsteady) is set as Fixed. For more information on Rotating and Rotor (unsteady), refer to Ice Crystal Impingement and Ice Accretion. This will set the axial velocity to zero while imposing the rotational frame velocity on the internal nodes such that the air in the internal cavity is rotating with the same speed as the cone itself. This would eventually be the case for this geometry where the boundary layer will entrain the rest of the air. However, it would take quite a few iterations to establish this condition if initialized otherwise.

Continue to the Boundaries panel. Select the inlet boundary BC_1001 and choose Subsonic under Type. Click the Import reference conditions button to set 265 K for Temperature and 100 m/s for Velocity X.

Select the inlet boundary BC_1002 and choose Subsonic in the boundary Type section. Set the Temperature to 400 K and the Velocity X component to -300 m/s. Set the inlet Reference frame to Relative. Setting the reference frame for an Inlet applies the currently imposed velocity components in that frame. If Relative frame is chosen, the zero components in the Y and Z directions will apply in the rotating frame of reference. This means that in the absolute frame, the flow will assume a swirl component, following the rotation of the nose cone.

Select wall boundary BC_2001 and set its Heat flux to 0 (Adiabatic). Set the Rotation to Counter-rotating. Counter-rotation makes a wall, that is normally stationary in the absolute frame, rotating in the opposite direction in the relative frame. This wall is the shroud and it is stationary in the absolute frame.

Select each wall boundary: 2010, 2040, 2050 and set their Temperature to 320 K. These walls will compute non-zero heat fluxes which will then be used with ICE3D and C3D for heat flux exchange during CHT3D computations.

Select wall boundary 2020 and set its Heat Flux to 0. Icing and heat transfer are disabled on this wall.

Select the outflow boundary BC_3000 and choose Subsonic under Type. Click Import reference conditions to set the exit pressure value from the reference conditions.

Switch to the Solver panel. Set the CFL number to 50. If running with a restart file, set the Maximum number of time steps to 0 and uncheck the Use variable relaxation option. Otherwise set it to 1000 and turn on the Use variable relaxation option. This will increase the CFL number from 1 to 50 in 300 steps linearly.

In the Out tab, set solution output frequency to 20 iterations.

Run this calculation on as many CPUs as possible.

Figure 8.20: Nose Cone: Shear Stress (Left) and Classical Heat Flux (Right)

Figure 8.21: Nose Cone: Mach Number Distribution on a Cross-Section, Showing the Jet on the Inside and the Air Exhausting on the Outside

Create a new DROP3D run and name it DROP3D_IPS_On.

Drag & drop the config icon of FENSAP_IPS_On onto the config icon of the new DROP3D_IPS_On run. This copies the reference conditions of the flow to the droplet run.

Double-click the DROP3D_IPS_On config icon to edit the input parameters.

Check that the following Droplets reference conditions in the Conditions panel have been set by default:

To reduce the amount of time required to run this tutorial, a converged droplet solution is provided. In the Droplets initial solution section, change the Velocity components option to Solution restart. Select the solution file nosecone_droplet_restart.

In the Boundaries panel, select BC_1001 (Inlet) and click Import reference conditions. For BC_1002 (Inlet) set LWC to 0 g/m³ and check the droplet velocity components and set its values to 0 m/s. The jet does not inject water into the nose cone.

In the Solver panel, keep the default CFL number at 20. If running with a restart file, set the Maximum number of time steps to 1. Otherwise set the Maximum number of time steps to 500.

Run the calculation on as many CPUs as possible.

Figure 8.22: Droplet Collection Efficiency and LWC Distribution

Create a new ICE3D run and name it ICE3D_IPS_On.

Drag & drop the config icon of DROP3D_IPS_On onto the config icon of this new run. This operation automatically links the air and droplet solutions, the grid of the external domain and the reference conditions into the ICE3D_IPS_On run.

Double-click the config icon to edit the input parameters.

In the Model panel, set the Icing model to Glaze - Advanced and select Classical for the Heat flux type.

In the Conditions panel, set the Recovery factor to 0.9. In the Advanced section, activate the Hot chamber reference evaporation conditions. Set the following conditions:

ICE3D computes the surface temperature of dry regions as the recovery temperature, using the recovery factor. Since the inside and the outside of the nose cone are in the same computational domain and are subject to very different thermal conditions, it is necessary to provide ICE3D with an additional set of reference conditions to apply when assuming a recovery temperature for the inner walls. ICE3D will be able to automatically detect the wall nodes that are under the influence of the hot air inside, and apply the second set of reference conditions.

In the Boundaries panel, disable icing on wall boundaries 2001 and 2020. Enable the water sink option by setting the Sink boundary condition flag to 2050. This will avoid water pooling at the edge of the holes. The current grid resolution is not enough to capture the exact behavior of the water film in these areas, and it is better to declare them as sink.

In the Solver panel, keep the Automatic time step option enabled. Set the Total time of ice accretion to 30 seconds.

Go to the run window and start the execution on as many CPUs as possible. There is no need to create a displaced grid for this case.

Figure 8.23: Initial Film Coverage and the 30-Second Ice Thickness on the Nose Cone

Create a new CHT3D run and name it CHT3D_nosecone. The CHT configuration window will ask for the type of CHT simulation desired. Select Electrothermal (1 fluid, 1 solid) in the Problem type pull-down menu. Choose Wet air and the Icing mode to Anti-icing. Press the OK button to continue with the set-up. Although there is no electro-thermal device used in this simulation, you selected the Electro-Thermal (1 fluid 1 solid) option as this CHT3D simulation only uses one solid and one fluid domain. A tree of coupled FENSAP (fluid_ext), ICE3D (ice_ext) and C3D (solid) runs will appear in the run window.

Drag & drop the config icon of FENSAP_IPS_On onto the config icon of the fluid_ext run in CHT3D_nosecone.

Double-click the config icon of the fluid_ext run to edit the input parameters. Go to the Solver panel. Ensure that the Use variable relaxation option is unchecked. Increase the CFL number to 20,000. Change the Maximum number of time steps to 50. CHT3D will perform 50 FENSAP iterations when updating the flow at every CHT loop. Change Convergence level from 1e-10 to 1e-12. This will ensure that the heat fluxes are properly converged on the flow domain at each CHT time step. Close and save the configuration file.

Drag & drop the config icon from ICE3D_IPS_On onto the config icon of the ice_ext of the CHT3D_nosecone run. Double-click the config icon of the ice_ext run to edit the input parameters. Go to the Solver panel and set the Total ice accretion time to 10 seconds. Close and save the configuration file. In each CHT3D iteration, ICE3D will only run for 10 seconds. This is a long enough time for the film to settle between each CHT loop.

To set up the solid run in CHT3D_nosecone, double-click the grid icon and select the file nosecone_solid.grid provided in the tutorials subdirectory ../workshop_input_files/Input_Grid/nosecone.

Double-click the config icon of the solid run to edit the input parameters for heat conduction. Go to the Settings panel. Set the Temperature value to 273.15 K. This will be the initial temperature throughout the solid domain.

Go to the Properties panel. Rename the default material to Duralumin. Specify the following characteristics for the material.

Go to the Boundaries panel. For BC_2060, select the Flux boundary condition and set the Heat flux to 0 W/m². This wall does not share an interface with the grid of the airflow and therefore it is considered as an adiabatic wall.

For boundary walls 2010, 2040, 2050, select the Nothing boundary condition. These walls share an interface with the grid of the airflow in the main CHT configuration panel.

Go to the Solver panel. Set the Maximum time step value to 0.2 seconds and the Total time to 5 seconds. At each CHT3D iteration, heat conduction will be computed for 5 seconds.

Close and save this configuration.

Double-click the main config icon of CHT3D_nosecone to set up the links between the solid and fluid domains.

Go to the Parameters panel. Set the Number of CHT iterations (loops) to 20 and the Solution output every 5 iterations. The simulation will perform 20 CHT loops and save solutions at every 5 loops.

Choose the Solve energy only option in the Flow solver mode - external pull-down menu. The flow dynamics will not change much when the surface temperatures change due to conduction. Simulating CHT3D cases using only the energy equation generally allows us to obtain a solution of sufficient accuracy while also saving significant computational time since only one out of six equations is solved.

Go to the Interfaces panel. There are 3 solid-fluid interfaces to assign. Set up the interfaces as follows:

Close and save the configuration.

Right-click the main config icon of CHT3D_nosecone, then select the Run option in the menu to launch the CHT3D calculation on as many CPUs as possible. The convergence graphs of minimum and maximum temperatures in the solid, as well as the airflow residuals are shown below.

Figure 8.24: Solid Minimum (Left) and Maximum (Right) Temperature

Figure 8.25: Average Residual of the Flow

The converged solid temperature distribution is shown below. The hottest region is the tip where the hot jet is concentrated and the water catch is low since the tip is in a shadow zone. Towards the back of the cone, the temperatures are cooler due to convective and evaporative cooling. The minimum temperature on the surface is above freezing, meaning that there is no ice formation at steady state IPS operation. This can be verified by loading the ICE3D solution and looking at the Instant. Ice Growth (kg/m^2s) solution field. There is no need to run ICE3D post CHT for this simulation since the ice accretion rate is zero.

Figure 8.26: Solid Temperature of the Nose Cone

Figure 8.27: Film Height (Left) and Ice Growth Rate (Right) on the Surface with Anti-Icing

Create a new project and name it nosecone_CHT_FLUENT. Make sure that the metric units system has been selected for this project. Close FENSAP-ICE.

Go to the project folder, nosecone_CHT_FLUENT, and create a new sub folder called INITIAL-AIR. Copy the provided Fluent case file, nosecone_external-FLUENT.cas.h5, from the tutorials subdirectory ../workshop_input_files/Input_Grid/nosecone to this new folder and launch Fluent.

Figure 8.32: Nose Cone Grid

In the Fluent Launcher window, select the Dimension as 3D, pick Double Precision under Options, and choose Parallel (Local Machine) under Processing Options. Assign a number of CPUs for this simulation, 2 to 4 CPUs, under Solver → Processes. Click Show More Options. Under General Options, set your Working Directory to the INITIAL-AIR directory. Press OK to close the Fluent Launcher.

Read the case file by going to the File → Read → Case menu and browse to and select the Fluent case file, nosecone_external-FLUENT.cas.h5, located inside the sub folder INITIAL-AIR.

From the top navigation menu, select Physics. Make sure the Solver is set to Pressure-Based, Absolute, and Steady. Click the Operating Conditions. Set the Operating Pressure to 0 Pa. Press OK to close the Operation Conditions window.

Expand the Setup → Material → Fluid from the side tree menu. Double-click air and modify its properties.

The table below describes the air properties to be imposed in this simulation.

Click the Change/Create button and Close this window to save the new air properties.

---

Note:  Thermal conductivity and viscosity are computed from these equations which are presented in the FENSAP-ICE User Manual and shown below:

In these equations, T_∞ refers to the ambient air static temperature, and C₁, T_ref and µ_ref are respectively equal to 0.00216176 W/m/K^3/2, 288 K and 17.9*10^-6 Pa s.

---

Expand and double-click Setup → Models → Energy from the side tree menu. Ensure that Energy is turned on.

Expand and double-click Setup → Models → Viscous from the side tree menu. Select the k-omega (2 eqn) option in the opened Viscous Model interface and then select k-omega Model → SST. Make sure that Viscous Heating and Production Limiter are activated in the Options section. In the Model Constants section, drag the scroll bar down and set Energy Prandtl Number and Wall Prandtl Number values to 0.9. Click OK to close this menu.

Expand and double-click Setup → Cell Zone Conditions from the side tree menu.

In the Task Page → Zone list, double-click live-external to open the Fluid window. Check the Frame Motion. Under Reference Frame panel, make sure absolute is selected in the drop-box of Relative To Cell Zone under Relative Specification. Set Rotation-Axis Origin coordinates (X, Y, Z) to (0, 0, 0) meter and Rotation-Axis Direction vector (X, Y, Z) to (1, 0, 0). Enter 144.44 rad/s as the Speed of Rotational Velocity, which corresponds to 1380 rpm. Click OK to close this window. Repeat this step for the live-internal zone.

Expand and double-click Setup → Boundary Conditions from the side tree men. In the Task Page, click Periodic Conditions. Select Specify Pressure Gradient and set the Flow Direction vector (X, Y, Z) to (1, 0, 0). Then close this window.

In the Task Page → Boundary Conditions list, click the pressure-far-field-1 and then select pressure-far-field from the drop-down list of Type. Click Edit to open the Pressure Far-Field window and set the following conditions for this boundary:

In the Task Page → Boundary Conditions list, click the pressure-inlet-5 and then select pressure-inlet from the drop-down list of Type. Click Edit to open the Pressure Inlet window and set the following conditions for this boundary:

In the Task Page → Boundary Conditions list, click the pressure-outlet-11 and then select pressure-outlet from the drop-down list of Type. Click Edit to open the Pressure Outlet window and set the following conditions for this boundary:

In the Task Page → Boundary Conditions list, double-click wall-7 to open the Wall window. In the Momentum panel, set the Wall Motion to Moving Wall and the Shear Condition to No Slip. Set the Motion to Absolute and Rotational; set the Speed to 144.44 rad/s, the Rotation-Axis Origin coordinates (X, Y, Z) to (0, 0, 0) meter and the Rotation-Axis Direction vector (X, Y, Z) to (1, 0, 0). In the Thermal panel, set the Thermal Conditions to a Temperature of 320 K. Click OK to close the window. Repeat this step for wall-8, wall-9, wall-10, and wall-17.

In the Task Page → Boundary Conditions list, double-click wall-6 to open the Wall window. In the Momentum panel, set the Wall Motion to Moving Wall and the Shear Condition to No Slip. Set the Motion to Absolute and Rotational; set the Speed to 0 rad/s, the Rotation-Axis Origin coordinates (X, Y, Z) to (0, 0, 0) meter and the Rotation-Axis Direction vector (X, Y, Z) to (1, 0, 0). In the Thermal panel, set the Thermal Conditions to a Heat Flux of 0 w/m^². Click OK to close the window.

Expand and double-click Setup → Reference Values from the side tree menu. In the Task Page, use the Compute from drop-down menu to select pressure-far-field-1.

Expand and double-click Solution → Methods from the side tree menu. Set the Pressure-Velocity Coupling section to Coupled. In the Spatial Discretization section, set the Gradient to Green-Gauss Node Based and the remaining settings to Second Order or Second Order Upwind. Make sure that the Pseudo Transient and High Order Term Relaxation options are disabled.

Expand and double-click Solution → Controls from the side tree menu. Under Solution → Controls, set the Density and Body Forces of Pseudo Transient Explicit Relaxation Factors to 0.5.

Expand and double-click Solution → Monitors → Residual from the side tree menu. In the Residual Monitors window, make sure that the Print to Console and the Plot are enabled. Disable all check box below the Convergence column. Close this window once this is done.

Expand and double-click Solution → Initialization from the side tree menu. In the opened Task Page, select Hybrid Initialization under Initialization Methods. Click the Initialize button to initialize the computational domain.

Expand and double-click Solution → Run Calculation from the side tree menu. Under Pseudo Transient Options, set the Time Step Method to Automatic and the Timescale Factor to 0.1. Enter 10000 as the total Number of Iterations and click Calculate to start the simulation.

Once the simulation is complete, save the new Fluent solution in File → Write → Case & Data. Name this simulation as nosecone_external-FLUENT.cas.h5/dat.h5.

The following figures compare the Fluent solution to a FENSAP kw-sst airflow solution.