Optimization example using PyMAPDL and SciPy Optimize library

Vishnu
Vishnu Member, Employee Posts: 222
100 Comments 100 Likes Second Anniversary Name Dropper
✭✭✭✭
edited February 7 in Structures

Below is a PyMAPDL example using FLUID218 elements (3D Hydrodynamic Bearing Element) on how you could find an optimized Inputs (Shaft center X and Y location) for a preferred FLUID forces in X and Y directions.

Please refer FLUID218 Documentation below
https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/v241/en/ans_elem/Hlp_E_FLUID218.html?q=FLUID218

Tagged:

Comments

  • Vishnu
    Vishnu Member, Employee Posts: 222
    100 Comments 100 Likes Second Anniversary Name Dropper
    ✭✭✭✭
    edited February 13

    The first step in optimization is formulating an objective function, in this case the function named extract_fluid_forces is evaluated by using PyMAPDL

    scipy.optimize ,minimize is used for optimization .please refer the documentation for it below:

    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

    Example code below

    # Import necessary libraries
    import matplotlib.pyplot as plt
    import numpy as np
    from ansys.dpf import core as dpf
    from scipy.optimize import minimize
    
    from ansys.mapdl.core import launch_mapdl
    
    # Define the path to the result file
    rstpath = r"D:\apdl\test\file.rst"
    
    # Launch MAPDL with specified run location and number of processors
    mapdl = launch_mapdl(run_location=r"D:\apdl\test",nproc=1,override=True)
    
    mapdl.clear("All")
    mapdl.prep7()
    mapdl.et('', 'FLUID218')
    mapdl.mp(lab='VISC', mat='1', c0='0.135')
    mapdl.mp(lab='DENS', mat='1', c0='910')
    mapdl.r(nset='1', r1='0', r2='10')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='0', rad2='', theta2='90', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='90', rad2='', theta2='180', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='180', rad2='', theta2='270', depth='50')
    mapdl.cyl4(xcenter='', ycenter='', rad1='10', theta1='270', rad2='', theta2='360', depth='50')
    mapdl.vdele(nv1='ALL', nv2='', ninc='', kswp='0')
    mapdl.run('CSYS,1')
    mapdl.asel('U', 'LOC', 'X', 9, 10)
    mapdl.adele('ALL')
    mapdl.allsel()
    mapdl.mshape(0, '2D')
    mapdl.mshkey(1)
    mapdl.lesize('ALL', size='', angsiz='',  ndiv=20)
    mapdl.amesh('ALL')
    mapdl.eplot()
    mapdl.nummrg('ALL')
    clearance = 0.1
    mapdl.rmodif(1, 8, clearance)
    mapdl.rmodif(1, 9, clearance)
    mapdl.rmodif(1, 10, clearance)
    mapdl.rmodif(1, 11, clearance)
    mapdl.nsel('S', 'LOC', 'Z', '', 0)
    mapdl.d('ALL', 'PRES', 0)
    mapdl.allsel()
    mapdl.nsel('S', 'LOC', 'Z', 50, 50)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.nsel('S', 'LOC', 'Y', -20, 20)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.nsel('S', 'LOC', 'Y', 160, 200)
    mapdl.d('ALL', 'PRES', 100)
    mapdl.allsel()
    mapdl.omega('', '', 10)
    mapdl.finish()
    
    # Define a function to extract fluid forces
    def extract_fluid_forces(xoffset,yoffset):
            mapdl.prep7()
            mapdl.rmodif(1, 3,xoffset)
            mapdl.rmodif(1, 4,yoffset)
            mapdl.solution()
            mapdl.solve()       
            # data source
            my_data_sources = dpf.DataSources(rstpath)
    
            # model
            my_model = dpf.Model(rstpath)
    
            # mesh
            my_mesh = my_model.metadata.meshed_region
    
            # Field Containers at all mesh
            # Get result operator
            result_operator = dpf.operators.result.__getattribute__("nmisc")
    
            result_field_container = result_operator(
                time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh,item_index=11
            ).outputs.fields_container.get_data()
    
            # First Field usually field corresponding to first time scoping
            resultfield = result_field_container[0]
    
            FXTOTAL=sum(resultfield.data_as_list)
    
            result_field_container = result_operator(
                time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh,item_index=12
            ).outputs.fields_container.get_data()
    
            # First Field usually field corresponding to first time scoping
            resultfield = result_field_container[0]
    
            FYTOTAL=sum(resultfield.data_as_list)
    
            fluid_forces = np.array([FXTOTAL,FYTOTAL])
            #plt.bar(range(len(fluid_forces)), fluid_forces)
    
            return fluid_forces
    
    
    # preferred forces
    preferred_fx = 10000  # replace with your preferred FXtotal
    preferred_fy = 90000  # replace with your preferred FYtotal
    
    # Define objective function
    def objective(offsets):
        try:
            [fx, fy] = extract_fluid_forces(offsets[0],offsets[1])
            return (((fx - preferred_fx)**2 + (fy - preferred_fy)**2))
        except:  # catch all exceptions
            return 9E50  # return a large value
    
    
    # initial guess for x and y offsets
    x0 = [0.0,0.0]
    
    # Minimize the objective function
    result = minimize(objective, x0,method='nelder-mead',options={'maxfev':1000,'maxiter':1000,'disp': True,'fatol' : 1e-10, 'xatol' : 1e-4})
    
    # Get the optimized x and y offsets
    x, y = result.x
    fx,fy = extract_fluid_forces(x, y)
    
    # Exit MAPDL
    mapdl.exit()
    
    # Print the optimized values
    print("Optimized X and Y Offsets are %f and %f"%(x,y))
    print("The Fx,Fy forces at these X and Y offsets are %f and %f "%(fx,fy))
    
    
    # Define the labels for the bars
    labels = ['FX', 'FY']
    
    # Create a figure and a set of subplots
    fig, ax = plt.subplots()
    
    # Define the width of the bars
    width = 0.35
    
    # Define the positions for the groups of bars
    x = np.arange(len(labels))
    
    # Plot the preferred and obtained values as bars next to each other
    ax.bar(x - width/2, [preferred_fx,preferred_fy], width, label='Preferred')
    ax.bar(x + width/2, [fx,fy], width, label='optimized')
    
    # Add some text for labels, title and custom x-axis tick labels, etc.
    ax.set_ylabel('Values')
    ax.set_title('Preferred and Obtained Values of FX and FY')
    ax.set_xticks(x)
    ax.set_xticklabels(labels)
    ax.legend()
    
    
    # Display the plot
    plt.show()
    
    
  • Vishnu
    Vishnu Member, Employee Posts: 222
    100 Comments 100 Likes Second Anniversary Name Dropper
    ✭✭✭✭
    edited February 6

    Output --->>>