Abstract:
Constraint-based methods for determining orientations of material physical properties using an isoparametric shape function are disclosed. In one embodiment, a method of defining an orientation of an material physical property includes defining nonlinear and/or discontinuous design constraints of design values in a geometric domain associated with one or more physical attributes of the material physical property, and translating the nonlinear and/or discontinuous design constraints into continuous, first order design constraints of the design values by applying an isoparametric shape function. The method further includes performing a topology optimization using the continuous, first order design constraints of the design values, and reverse-translating results of the topology optimization back into the geometric domain using the isoparametric shape function. The results of the topology optimization in the geometric domain are indicative of the orientation of the material physical property.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    The present application is a continuation of U.S. patent Ser. No. 13/747,955, Jan. 23, 2013. 
     
    
     TECHNICAL FIELD 
       [0002]    The present specification generally relates to methods for determining orientations of material physical properties and, more particularly, constraint-based methods for determining orientations of material physical properties using an isoparametric shape function. 
       BACKGROUND 
       [0003]    Computer aided engineering analysis, for example, topology optimization techniques such as finite element models (FEM), incorporate computational techniques that are used to find approximate solutions to engineering problems. In general, a model representing the geometry of interest is discretized into a plurality of elements that collectively represent the entire geometry. Because of the elements&#39; reduced size, computers executing the FEM analysis are capable of solving the partial differential equations that govern the physical properties of each of the elements and, therefore, provide an approximate solution to the physical problem for the entire geometry of interest. 
         [0004]    Defining in the FEM model the orientation of the material physical property or properties (e.g., magnetic north of a permanent magnet) in polar coordinates (i.e., radius and angle of orientation) imposes a non-continuous design constraint on the material physical property because the angle of orientation is between 0° and 360°. Similarly, limiting the orientation of the material physical property using a Cartesian coordinate system imposes a nonlinear design constraint on the property because the position of the property must satisfy a second-order condition (i.e., (x 2 +y 2 ) 1/2 &lt;1). 
         [0005]    Another difficulty in optimizing the orientation of material physical properties includes application of design constraints that are non-continuous. For example, properties may be limited in orientation to a pre-determined number of directions in the physical domain. For example, manufacturing constraints may limit the orientation of the property, such as the direction of magnetic fields of magnets in a magnetic system, into four, six, or eight directions. 
         [0006]    Using traditional optimization techniques, the material properties of each of the finite elements would have non-continuous design constraints applied, which increases optimization model complexity and computational time to complete convergence of the optimization model. Further, incorporating such a model with non-continuous design constraints into an optimization routine may be difficult or prohibitive, as the optimization model having non-continuous design constraints may lack the flexibility to make incremental changes to the orientation of the material physical properties of the elements while allowing the optimization model to converge. 
         [0007]    Accordingly, a need exists for alternative methods for convergence of nonlinear programming problems that involve angles as material property design variables, and promoting convergence of the material property design variable angles into a set of specified desired angles. 
       SUMMARY 
       [0008]    In one embodiment, a method of defining an orientation of a material physical property includes defining nonlinear and/or discontinuous design constraints of design values in a geometric domain associated with one or more physical attributes of the material physical property, and translating the nonlinear and/or discontinuous design constraints into continuous, first order design constraints of the design values by applying an isoparametric shape function. The method further includes performing a topology optimization using the continuous, first order design constraints of the design values, and reverse-translating results of the topology optimization back into the geometric domain using the isoparametric shape function. The results of the topology optimization in the geometric domain are indicative of the orientation of the material physical property. 
         [0009]    In another embodiment, a method of determining an arrangement of a plurality of material physical properties of a device includes providing initial design values ξ i  and η i  for each material physical property on a normalized coordinate system, wherein −1&lt;ξ&lt;1 and −1&lt;η&lt;1 for i=1, 2, . . . , n. The method further includes projecting the initial design values (ξ i , η i ) to projected initial design values (x i , y i ) on a Cartesian coordinate system using an isoparametric shape function Nx(ξ, η), Ny(ξ, η) having a plurality of nodes such that x i =Nx(ξ i , η i ), y i =Ny(ξ i , η i ), and √(x i   2 +y i   2 )&lt;1. The method further includes obtaining convergence in an objective function by execution of a topology optimization of x i , y i  and determining an orientation of each material physical property based on final design variables x i , y i  for each material physical property i of the plurality of material physical properties. 
         [0010]    These and additional features provided by the embodiments described herein will be more fully understood in view of the following detailed description, in conjunction with the drawings. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0011]    The embodiments set forth in the drawings are illustrative and exemplary in nature and not intended to limit the subject matter defined by the claims. The following detailed description of the illustrative embodiments can be understood when read in conjunction with the following drawings, where like structure is indicated with like reference numerals and in which: 
           [0012]      FIG. 1  schematically depicts an exemplary magnetic fluid cooling device according to one or more embodiments described and illustrated herein; 
           [0013]      FIG. 2  graphically depicts simulated magnetic field patterns of the magnetic fluid cooling device depicted in  FIG. 1  according to one or more embodiments described and illustrated herein; 
           [0014]      FIG. 3  graphically depicts a simulated temperature profile of the magnetic fluid cooling device depicted in  FIG. 1  and the magnetic field patterns depicted in  FIG. 2  according to one or more embodiments described and illustrated herein; 
           [0015]      FIG. 4  graphically depicts simulated magnetic fluid velocity magnitude contours with the magnetic fluid cooling device depicted in  FIG. 1  and the magnetic field patterns depicted in  FIG. 2  according to one or more embodiments described and illustrated herein; 
           [0016]      FIG. 5  schematically depicts an arrangement of magnet elements having continuous angles to produce the magnetic field patterns depicted in  FIG. 2  according to one or more embodiments described and illustrated herein; 
           [0017]      FIG. 6A  graphically depicts coordinates in polar form; 
           [0018]      FIG. 6B  graphically depicts coordinates in Cartesian form; 
           [0019]      FIG. 7A  graphically depicts a normalized coordinate system with side constraints according to one or more embodiments described and illustrated herein; 
           [0020]      FIG. 7B  graphically depicts an isoparametric shape function used to map ξ, η to x, y according to one or more embodiments described and illustrated herein; 
           [0021]      FIG. 8A  graphically depicts squares drawn on a normalized ξ, η coordinate system according to one or more embodiments described and illustrated herein; 
           [0022]      FIG. 8B  graphically depicts a four quadratic eight node isoparametric shape function according to one or more embodiments described and illustrated herein; 
           [0023]      FIG. 9  graphically depicts the isoparametric shape function depicted in  FIG. 8B  with a penalty factor applied thereto according to one or more embodiments described and illustrated herein; 
           [0024]      FIGS. 10A-10D  depict isoparametric shape functions having various shapes according to embodiments described and illustrated herein; 
           [0025]      FIG. 11  is a flowchart depicting a process to obtain convergence in an objective function and convergence in discrete angle by constraint transformation using an isoparametric function and penalty factor according to one or more embodiments described and illustrated herein; 
           [0026]      FIG. 12A  schematically depicts a magnetic fluid cooling device with magnet elements arranged at an initial position prior according to one or more embodiments described and illustrated herein; 
           [0027]      FIG. 12B  graphically depicts a simulated temperature profile of the magnetic fluid cooling device with the magnet elements arranged as depicted in  FIG. 12A  according to one or more embodiments described and illustrated herein; 
           [0028]      FIG. 12C  graphically depicts simulated magnetic fluid velocity magnitude contours within the magnetic fluid cooling device with the magnet elements arranged as depicted in  FIG. 12A  according to one or more embodiments described and illustrated herein; 
           [0029]      FIG. 12D  graphically depicts simulated magnetic field patterns of the magnetic fluid cooling device depicted in  FIG. 12A  according to one or more embodiments described and illustrated herein; 
           [0030]      FIG. 12E  graphically depicts an objective function value of the magnetic fluid cooling device depicted in  FIG. 12A  and a penalty factor according to one or more embodiments described and illustrated herein; 
           [0031]      FIG. 12F  graphically depicts a maximum temperature curve of the magnetic fluid cooling device depicted in  FIG. 12A  after several iterations of a topology optimization according to one or more embodiments described and illustrated herein; 
           [0032]      FIG. 13A  schematically depicts a magnetic fluid cooling device with magnet elements arranged in discrete angles after convergence of an objective function and convergence in discrete angle according to one or more embodiments described and illustrated herein; 
           [0033]      FIG. 13B  graphically depicts a simulated temperature profile of the magnetic fluid cooling device with the magnet elements arranged as depicted in  FIG. 12A  according to one or more embodiments described and illustrated herein; 
           [0034]      FIG. 13C  graphically depicts simulated magnetic fluid velocity magnitude contours within the magnetic fluid cooling device with the magnet elements arranged as depicted in  FIG. 13A  according to one or more embodiments described and illustrated herein; 
           [0035]      FIG. 13D  graphically depicts simulated magnetic field patterns of the magnetic fluid cooling device depicted in  FIG. 13A  according to one or more embodiments described and illustrated herein; 
           [0036]      FIG. 14A  graphically depicts an eight pointed star-shaped isoparametric shape function according to one or more embodiments described and illustrated herein; 
           [0037]      FIG. 14B  schematically depicts a magnetic fluid cooling device with magnet elements arranged in accordance with the eight pointed star-shaped isoparametric shape function depicted in  FIG. 14A ; 
           [0038]      FIG. 15A  graphically depicts a six pointed star-shaped isoparametric shape function according to one or more embodiments described and illustrated herein; 
           [0039]      FIG. 15B  schematically depicts a magnetic fluid cooling device with magnet elements arranged in accordance with the six pointed star-shaped isoparametric shape function depicted in  FIG. 15A ; 
           [0040]      FIG. 16A  graphically depicts a four pointed star-shaped isoparametric shape function according to one or more embodiments described and illustrated herein; 
           [0041]      FIG. 16B  schematically depicts a magnetic fluid cooling device with magnet elements arranged in accordance with the four pointed star-shaped isoparametric shape function depicted in  FIG. 16A ; 
           [0042]      FIG. 17A  graphically depicts a π/4 rotated four pointed star isoparametric shape function according to one or more embodiments described and illustrated herein; and 
           [0043]      FIG. 17B  schematically depicts a magnetic fluid cooling device with magnet elements arranged in accordance with the π/4 rotated four pointed star isoparametric shape function depicted in  FIG. 17A . 
       
    
    
     DETAILED DESCRIPTION 
       [0044]    Referring generally to the appended figures, embodiments of the present disclosure are directed to nonlinear programming methods for arranging material physical properties of devices or systems. As used here, the phrase “material physical property” means a property of the device or system that is being optimized, such as the arrangement of a magnet and the corresponding direction of its magnetic field. Embodiments ease handling and convergence of nonlinear programming problems with a set of pairs or triplets of design values, which usually represent vector values in Cartesian coordinates involving magnitude and direction. Embodiments also utilize an isoparametric shape function with application of a penalty factor to promote convergence of such directional material physical property design values into a set of specified desired options (e.g., angles of permanent magnets or arrangement of fibers of layers of a composite structure) to avoid unwanted intermediate angles in the final design. Various embodiments of nonlinear programming methods will be described in more detail herein. 
         [0045]    Referring now to  FIG. 1 , embodiments described herein may be utilized to calculate discrete angles of magnet elements  111  of a magnetic fluid cooling device  100 . Although embodiments described herein are described in the context of arranging magnet elements  111  of a magnetic fluid cooling device  100 , embodiments are not limited thereto. Embodiments may be used to determine the orientation of physical material properties of any device or system. As an example and not a limitation, embodiments described herein may be used to determine the direction of fibers in various layers of a composite material structure. 
         [0046]    As described in more detail below, the magnetic fluid cooling device  100  is configured to generate a magnetic field that disrupts thermal boundary layers that may develop within coolant fluid located in a magnetic fluid chamber  112 . Embodiments utilize a magnetic fluid in a thermo-magnetic cooling system to reduce the temperature field due to non-uniform heat flux caused by a heat generating device  120 , such as a semiconductor device. More particularly, embodiments comprise a magnetic fluid chamber  112  in which a magnetic fluid flows, and magnetic field generating devices (i.e., magnet elements  111 ) capable of generating magnetic fields that disrupt the flow of the magnetic fluid and therefore preventing, disrupting and/or altering natural fluid and thermal boundary layers within the magnetic fluid. Disruption of the fluid and thermal boundary layers may increase the performance of the cooling structure. Exemplary magnetic fluid cooling devices  100  are further described in U.S. patent application Ser. No. 13/316,954 entitled “Magnetic Fluid Cooling Devices and Power Electronics Assemblies” filed on Dec. 12, 2011, which is hereby incorporated by reference in its entirety. 
         [0047]      FIG. 1  schematically depicts a magnetic fluid cooling device  100  according to one embodiment. The magnetic fluid cooling device  100  may be coupled to a heat generating device  120  (e.g., a power semiconductor device, a motor, etc.) to remove heat flux q generated by the heat generating device. The magnetic fluid cooling device  100  generally comprises a magnetic fluid chamber  112  positioned between a first permanent magnet array  110 A and a second permanent magnet array  110 B. The first and second permanent magnet arrays  110 A comprise magnet elements  111  configured as individual permanent magnets. As described in detail below, each magnet element  111  is arranged at a particular angle within the array to achieve a desired magnetic field pattern to optimally cool the heat generating device  120 , which may be thermally coupled to a heat receiving surface of the magnetic fluid cooling device  100 . It should be understood that alternative embodiments may include more or fewer permanent magnet arrays, as other configurations are also possible. 
         [0048]    Magnetic fluid is provided in the magnetic fluid chamber  112 . The magnetic fluid may comprise any liquid that is susceptible to magnetic fields, and is capable of disrupting thermal boundary layers that are formed within the magnetic fluid by the application of the magnetic fields generated by the magnet elements  111 . The magnetic fluid may include, but is not limit to, magneto-rheological fluids, ferrofluids, and liquid metals. The mechanism through which the thermal and magnetic fields are coupled is known as the Curie effect, where the susceptibility, X, of the magnetic fluid increases as a function of its temperature, T. 
         [0049]      FIGS. 2-4  depict results of a simulation wherein a heat generating device  120  is cooled by magnetic fluid within a magnetic fluid chamber  112  as depicted in  FIG. 1 .  FIG. 2  depicts magnetic field lines defining a magnetic field pattern (shown as magnetic vector potential with units of Wb/m) generated by the first and second permanent magnet arrays  110 A,  110 B of a simulated magnetic fluid cooling device  100 , while  FIG. 3  depicts the resulting temperature profile (° C.).  FIG. 4  depicts the velocity magnitude contours of the magnetic fluid within the magnetic fluid chamber  112  resulting from the magnetic field patterns depicted in  FIG. 2  (m/s 2 ), and the heat flux generated by the heat generating device  120 . Colder magnetic fluid is drawn toward the heat generating device  120 , while warmer magnetic fluid is pushed away from the heat generating device  120 , thereby forming enhanced recirculation zones or flows within the magnetic fluid. These enhanced recirculation zones disrupt natural thermal boundary layers that typically form in liquid-cooled packages and improve the removal of heat flux from the heat generating device(s). 
         [0050]    The magnetic field patterns depicted in  FIG. 2  represent the optimal magnetic field patterns to provide the minimum surface temperatures of the temperature profile depicted in  FIG. 3  as determined by simulation. To achieve the optimal magnetic field patterns depicted in  FIG. 2 , the individual magnet elements  111  of the first and second permanent magnet arrays  110 A,  110 B should be arranged such that the magnetic north of each magnet element  111  is pointed at a proper angle.  FIG. 5  depicts an arrangement of the magnet elements  111  that yield the magnetic field patterns depicted in  FIG. 2  as determined by simulation. The arrangement depicted in  FIG. 5  is the result of gradient-based structural optimization. As shown in  FIG. 5 , the magnet elements  111  are arranged in continuous angles. However, such a magnetic fluid cooling device  100  with continuously angled magnet elements may be difficult to manufacture. It may be advantageous for manufacturing purposes to limit the number of possible angles for the permanent magnets to discrete, desirable angles (e.g., {0, π/2, π, 3 π/2} or {0, π/4, π/2, 3 π/4, π, 5 π/4, 3 π/2, 7 π/4}) (i.e., discrete angle design). However, it may be difficult to achieve convergence with polar based design values and penalty term approach. 
         [0051]    In the above-described design problem, the ideal magnitude and direction of the permanent magnet should be determined in two dimensions. There are two potential options in the selection of design values as shown in  FIGS. 6A and 6B . A first option, as shown in  FIG. 6A , is to use magnitude and angle as design values (polar form). A second option, as shown in  FIG. 6B , is to use vector component x and y as design values (Cartesian form). 
         [0052]    Polar form presents issues about the bound of angle. If the angle is bounded by 0 to 2 π radian, as shown in  FIG. 6A , the range between 0 and 2 π in a clockwise direction will be discontinuous even physically it is continuous. This may be an obstacle for gradient-based optimization because a variable on the bounds cannot move to a value across the bound. With Cartesian form, such discontinuity does not appear. Instead, a nonlinear design constraint must be introduced to keep the amplitude of the vector within given maximum value (e.g. unity, which is written as √(x2+y2)&lt;1, i.e., quadratic constraint). Because this nonlinear constraint has to be enforced to each pair of design values independently, it is difficult to implement when the problem had a large number of design values. 
         [0053]    As described in more detail below, in embodiments of the present disclosure, the latter Cartesian formulation is converted into a simpler problem with only side constraints. The x and y design values in the original problem are derived from newly introduced design values ξ and η on a normalized coordinate system with side constraints. An isoparametric shape function is used to map from ξ and η to x and y ( FIGS. 7A and 7B ). The nonlinear constraint to x and y will automatically be satisfied by choosing appropriate shape function.  FIGS. 7A and 7B  show that normalized design values ξ and η are projected onto a circular isoparametric shape function Nx(ξ, η), Ny(ξ, η). 
         [0054]    Embodiments also use isoparametric projection for promoting angular design values of material physical properties to converge into the set of desired angles. For example, the isoparametric projection shape function may be manipulated (via penalty factor) such that the magnetic elements  111  (or other components depending on the structure being designed) may be arranged at the desired angles ({0, 7 π/2, π, 3 π/2} or {0, π/4, π/2, 3 π/4, π, 5 π/4, 3 π/2, 7 π/4}). 
         [0055]    Referring now to  FIGS. 8A and 8B , squares drawn on a normalized coordinate system ξ, η (box constraint), and an isoparametric shape function to map the design values ξ, η to the magnetization space x, y (circular constraint) are respectively illustrated. The geometric transform from a box to a circle is performed by the isoparametric shape function. 
         [0056]    The isoparametric shape function depicted in  FIG. 8B  uses four quadratic eight node elements (nodes  130 - 137 ) to cover four quadrants A-D. Nodes  130 ,  132  and  134  represent discrete, desired angles for the magnetic elements. The illustrated isoparametric shape function provides for eight desired angles {0, π/4, π/2, 3 π/4, π, 5 π/4, 3 π/2, 7 π/4}. Intermediate nodes  131  and  133  represent undesirable angles between desirable angles provided by nodes  130 ,  132 , and  134 . Accordingly, the nodes alternate between discrete, desirable angle nodes and intermediate nodes. In some embodiments, the number of nodes is at least two times the number of discrete, desirable angles. 
         [0057]    To encourage the magnet elements  111  to be positioned in at least one of the desired angles (i.e., ξ, η converged into −1, 1, or 0), direct penalization by geometry transform of the isoparametric shape function is provided. Referring now to  FIG. 9 , the circular isoparametric shape function depicted in  FIG. 8B  is transformed into a star-shaped isoparametric shape function by application of a penalty factor c. With the star-shaped function, the amplitude of the vector (x, y) will have a high value at the vertices of the star when (ξ, η) is mapped to (x, y). Therefore, if the objective function has a characteristic to show better value (e.g., smaller) with the design values with larger amplitude in general, the optimizer will try to utilize that part. As a result, the design values converge at the vertices of the star. Thus, the vector design values x, y will have an angle close to the vertices of the star. When the penalty factor c=1, there is no penalty applied. A smaller c value yields a stronger penalty. Application of a smaller c value moves intermediate nodes  131  and  133  closer to (0, 0), as depicted by shifted intermediate notes  131 ′ and  133 ′.  FIGS. 10A-10D  depict isoparametric shape functions with c=1, c=0.9, c=0.8, and c=0.7, respectively. During topology optimization, which may be implemented by finite element methods, the magnet elements  111  will be encouraged to be oriented at one of the desired angles. 
         [0058]    Referring now to  FIG. 11 , a flowchart describing an exemplary process for solving a structural optimization problem (e.g., arrangement of magnet elements  111  to provide for magnetic field patterns that yield a desirable temperature profile of a magnetic fluid cooling device  100 ) is depicted. Generally, initial design values ξ, η are mapped to x, y using an isoparametric shape function. A topology optimization is performed to solve for the physical state of the system under analysis. If convergence in the objective function (e.g., a difference in temperature between iterations of the topology optimization) is not achieved, a gradient-based analysis is used to update the design values ξ, η, which are again mapped to x, y using the isoparametric shape function. The topology optimization is again performed and the objective function calculated. The process is repeated until convergence is achieved in the objective function. After convergence in the objective function occurs, the process checks to see if there is convergence in discrete angle (i.e., if each and every material physical property (e.g., magnet element  111 ) is arranged in one of the desired angles). If not, the shape of the isoparametric shape function is manipulated using penalization and the topology optimization is repeated until there is convergence in both the objective function and in discrete angle. 
         [0059]    The flowchart of  FIG. 11  will now be described in detail. At block  201 , initial material physical property design values ξ i  and η i  representing the direction of magnet elements  111  on the normalized coordinate system are received, and the coordinates for the isoparametric nodes are defined on a circular isoparametric shape function. At block  202 , the design values τ i  and η i  are projected to x i  and y i  using the isoparametric shape function in preparation of execution of a topology optimization in x, y space. Referring briefly to  FIG. 12A , the magnet elements  111  of the first and second permanent magnet arrays  110 A,  110 B are arranged at an initial angle (e.g., π/4). The orientation (i.e., direction of magnetization) of the magnet elements  111  is defined by ξ i  and η i , where i=0 to the total number of magnet elements  111  (or other material physical properties in applications different from the exemplary magnetic fluid cooler application described herein). 
         [0060]    At block  203 , a topology optimization is executed using the projected design values to x i , y i  to solve for the physical state of the system. In the present application, the physical state of the magnetic fluid cooling device  100  may include, but is not limited to, the magnetic field patterns generated by the current arrangement of the magnet elements  111  depending on the present iteration of the topology optimization, the magnetic fluid flow within the magnetic fluid chamber  112  resulting from the magnetic field patterns and the heat flux generated by the heat generating device  120 , and the temperature profile of the magnetic fluid cooling device  100  resulting from the magnetic fluid flow within the magnetic fluid chamber  112 . 
         [0061]      FIGS. 12B-12D  graphically depict, after one iteration of the topology optimization (i.e., the magnet elements  111  are arranged in the initial configuration depicted in  FIG. 12A ), a temperature profile (with temperature contour lines) of the magnetic fluid cooling device  100  ( FIG. 12B ) (° C.), magnetic fluid flow velocity magnitude contours within the magnetic fluid chamber  112  ( FIG. 12C ) (m/s 2 ), and magnetic field patterns  140 A,  140 B produced by the first and second permanent magnet arrays  110 A,  110 B ( FIG. 12D ) (wb/m). 
         [0062]    Referring once again to  FIG. 11 , at block  204  the objective function is calculated based on the previous iteration of the topology optimization. In the magnetic fluid cooling device application, the objective function may be based on the temperature of the magnetic fluid cooling device  100 . As an example and not a limitation, the objective function may be based on a difference in temperature between the present iteration of the topology optimization and the prior iteration. The objective function may be, without limitation, a different in the average temperature of the magnetic fluid cooling device or a difference in the maximum temperature of the magnetic fluid cooling device. 
         [0063]    A determination of whether there is convergence in the objective function is made at decision block  205 . Convergence in the objective function may be, for example, when the difference between the present calculation of the objective function and the prior calculation of the objective function is less than (or greater than, depending on the application) a threshold. As an example and not a limitation, convergence may be defined as being achieved when the calculation of the prior objective function and the present calculation of the objective function (e.g., the difference between the previous average temperature and the present average temperature) is within 0.5%. It should be understood that other metrics/thresholds may be used depending on the particular application. 
         [0064]    If there is not convergence in the objective function, a determination as to how to update the design values is made at blocks  206  and  207  such that the topology optimization may be executed again using updated design variables. At block  206 , sensitivities of the objective function with respect to x and y are calculated. The sensitivities of the objective function represent the slope of the previous objective function calculations. The sensitivities (i.e., slope) in x i  and y i  is converted back (i.e., reverse-translated) into ξ i  and η i  using the isoparametric shape function at block  207 . This retranslates the sensitivities of the design variables back onto the normalized coordinate system. 
         [0065]    At block  208 , ξ i  and η i  are updated using a gradient-based analysis. A determination of how much to change ξ i  and η i  to achieve convergence in the objective function is made based on the slope calculated at blocks  206  and  207 . A math programming method may be used to determine how to update ξ i  and η i  to achieve convergence in the objective function, such as the method of moving asymptotes (“MMA”). Accordingly, updated design variables ξ i  and η i  are created at block  208 . The process moves back to block  202 , where the updated design variables ξ i  and η i  are projected (i.e., translated) to x i  and y i  using the isoparametric shape function. Blocks  202 - 208  are repeated until convergence in the objective function is achieved at decision block  205 . 
         [0066]    An exemplary process depicted by blocks  202 - 208  will now be described mathematically. As stated above, vector design variables with Cartesian vector components are used. The following nonlinear programming problem, assumed to be solved by gradient based methods, is provided below: 
         [0067]    Minimize: f(p) 
         [0068]    Subject to:
       p={(x 1 ,y 1 ), (x 2 ,y 2 ), . . . , (x i , y i ), . . . , (x n , y n )}, and   √(x i   2 +y i   2 )&lt;1 for i=1, 2, . . . , n.
 
where p, x i , y i  are respectively the design variable set, and x and y component of each design variable. Note that the nonlinear constraints √(x i   2 +y i   2 )&lt;1 can be any other formula describing the bound which pair x i , y i  should reside within.
       
 
         [0071]    For partial differential equation (“PDE”) constrained optimization problems, such as structural optimization problems, it can be written as follows: 
         [0072]    Minimize: f(u,p)
       Subject to:   K(p)u=0,   p={p i |(x 1 ,y 1 ), (x 2 ,y 2 ), . . . , (x i , y i ), . . . , (x n , y n )}, and   √(x i   2 +y i   2 )&lt;1 for i=1, 2, . . . , n.
 
where K(p) is design dependent physical system matrix derived from PDE, such as stiffness matrix obtained by the finite element method (i.e., topology optimization) and u is an unknown variable vector of the equation that usually holds physical states, such as temperature or displacement.
       
 
         [0077]    Rather than solving the above problem directly, the design variables ξ, η are introduced on a normalized coordinate system such that q={(ξ 1 , η 1 ), (ξ 2 , η 2 ), . . . , (ξ i , η i ), . . . , (ξ n , η n )}; p is formulated with q using the following relationship: 
         [0078]    x i =Nx(ξ i , η i ), and 
         [0079]    y i =Ny(ξ i , η i ), 
         [0000]    where Nx(ξ i , η i ) and Ny(ξ i , η i ) are the shape functions of the isoparametric shape function used in the finite element method analysis. By choosing the appropriate isoparametric shape function Nx(ξ i , η i ) and Ny(ξ i , η i ), √(x i   2 +y i   2 )&lt;1 for i=1, 2, . . . , n will be automatically satisfied with the following simple side constraints to ξ and η: 
         [0080]    −1&lt;ξ&lt;1, and 
         [0081]    −1&lt;η&lt;1. 
         [0082]    The entire optimization problem may be reformulated as follows: 
         [0083]    Minimize: f(u,q) 
         [0084]    Subject to:
       K(p)u=0,   p=N(q),   q={(ξ 1 ,η 1 ), (ξ 2 ,η 2 ), . . . , (ξ i ,η i ), . . . , (ξ n ,η n )},   −1&lt;ξ&lt;1, and   −1&lt;η&lt;1.       
 
         [0090]    The above problem may be solved by a nested approach, e.g., solving K(p)u=0 separately to obtain values of objective and constraint functions and their derivatives. Calculation cost of p=N(q) is almost ignorable and derivatives such as df/dq may be calculated using the chain rule df/dq=df/dp*dp/dq. The calculation cost of dp/dq is also ignorable. Thus, the original problem is transformed into a simpler problem without n of nonlinear constraints with only ignorable cost. 
         [0091]    As an example and not a limitation, an isoparametric shape function of eight node quadratic elements as described above may be used for the above problem with quadratic constraints. The quadratic constraints, √(x i   2 +y i   2 )&lt;1 for i=1, 2, . . . , n will be automatically satisfied by putting eight nodes along a unit circle to make the element boundary circle. For example, (cos(n π/4), sin(n π/4)) for n=0, 1, . . . , 7 may be an appropriate choice. 
         [0092]    Referring once again to  FIG. 11 , if there is convergence in the objective function at decision block  205 , the process checks to determine if there is convergence in discrete angle at decision block  209 . Convergence in discrete angle is defined as all of the material physical properties (e.g., the magnet elements  111 ) arranged in the discrete, desired angles described above. Convergence in discrete angle occurs when ξ, η=−1, 0, or 1 (see  FIG. 8A ). 
         [0093]    If convergence in discrete angle is not achieved, the isoparametric shape function is manipulated by applying a penalty factor c that moves the intermediate nodes  131 ,  133  (see  FIG. 8B ) to encourage the properties to be oriented in the desired angles. The penalty factor c changes the shape of the isoparametric shape function from circular as originally provided to a star shape. The shape of the modified isoparametric shape function depends on how many nodes are provided. 
         [0094]    After modification of the isoparametric shape function, the process returns to block  202  where the current design variables ξ, η are projected onto the modified isoparametric shape function. The topology optimization is again repeated until convergence in the objective function is achieved. If there is no convergence in the discrete angle once again at decision block  209 , the isoparametric shape function is modified once again and the entire process repeated until there is convergence in both the objective function and discrete angle, where the process ends at  211 . The result is an arrangement of the magnet elements  111  that provides for enhanced cooling performance while also easing assembly of the magnetic fluid cooling device  100 . 
         [0095]    Referring once again to  FIG. 5 , a magnetic fluid cooling device  100  having first and second permanent magnet arrays  110 A,  110 B with magnet elements  111  arranged such that there is convergence in the objective function (e.g., maximum or average temperature) after  20  iterations of the topology optimization but not convergence in discrete angle.  FIG. 12E  plots the objective function value as curve  300  and the penalty factor as curve  302 . Convergence in the objective function is achieved after 20 iterations in the present example. As shown in curve  304  of  FIG. 12F , the minimum maximum temperature of the magnetic fluid cooling device  100  (the basis for the objective function) occurs after about 20 iterations of the topology optimization. Although the arrangement of the magnet elements  111  depicted in  FIG. 5  may provide for the optimal solution, fabrication of such a magnetic fluid cooling device  100  may be difficult because many of the magnet elements  111  are not arranged in the desired angles. 
         [0096]    Thus, a penalty factor is applied to the isoparametric shape function to encourage convergence in discrete angle. In the example illustrated in  FIG. 12E , a penalty factor (c=0.9) is applied after 20 iterations. The penalty factor is depicted by curve  302 . As shown in  FIGS. 12E and 12F , there is a small performance drop with the application of the penalty factor. Embodiments of the present disclosure provide a balance between the optimal solution and ease of fabrication of the magnetic fluid cooling device  100 . 
         [0097]      FIG. 13A  schematically depicts a magnetic fluid cooling device  100  after one hundred iterations and application of a penalty factor (c=0.9). As shown in  FIG. 13A , convergence in discrete angle is achieved—each magnet element  111  is positioned in one of the discrete, desired angles.  FIGS. 13B-13D  depict the temperature profile in ° C. (with contour lines), the flow velocity magnitude contours of the magnetic fluid in m/s 2 , and the magnetic field patters  140 A′,  140 B′ resulting from the magnet element arrangement depicted in  FIG. 13A  (weber), respectively. There is only a small performance drop of the arrangement depicted in  FIG. 13A  compared to the arrangement depicted in  FIG. 5 . 
         [0098]    Accordingly, embodiments of the present disclosure enable the determination of an arrangement of material physical properties (e.g., magnetic fields of magnet elements, carbon fibers, etc.) for high performance while also considering manufacturing considerations (e.g., discrete angle design).  FIGS. 14A-17B  provide examples of variously shaped isoparametric shape functions and the arrangement of magnet elements that they may yield.  FIG. 14A  depicts an eight pointed star-shaped isoparametric shape function as described above, and  FIG. 14B  depicts the resulting arrangement of magnet elements  111 .  FIG. 15A  depicts a six pointed star shaped isoparametric shape function, and  FIG. 15B  depicts the arrangement of magnet elements  111  resulting from the isoparametric shape function depicted in  FIG. 15A .  FIGS. 16A and 16B  depict a four pointed star shaped isoparametric shape function and the resulting magnet element  111  arrangement, respectively.  FIGS. 17A and 17B  depicts a π/4 rotated four pointed star isoparametric shape function and the resulting magnet element  111  arrangement respectively. 
         [0099]    It should now be understood that embodiments described herein provide for conversion of nonlinear bound of vector design variables into side bounds using an isoparametric shape function. Embodiments also give a penalty to design variables arranged at undesirable angles by using a star-shaped isoparametric shape function. Embodiments may be utilized to determine an arrangement of material physical properties within a design. In one example, the methods described herein may be utilized to determine discrete angles of permanent magnets arranged in one or more arrays within a magnetic fluid cooling device. 
         [0100]    While particular embodiments have been illustrated and described herein, it should be understood that various other changes and modifications may be made without departing from the spirit and scope of the claimed subject matter. Moreover, although various aspects of the claimed subject matter have been described herein, such aspects need not be utilized in combination. It is therefore intended that the appended claims cover all such changes and modifications that are within the scope of the claimed subject matter.