User-Defined Wind Simulation

Running a steady-flow wind simulation on a closed box building with one load case defined per wind direction: Activate the wind simulation add-on and orient the global Z-axis upwards Model a closed box building of four walls and a flat roof on a fixed base Define a horizontal terrain plane as the floor of the CFD wind tunnel Create the shrink wrapping, surface roughness and a user-defined wind profile Define steady-flow analysis settings and one load case per wind direction Calculate the CFD simulation and read the wind force resultants, support forces and deformations Keywords: wind simulation CFD wind profile terrain wind tunnel wind load case

PythonC#

from math import inf, radians

from dlubal.api import rfem, common

# -------------------------------------------------------
# This example runs the wind simulation add-on on a simple closed box building
# standing on a horizontal terrain plane. RFEM hands the model over to the CFD
# solver, once per wind direction, and the resulting surface pressures are
# applied to the FE model and solved by the following static analysis.
#
# Every wind direction is a separate CFD run. With the settings below one
# direction takes roughly 3-4 minutes on a regular workstation, so the whole
# script runs for several minutes.
# -------------------------------------------------------

# Building dimensions [m]
LENGTH = 10.0     # along global X
WIDTH = 8.0       # along global Y
HEIGHT = 6.0      # along global Z
THICKNESS = 0.2   # wall and roof thickness

# Wind directions around the global Z-axis, clockwise, one load case each [deg]
WIND_DIRECTIONS = (0.0, 90.0)


def define_structure() -> list:
    """Define a closed box building - four walls and a flat roof."""

    return [

        # Material
        rfem.structure_core.Material(
            no=1,
            name='C30/37 | EN 1992-1-1:2004/A1:2014',
        ),

        # Nodes - base
        rfem.structure_core.Node(no=1, coordinate_1=0.0, coordinate_2=0.0, coordinate_3=0.0),
        rfem.structure_core.Node(no=2, coordinate_1=LENGTH, coordinate_2=0.0, coordinate_3=0.0),
        rfem.structure_core.Node(no=3, coordinate_1=LENGTH, coordinate_2=WIDTH, coordinate_3=0.0),
        rfem.structure_core.Node(no=4, coordinate_1=0.0, coordinate_2=WIDTH, coordinate_3=0.0),

        # Nodes - top
        rfem.structure_core.Node(no=5, coordinate_1=0.0, coordinate_2=0.0, coordinate_3=HEIGHT),
        rfem.structure_core.Node(no=6, coordinate_1=LENGTH, coordinate_2=0.0, coordinate_3=HEIGHT),
        rfem.structure_core.Node(no=7, coordinate_1=LENGTH, coordinate_2=WIDTH, coordinate_3=HEIGHT),
        rfem.structure_core.Node(no=8, coordinate_1=0.0, coordinate_2=WIDTH, coordinate_3=HEIGHT),

        # Lines - base
        rfem.structure_core.Line(no=1, definition_nodes=[1, 2]),
        rfem.structure_core.Line(no=2, definition_nodes=[2, 3]),
        rfem.structure_core.Line(no=3, definition_nodes=[3, 4]),
        rfem.structure_core.Line(no=4, definition_nodes=[4, 1]),

        # Lines - vertical
        rfem.structure_core.Line(no=5, definition_nodes=[1, 5]),
        rfem.structure_core.Line(no=6, definition_nodes=[2, 6]),
        rfem.structure_core.Line(no=7, definition_nodes=[3, 7]),
        rfem.structure_core.Line(no=8, definition_nodes=[4, 8]),

        # Lines - top
        rfem.structure_core.Line(no=9, definition_nodes=[5, 6]),
        rfem.structure_core.Line(no=10, definition_nodes=[6, 7]),
        rfem.structure_core.Line(no=11, definition_nodes=[7, 8]),
        rfem.structure_core.Line(no=12, definition_nodes=[8, 5]),

        # Surfaces - four walls and the roof
        rfem.structure_core.Surface(no=1, boundary_lines=[1, 6, 9, 5]),
        rfem.structure_core.Surface(no=2, boundary_lines=[2, 7, 10, 6]),
        rfem.structure_core.Surface(no=3, boundary_lines=[3, 8, 11, 7]),
        rfem.structure_core.Surface(no=4, boundary_lines=[4, 5, 12, 8]),
        rfem.structure_core.Surface(no=5, boundary_lines=[9, 10, 11, 12]),

        # Thickness
        rfem.structure_core.Thickness(
            no=1,
            material=1,
            uniform_thickness=THICKNESS,
            assigned_to_surfaces=[1, 2, 3, 4, 5],
        ),

        # Fixed line support along the base
        rfem.types_for_lines.LineSupport(
            no=1,
            lines=[1, 2, 3, 4],
            spring=common.Vector3d(x=inf, y=inf, z=inf),
        ),

        # Static analysis settings of the wind load cases
        rfem.loading.StaticAnalysisSettings(
            no=1,
            analysis_type=rfem.loading.StaticAnalysisSettings.ANALYSIS_TYPE_GEOMETRICALLY_LINEAR,
        ),
    ]


def define_terrain() -> rfem.base_data_objects.Terrain:
    """Define a horizontal terrain plane - the floor of the CFD wind tunnel."""

    return rfem.base_data_objects.Terrain(
        no=1,
        type=rfem.base_data_objects.Terrain.TYPE_HORIZONTAL_PLANE,
        center_of_terrain_z=0.0,
        bounding_box_offset_x=30.0,
        bounding_box_offset_y=30.0,
    )


def define_wind_simulation() -> list:
    """Define the wind simulation objects.

    The CFD model is built from a shrink-wrapped envelope of the structure, so
    a shrink wrapping is required for the main structure and one for the
    surrounding objects. Both are referenced by the mesh settings further down.
    The wind profile describes the atmospheric boundary layer at the tunnel
    inlet, and one load case is created per wind direction.
    """

    return [

        # Shrink wrapping - main structure
        rfem.types_for_wind_simulation.ShrinkWrapping(
            no=1,
            simplification_defined_by=rfem.types_for_wind_simulation.ShrinkWrapping.SIMPLIFICATION_DEFINED_BY_LEVEL_OF_DETAILS,
            level_of_detail=2,
            small_openings_closure_type=rfem.types_for_wind_simulation.ShrinkWrapping.SMALL_OPENINGS_CLOSURE_TYPE_PERCENT_OF_MODEL_DIAMETER,
            closure_relative_to_model_parameter=0.05,
            deactivate_shrink_wrapping=False,
        ),

        # Shrink wrapping - surrounding objects
        rfem.types_for_wind_simulation.ShrinkWrapping(
            no=2,
            simplification_defined_by=rfem.types_for_wind_simulation.ShrinkWrapping.SIMPLIFICATION_DEFINED_BY_LEVEL_OF_DETAILS,
            level_of_detail=0,
            small_openings_closure_type=rfem.types_for_wind_simulation.ShrinkWrapping.SMALL_OPENINGS_CLOSURE_TYPE_PERCENT_OF_MODEL_DIAMETER,
            closure_relative_to_model_parameter=0.2,
            deactivate_shrink_wrapping=False,
        ),

        # Surface roughness of the model in the wind tunnel
        rfem.types_for_wind_simulation.RoughnessAndPermeability(
            no=1,
            type_of_surface=rfem.types_for_wind_simulation.RoughnessAndPermeability.TYPE_OF_SURFACE_SMOOTH,
        ),

        # Wind profile at the tunnel inlet
        rfem.load_wizards.WindProfile(
            no=1,
            type=rfem.load_wizards.WindProfile.TYPE_USER_DEFINED,
            user_defined_input_type=rfem.load_wizards.WindProfile.USER_DEFINED_INPUT_TYPE_CONSTANT,
            user_defined_inlet_variables=rfem.load_wizards.WindProfile.USER_DEFINED_INLET_VARIABLES_I_TUL,
            user_defined_wind_profile=rfem.load_wizards.WindProfile.UserDefinedWindProfileTable(
                rows=[
                    # height [m], velocity [m/s], turbulence intensity [-]
                    rfem.load_wizards.WindProfile.UserDefinedWindProfileRow(
                        no=1, height=0.0, velocity=15.0, turbulence_intensity=0.20),
                    rfem.load_wizards.WindProfile.UserDefinedWindProfileRow(
                        no=2, height=5.0, velocity=20.0, turbulence_intensity=0.18),
                    rfem.load_wizards.WindProfile.UserDefinedWindProfileRow(
                        no=3, height=10.0, velocity=23.0, turbulence_intensity=0.16),
                    rfem.load_wizards.WindProfile.UserDefinedWindProfileRow(
                        no=4, height=20.0, velocity=26.0, turbulence_intensity=0.14),
                    rfem.load_wizards.WindProfile.UserDefinedWindProfileRow(
                        no=5, height=40.0, velocity=29.0, turbulence_intensity=0.12),
                ]
            ),
        ),

        # Steady-flow CFD simulation with a reduced iteration count
        rfem.loading.WindSimulationAnalysisSettings(
            no=1,
            simulation_type=rfem.loading.WindSimulationAnalysisSettings.SIMULATION_TYPE_STEADY_FLOW,
            turbulence_model_type=rfem.loading.WindSimulationAnalysisSettings.TURBULENCE_MODEL_TYPE_EPSILON,
            finite_volume_mesh_density=0.2,
            minimum_number_of_iterations=60,
            maximum_number_of_iterations=300,
            residual_type=rfem.loading.WindSimulationAnalysisSettings.RESIDUAL_TYPE_PRESSURE,
            residual_pressure=0.01,
            consider_turbulence=True,
            snap_to_model_edges=True,
            use_potential_flow_solver_for_initial_condition=True,
            member_load_distribution=rfem.loading.WindSimulationAnalysisSettings.MEMBER_LOAD_DISTRIBUTION_CONCENTRATED,
        ),

        # One load case per wind direction - the wind tunnel dimensions are
        # derived automatically from the bounding box of the model
        *[
            rfem.loading.LoadCase(
                no=index,
                name=f'Wind simulation | {direction:.0f} deg',
                analysis_type=rfem.loading.LoadCase.ANALYSIS_TYPE_WIND_SIMULATION,
                action_category=rfem.loading.LoadCase.ACTION_CATEGORY_WIND_QW,
                static_analysis_settings=1,
                wind_simulation_wind_profile=1,
                wind_simulation_analysis_settings=1,
                wind_simulation_wind_direction_angle=radians(direction),
                to_solve=True,
            )
            for index, direction in enumerate(WIND_DIRECTIONS, start=1)
        ],
    ]


with rfem.Application() as rfem_app:

    # Create an empty model
    rfem_app.close_all_models(save_changes=False)
    rfem_app.create_model(name='wind_simulation_user_defined')

    # Activate the wind simulation add-on and orient the global Z-axis upwards
    base_data: rfem.BaseData = rfem_app.get_base_data()
    base_data.addons.wind_simulation_active = True
    base_data.general_settings.global_axes_orientation = \
        rfem.BaseData.GeneralSettings.GLOBAL_AXES_ORIENTATION_ZUP
    rfem_app.set_base_data(base_data=base_data)

    # Clean up the objects pre-created by the wind simulation add-on
    rfem_app.delete_all_objects()

    # Create the structure and the wind simulation objects
    rfem_app.create_object_list(
        objs=define_structure() + define_wind_simulation()
    )

    # Terrain no. 1 always exists in the model and is not removed by
    # 'delete_all_objects', so it is updated instead of created
    rfem_app.update_object(
        obj=define_terrain()
    )

    # Mesh settings of the wind simulation - shrink wrapping, terrain, solver run
    mesh_settings = rfem_app.get_mesh_settings()
    mesh_settings.wind_simulation.consider_terrain_enabled = True
    mesh_settings.wind_simulation.shrink_wrapping_main_structure = 1
    mesh_settings.wind_simulation.shrink_wrapping_surrounding_objects = 2
    mesh_settings.wind_simulation.consider_surface_thickness_above_enabled = True
    mesh_settings.wind_simulation.consider_surface_thickness_above_value = 0.005
    mesh_settings.wind_simulation.run_rwind_in_background_enabled = True
    rfem_app.set_mesh_settings(mesh_settings=mesh_settings)

    for index, direction in enumerate(WIND_DIRECTIONS, start=1):
        load_case = rfem_app.get_object(rfem.loading.LoadCase(no=index))
        print(f'LC{index} | {direction:.0f} deg | wind tunnel '
              f'{load_case.wind_simulation_wind_tunnel_depth:.1f} x '
              f'{load_case.wind_simulation_wind_tunnel_width:.1f} x '
              f'{load_case.wind_simulation_wind_tunnel_height:.1f} m')

    # Run the CFD simulation of every wind load case and the static analysis
    print('\nCalculating - this takes a few minutes per wind direction ...')
    print(rfem_app.calculate_all(skip_warnings=True).message)

    # Retrieve results per wind direction. The result tables are wider than the
    # terminal, so they are printed with 'to_string()' to keep them untruncated.
    for index, direction in enumerate(WIND_DIRECTIONS, start=1):

        loading = rfem.ObjectId(no=index, object_type=rfem.ObjectType.OBJECT_TYPE_LOAD_CASE)
        print(f'\n=== LC{index} | wind direction {direction:.0f} deg ===')

        # Solver messages, e.g. on the convergence of the residual pressure
        errors_and_warnings = rfem_app.get_result_table(
            table=rfem.results.ResultTable.ERRORS_AND_WARNINGS_TABLE,
            loading=loading,
        )
        print(f'\nErrors and Warnings:\n{errors_and_warnings.data.to_string(index=False)}')

        # Summary of the static analysis - the wind force resultants handed over
        # by the CFD solver, checked against the support forces. The 'value'
        # column holds SI base units (N, Nm, m).
        summary = rfem_app.get_result_table(
            table=rfem.results.ResultTable.STATIC_ANALYSIS_SUMMARY_TABLE,
            loading=loading,
        )
        totals = summary.data[summary.data['description'].str.contains(
            'sum of loads|sum of support forces', na=False)]
        print(f'\nWind Force Resultants and Support Forces [N]:\n'
              f'{totals[["description", "value"]].to_string(index=False)}')

    # Support forces of all wind load cases along the base lines [N/m, Nm/m]
    support_forces = rfem_app.get_results(
        results_type=rfem.results.ResultsType.STATIC_ANALYSIS_LINES_SUPPORT_FORCES
    )
    print(f'\nLines Support Forces:\n{support_forces.data.head(10).to_string(index=False)}')

    # Deformations of the walls and the roof in the mesh nodes [m, rad]
    deformations = rfem_app.get_results(
        results_type=rfem.results.ResultsType.STATIC_ANALYSIS_SURFACES_GLOBAL_DEFORMATIONS_MESH_NODES
    )
    print(f'\nSurfaces Global Deformations:\n{deformations.data.head(10).to_string(index=False)}')

    # Governing deformation per wind direction
    governing = deformations.data.loc[
        deformations.data.groupby('loading')['u_abs'].idxmax(),
        ['loading', 'surface_no', 'mesh_node_no', 'u_abs', 'u_x', 'u_y', 'u_z'],
    ]
    print(f'\nMaximum Absolute Deformation per Load Case [m]:\n{governing.to_string(index=False)}')