Numerical Modelling - Systems Simulation - Case Study Example

ASSIGNMENT 2: NUMERICAL MODELLING (SYSTEMS MODELLING AND SIMULATION COMPUTATIONAL FLUID DYNAMICS (CFD) Computational Fluid Dynamics (CFD) canbe used to study the processes of lift, drag and thrust and the air motion involved in each. A Rolls-Royce Trent 900 Engine, with a three shaft turbofan containing rotors with three different levels of pressure, can be analysed with CFD using the method of meshing. The meshing method first breaks the solution domain into small elements or control volumes (constituting a mesh). So for example, when calculating the air flow around an aerofoil, the surrounding space is drawn with a triangular mesh. Each of these control volumes or elements must be in balance. The CFD solution contains the flow variables, or velocities and pressure of the air flow, at each point in the mesh. Conservation of mass and momentum equations are used to obtain these values of the flow variables. Each control volume is analyzed for the mass that flows into it and out of it, which creates a balanced equation. The flow variables along an edge of a control volume are averaged between the values of the mesh points on the endpoints of the edge. The CFD program FLUENT calculates the air flow around an aerofoil. Fluent convergence histories graphs and calculated streamlines for different incidence angles of the air flow are compared to experimental findings from a wind tunnel. In our CFD lecture, the pressure coefficient on the aerofoil surface is also compared between the calculated and measured findings for incidence angle of 17 and 15 degrees. Pressure fields and calculated streamlines were made for incidence angles of 10 and 17 degrees. The lift-to-drag ratio varies with incidence, and this affects when the aerofoil stalls. There are many more streamlines for angle of incidence of 17 degrees compared to the 10 degree case, and the aerofoil for 17 degrees incidence is stalled. The angle of incidence is also called the angle of attack that the wing makes; the lift force is perpendicular to the wing. The lift force is then directed both up and to the back. The backward component of this lift is called drag. A higher angle of incidence creates greater drag; so that the lift-to-drag ratio decreases as the angle of incidence increases. Examples of turbulence were then considered: these examples include tornado funnels; water vapor vortices from aircraft; a broken bridge due to a vortex system at the base; and thermal convection on a heated plate. This turbulence can be characterized by the Reynolds number, which is low for smooth laminar flow and high for turbulent flow. The regions of laminar and turbulent boundary layer flows were examined, as well as the transition boundary layer between them. The definition of turbulence and some of its effects were considered. Turbulent systems have large scale eddies which interact with the mean flow and feed smaller and smaller scale eddies until viscous dissipation converts the energy of the turbulence into heat. In order to model this turbulence, mean values and fluctuating components are used. Case Study 1 In the following, research for two short case studies illustrating the use of CFD is presented in the application areas of engines and external flow. The case studies used the CFD software Techplot 360. The case study of an engine used data from an Excel spreadsheet (shown in figure 2), which had columns for RPM (revolutions per minute); BMEP (brake mean effective pressure); and BSFC (brake specific fuel consumption). The data of figure 2 is graphed in figure 1 using Techplot 360; the graph shows a direct relationship between the brake mean effective pressure and the revolutions per minute. The BMEP reaches a maximum and falls to zero in a sawtooth pattern; this pattern repeats itself several times, with the BMEP maximum getting slightly larger. Fig.1-Graph of CFD relationships from engine data. Fig. 2- Engine data. This relationship can be studied further by varying the RPM, BSFC and BMEP data to see if the pattern changes in the BMEP vs. RPM plots. The goal is to determine what the optimal values are for all three variables within an engine, in order to obtain the most efficient engine performance. For instance, the goal may be to create a higher brake mean effective pressure at lower revolutions per minute. The three variables RPM, BMEP and BSFC can be plotted simultaneously in three dimensions in order to create a solid surface, which can be analysed by a CFD mesh. Using optimization techniques on the elements of the mesh, the desired pattern between variables can be found. Case Study 2 In the second case study, gas burner data was used. A 3D plot of this data is shown in figure 3, which the fluid is meshed into bands. In figure 4, the data is plotted using a finer mesh, created by triangularising the bands. The CFD software can be used to analyse where the maximum peaks of the gas are in the z-direction, and how these evolve in time. The graphs illustrate the wave nature of the gas, and are able to display the different contours of the fluid. Fluid variables can be adjusted so that the number of peaks in the z-direction is varied, and relationships between these different fluid variables and the fluid height can be found. Fig. 3- Gas burner CFD 3D plot. Fig.4- Gas burner showing the triangular elements for a CFD analysis. To do a CFD analysis of the gas burner fluid, each of the elements in figure 4 must be analysed so that the element is in balance, with the value of the variables entering the element equaling the values leaving the element. For instance, the momentum of the fluid (mass multiplied by velocity) entering the element must equal the momentum leaving the element. The mass entering and leaving the element must also be equal. The momentum conservation equations can then be averaged. There will be terms that contain the square of the fluctuating components. Averaging these equations again gives the cube of the fluctuating components. This process can be continued indefinitely, to get a better approximation. This is an example of modeling turbulence, which always contains approximation. 2. FINITE ELEMENT ANALYSIS (FEA) Finite element analyses can be categorized as structural or field analyses: structural includes static, modal and transient analyses; field includes thermal, acoustics, and electromagnetic analyses. The principle of stationary total energy states that a loaded structure will deform so that its shape minimises the total energy of the system. For static loading, the total energy equals the strain energy U plus the potential energy of the loads V: = U + V (static loading) The purpose of FEA code is to find the displacement of a node (= [a]) that minimises the total energy. A mesh design’s second order elements are more accurate than first order elements because the second order elements take into consideration nodes that are half way along the edges of the first order elements. In order to determine the number of elements that should be used, the mesh should be refined until the solution is converged upon. In the following, four short case studies illustrating the use of FEA in the areas of static stress analysis, modal analysis, and thermal analysis are described. For the first case study, the Carnegie Mellon University mini-FEA program was used. Case 1: Structural Deformations due to Forces and Displacements In this case, a rectangular object 10” by 4” is used, with material properties E = 1.5E5 psi and v = 0.3. The mesh will be made of 10 by 4 elements. A force of Fx = -100, 000, Fy = 0 is applied to the object on the middle node of the left hand side. The upper right and lower right node have displacements Ux = 0.2, Uy = 0 (horizontal motion to the right). After obtaining the solution, forces on the object can be drawn, and three equations of equilibrium can be shown to hold. The mesh is maximally deformed, and compared to the original mesh. Then the displacements at the point of the specified force are found; likewise, the forces at the specified displacements are found. The maximally deformed mesh has upper right hand corner forces: Fx = 5E4, Fy = 2.7E4; at the lower right hand corner forces are: Fx = 5E4 and Fy = -2.7E4. At the node on the left side where the force Fx = -1E5 was applied, the displacements are: Ux = -3.7, Uy = 1.5E-13. Case 2: Transient Heat Transfer in a One-Dimensional Rod In this case study, the FEA program WinFElt was used. The rod is modeled using two elements, so that there are three nodes. The initial temperature across the rod is 20 degrees. The left end of the rod (node one) is heated to 100 degrees and the right end of the rod is insulated; convection occurs across the length of the rod. The variables are given in figure 5. This is followed by the tabulated results, which shows how the temperature of the node on the right and the node in the middle of the element change over time. In the problem description, “values (1,150) (2,20)” means that there is a convection coefficient of 150 and a free stream temperature of 20. Since the node points 1 and 2 are used, the convection is over the surface area of the rod. If the same node point was used, then the convection would be over just one exposed node. Fig. 5: Variables for Transient Heat Analysis of a One-Dimensional Rod problem description title="One-Dimensional Transient Heat Analysis" nodes=3 elements=2 analysis=transient-thermal analysis parameters duration=25 dt=1 alpha=.66 mass-mode=consistent nodes=[2,3] nodes 1 x=0.0 constraint=given 2 x=0.1 constraint=free 3 x=0.2 rod elements 1 nodes=[1,2] material=copper load=surface 2 nodes=[2,3] load=surface material properties copper Ix=400 nu=375 rho=8900 Kx=401 c=24.44 A=0.2*0.2*3.14159 distributed loads surface values=(1,150) (2,20) constraints given Tx=100 iTx=20 free Tx=u iTx=20 end Figure 6- Results for the 1D Rod ---------------------------------------------------------------- time Tx(2) Tx(3) ------------------------------------------------------------------ 0 20 20 1 15.971 20.404 2 33.931 16.822 3 43.352 22.768 4 50.032 30.338 5 55.515 37.614 6 60.252 44.19 7 64.407 50.036 8 68.071 55.21 9 71.305 59.781 10 74.16 63.819 11 76.682 67.385 12 78.91 70.535 13 80.877 73.317 14 82.614 75.774 15 84.148 77.944 16 85.504 79.86 17 86.7 81.553 18 87.758 83.048 19 88.691 84.368 20 89.516 85.534 21 90.244 86.564 22 90.887 87.474 23 91.455 88.277 24 91.957 88.986 25 92.4 89.613 Case 3: Modal Analysis of a Two-Dimensional Frame This modal analysis studies the deformations of a frame consisting of a high beam and a middle beam supported by two columns. In this case, the WinFElt program was used, in order to calculate the eigenvalues, or natural frequencies, and eigenvectors, or mode shapes. An input file is created which contains the input variables; in WinFElt, this is also called the problem description. The output is a table of eigenvalues and eigenvectors. Using a tool called Velvet, the mode shapes can also be drawn. The material properties listed in the problem description are: the cross-sectional area (A), the Ixx moment of inertia (Ix), Young’s modulus (E), and density (rho). When the contraints are unconstrained they equal “u”; constrained is denote “c”. The Tx constraint is the translation in the x direction; the Rx constraint is the rotation in the x-direction. When 6 elements and 6 nodes were used, 12 frequencies were found, ranging from 6 to 227 rad/s. When the number of nodes and elements were increased to 8 each, so that there was an extra node on each of the lower parts of the columns, the results show 18 different frequencies ranging from 6 to 488 rad/s. Figures 4 and 5 show the output for 8 elements/nodes case, and figure 6 shows the configuration of the frame and its deformation. Figure 7- Input file for the modal analysis of a 2D frame. problem description title = "Modal Analysis Part 2" analysis=modal nodes=8 elements =8 analysis parameters mass-mode = lumped nodes 1 x=0.0 y=0.0 constraint=fixed 2 x=325 y=0 3 x=0 y=150 constraint=free 4 x=325 5 x=0 y=300 6 x=325 7 x=0 y=75 8 x=325 beam elements 1 nodes=[1,7] material=column 2 nodes=[7,3] 3 nodes=[3,5] 4 nodes=[2,8] 5 nodes=[8,4] 6 nodes=[4,6] 7 nodes=[3,4] material=beam 8 nodes=[5,6] material properties column A=20 Ix=450 E=3e6 rho=0.005 beam A=10 Ix=500 E=3e6 rho=0.005 constraints fixed Tx=c Ty=c Rz=c free Tx=u Ty=u Rz=u end Figure 8-Solutions for the modal analysis of the 2D frame ** Modal Analysis Part 2 ** Modal frequencies (rad/sec) ------------------------ Mode # Frequency ------------------------ 1 6.2582 ( 0.99602 Hz) 2 17.451 ( 2.7774 Hz) 3 20.719 ( 3.2975 Hz) 4 28.476 ( 4.5321 Hz) 5 31.168 ( 4.9605 Hz) 6 37.998 ( 6.0476 Hz) 7 86.513 ( 13.769 Hz) 8 94.063 ( 14.971 Hz) 9 94.396 ( 15.024 Hz) 10 107.61 ( 17.126 Hz) 11 110.16 ( 17.533 Hz) 12 117.61 ( 18.719 Hz) 13 204.58 ( 32.561 Hz) 14 205.18 ( 32.655 Hz) 15 231.09 ( 36.78 Hz) 16 231.21 ( 36.798 Hz) 17 488.48 ( 77.745 Hz) 18 488.49 ( 77.746 Hz) Mode shapes ------------------------------------------------------------------------------ Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 Mode 6 ------------------------------------------------------------------------------ 1 1 1 1 1 1 0.0052459 -0.0058951 2.3861e-14 -0.00098816 -1.0502e-14 0.032196 -0.0072943 -0.00099593 -0.22453 0.023137 0.76352 -0.037547 1 1 -1 1 -1 1 -0.0052459 0.0058951 -1.5667e-14 0.00098816 -2.2704e-15 -0.032196 -0.0072943 -0.00099593 0.22453 0.023137 -0.76352 -0.037547 2.2068 -0.44774 -0.60259 -0.48338 -2.2357 -1.5708 0.007387 -0.012135 3.7279e-14 0.007921 -1.8428e-14 0.04601 -0.0050446 0.014381 0.54369 -0.020372 0.33186 -0.026615 2.2068 -0.44774 0.60259 -0.48338 2.2357 -1.5708 -0.007387 0.012135 -2.9951e-14 -0.007921 -4.667e-15 -0.04601 -0.0050446 0.014381 -0.54369 -0.020372 -0.33186 -0.026615 0.36463 0.49608 -3.8722 1.0141 16.369 -0.2375 0.0026234 -0.0029517 7.0963e-15 -0.00049596 -1.5722e-14 0.016208 -0.0081843 -0.0098241 0.04662 -0.016103 -0.20576 -0.0006356 0.36463 0.49608 3.8722 1.0141 -16.369 -0.2375 -0.0026234 0.0029517 -1.4044e-14 0.00049596 -2.0847e-14 -0.016208 -0.0081843 -0.0098241 -0.04662 -0.016103 0.20576 -0.0006356 Mode 7 Mode 8 Mode 9 Mode 10 Mode 11 Mode 12 ------------------------------------------------------------------------------ 1 1 1 1 1 1 7.3134e-16 -3.736e+16 -230.12 0.0083927 -1.5078e-14 3.7099e-18 -0.003114 -0.10997 -0.053214 0.0089176 0.021111 0.0024607 -1 -57.003 1 1 -1 -1 7.4925e-16 -3.736e+16 230.12 -0.0083927 -7.8531e-15 -1.6579e-16 0.003114 0.98691 -0.053214 0.0089176 -0.021111 -0.0024607 0.084919 8.6429 -0.41338 -0.0094675 34.476 -0.15086 9.2526e-16 -5.7091e+16 -351.52 0.015418 -1.6969e-14 2.0836e-16 -0.00064756 0.69196 -0.087813 -0.00021598 0.012514 -0.00032499 -0.084919 78.357 -0.41338 -0.0094675 -34.476 0.15086 9.7601e-16 -5.7091e+16 351.52 -0.015418 -9.5248e-15 -6.7464e-16 0.00064756 0.3464 -0.087813 -0.00021598 -0.012514 0.00032499 1.6411 52.423 -3.8343 -5.1026 -4.8379 -1.5566 2.8328e-16 -1.9488e+16 -120.08 0.0044372 -1.0245e-14 -1.5631e-16 -0.011283 -0.16872 0.0042221 -0.017049 -0.02171 -0.016028 -1.6411 -71.866 -3.8343 -5.1026 4.8379 1.5566 -5.6876e-17 -1.9488e+16 120.08 -0.0044372 -5.7106e-15 -2.3052e-16 0.011283 0.41559 0.0042221 -0.017049 0.02171 0.016028 Mode 13 Mode 14 Mode 15 Mode 16 Mode 17 Mode 18 ------------------------------------------------------------------------------ 1 1 1 1 1 1 -0.001901 2.555e-16 -2.3825e+16 9014.1 3.0647e+15 4.5195e+05 0.0057414 0.0044222 0.034787 0.34086 0.23022 3.72 1 -1 -1.0675 1 1.2686 1 0.001901 3.7763e-16 -2.3825e+16 -9014.1 3.0647e+15 -4.5195e+05 0.0057414 -0.0044222 0.045502 0.34086 0.26691 3.72 -0.0042507 -0.0069106 -2.9199 0.006563 -0.45638 0.31471 0.002957 -2.0746e-16 2.1936e+16 -8303.4 -3.6831e+14 -54330 -7.3242e-05 -7.952e-05 -0.078542 -0.31238 -0.036164 -0.45152 -0.0042507 0.0069106 3.2411 0.006563 -0.95243 0.31471 -0.002957 -3.973e-16 2.1936e+16 8303.4 -3.6831e+14 54330 -7.3242e-05 7.952e-05 0.11606 -0.31238 -0.017646 -0.45152 -0.19682 -0.18737 2.3919 -1.6327 -0.37416 -3.1498 -0.0011825 -1.5715e-16 -1.589e+16 6014.1 -1.2929e+16 -1.9062e+06 0.52356 0.39986 0.30353 0.31203 0.089339 0.19479 -0.19682 0.18737 -1.9824 -1.6327 0.58262 -3.1498 0.0011825 -1.3378e-15 -1.589e+16 -6014.1 -1.2929e+16 1.9062e+06 0.52356 -0.39986 -0.0044592 0.31203 0.026889 0.19479 Fig. 9: Drawing of the deformation of the frame (high and middle beam supported by two columns). Case 4: Static Structural Analysis of a Cantilever Beam This case study is an analysis of a cantilever beam, consisting of a steel beam projecting out horizontally, supported only on the left side, and under the load of its own weight. The beam is under an end load on the right side, which is a downward vertical force (Fy is negative). The beam contains one element, and two end nodes. The weight is distributed over the single element. These features are the input to the WinFElt program, as shown in figure 10. Fig. 10- Input file for structural analysis of the cantilever beam. problem description title = "Cantilever Beam" analysis=static nodes=2 elements=1 nodes 1 x=0.0 y=0.0 constraint=fixed 2 x=10.0 y=0.0 constraint=free force=end_load beam elements 1 nodes=[1,2] material=steel load=self_weight material properties steel A=10.0 E=30e6 Ix=357 forces end_load Fy=-1000 distributed loads self_weight direction=perpendicular values=(1,2000) (2,2000) constraints fixed Tx=c Ty=c Rz=c free Tx=u Ty=u Rz=u end From the solutions in fig. 11, the reaction forces include a translation force in the y-direction of node 1 (the fixed node) of Ty= -19000, which is a downward translation; there is also a rotation in the z-direction of Rz = -90000 at node 1. Node 2 is making a nodal displacement in the second and sixth degrees of freedom: DOF2 (= Ty, translation in the y-direction) and DOF6 (=Rz, rotation in the z-direction) degrees of freedom only. There are six degrees of freedom at each node (3 translation directions and 3 rotation directions). Node 1 makes no displacements at all, since it is fixed. The entire element is stressed in the z and negative y- rotational directions. Fig.11- Solution for the cantilever beam FEA. ** Cantilever Beam ** Nodal Displacements ----------------------------------------------------------------------------- Node # DOF 1 DOF 2 DOF 3 DOF 4 DOF 5 DOF 6 ----------------------------------------------------------------------------- 1 0 0 0 0 0 0 2 0 0.0002023 0 0 0 2.6455e-05 Element Stresses ------------------------------------------------------------------------------- 1: 0 -19000 -90000 0 -1000 7.276e-12 Reaction Forces ----------------------------------- Node # DOF Reaction Force ----------------------------------- 1 Tx 0 1 Ty -19000 1 Rz -90000 Finally, this case study was done for a wood cantilever beam. The material property changed was Young’s modulus, E=1.6 E6. The self-weight was also cut in half. In this case, the displacements of the second node were an order of magnitude larger. Stress on the element was to cause a rotation in the opposite z-direction. Reaction forces on the first node were only half as large. References Crocombe, Andrew. "Computational Fluid Dynamics." PowerPoint presentation for SE5305 Systems Modelling and Simulation. Crocombe, Andrew. "Introduction to FEA." PowerPoint presentation for SE5305 Systems Modelling and Simulation. Steif, Paul S. “Introduction to CMU Mini-FEA.” Retrieved 27 Apr 2007 from http://engineering-education.com/miniFEA/index.htm Techplot, Inc. "Techplot 360: The Complete CFD Visualization Tool." Retrieved 27 Apr 2007 from http://www.tecplot.com/products/360/360_main.htm “The FElt Home Page.” Retrieved 27 Apr 2007 from http://felt.sourceforge.net/ Read More