tag:blogger.com,1999:blog-13081399605160804122020-08-12T14:22:52.775-05:00Christopher SimpsonChristopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.comBlogger8125ChristopherSimpsonhttps://feedburner.google.comtag:blogger.com,1999:blog-1308139960516080412.post-1560618655606029462019-09-18T01:30:00.000-05:002019-10-03T14:33:31.462-05:00Opening the Space BoxRecently, I solved the high-value problem posed by Lockheed Martin's Space Box Challenge at my university. I was not prepared for what came next.<br /><br /><h2>What is the Space Box?</h2>Physically, it is a big, shiny, black box that looks like it came from a science fiction movie. The purpose of the box is to give university students a chance to demonstrate their problem-solving abilities. To get inside, you must solve one of four problems. There are three questions that are challenging, but they can be solved relatively easily with some thought and a little bit of math. Those three questions change daily for the three days the box was on campus. Upon solving a problem correctly, participants are rewarded by getting to go inside the box to view a cool space animation that surrounds you on all sides. Upon exiting the box, problem solvers receive a replica of the albums sent into space on the Voyager spacecraft. The fourth question, however, was much more difficult to solve than the other three. The reward was also much more substantial: A job offer from Lockheed Martin.<br /><div><br /><h2>Opening the Box</h2>As I mentioned, I solved the high-value problem, and I have since been offered a position working for Lockheed Martin. I am still in disbelief as I think about it right now. Many, many thanks to Lockheed Martin for bringing this challenge to my university and giving me this amazing opportunity. I am excited about the future.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-sNl_MEM_DGA/XZZMzlD3tYI/AAAAAAAACzg/QohKodJ4u5g4wUW0heN1IKQ25dSgQWrmwCLcBGAsYHQ/s1600/LM-Accept-OrionEarth%255B2%255D.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1000" data-original-width="1000" height="320" src="https://1.bp.blogspot.com/-sNl_MEM_DGA/XZZMzlD3tYI/AAAAAAAACzg/QohKodJ4u5g4wUW0heN1IKQ25dSgQWrmwCLcBGAsYHQ/s320/LM-Accept-OrionEarth%255B2%255D.png" width="320" /></a></div><br /></div><div><h2>The Problem and Solution</h2>Many people have been asking me about the problem I solved and how I solved it. I usually just answer, "with some math." For the more curious, here is the problem I solved:</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-TI-Orhc9vyg/XYBNbd3wCNI/AAAAAAAACzI/L64IM8iC2J01syDcUWRHslx7pTuNx-J4wCLcBGAsYHQ/s1600/spaceBoxChallengeQuestion.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="900" data-original-width="1600" height="360" src="https://1.bp.blogspot.com/-TI-Orhc9vyg/XYBNbd3wCNI/AAAAAAAACzI/L64IM8iC2J01syDcUWRHslx7pTuNx-J4wCLcBGAsYHQ/s640/spaceBoxChallengeQuestion.jpg" width="640" /></a></div><div><br /></div><div><br />Right away, this problem looks very intimidating. It seems to combine aspects of orbital mechanics and electromagnetism into a single problem. However, after you start to look at the problem more closely, it is only asking us to find a relationship between two reference frames. So let's look at how to tackle this problem.<br /><br />The problem basically defines a magnetic reference frame, and then asks us to find an expression for the angle at which the angular rate of rotation of the magnetic reference frame is double the orbital rate of the satellite (read: angular velocity).</div><div><br /></div><div>The only outside information needed to solve this problem is a dipole model of the Earth's magnetic field, which can be found in multiple forms across the internet, including <a href="https://en.wikipedia.org/wiki/Dipole_model_of_the_Earth%27s_magnetic_field">Wikipedia</a>. Because the expressions get quite long - and because I am lazy - I use the following shorthand:</div><div>\[s_\alpha=\sin\left(\alpha\right)\]</div><div>\[c_\alpha=\cos\left(\alpha\right)\]</div><div>Since the problem uses a reference angle measured from the equator, this is the form we want.</div><div>\[B_r=-\frac{2 B_0}{R^3}s_\alpha\]</div><div>\[B_\theta = \frac{B_0}{R^3}c_\alpha\]</div><div>\[|B| = \frac{B_0}{R^3} \sqrt{1+3s_\alpha^2}\]</div><div>where \(R\) is the distance from the center of the orbit divided by the radius of Earth. Note that these equations define a polar coordinate system which depends on the angular position within the orbit and both \(B_0\) and \(R\), neither of which are given. As it turns out, we don't need those values to solve the problem. I will define a set of unit vectors fixed on the satellite body, with \(\hat{e}_r\) pointing radially outward from the orbit and \(\hat{e}_\theta\) pointing along the trajectory of the satellite. The magnetic field is then given by</div><div>\[\vec{B} = \frac{B_0}{R^3}\left(-2s_\alpha \hat{e}_r + c_\alpha \hat{e}_\theta\right)\]</div><div>Using the given definition of the magnetic reference frame, we can compute the unit vector<br />\[\hat{x}_{mz} = \frac{\vec{B}}{|\vec{B}|} = \frac{-2s_\alpha}{\sqrt{1+3s_\alpha^2}}\hat{e}_r + \frac{c_\alpha}{\sqrt{1+3s_\alpha^2}}\hat{e}_\theta\]</div><div>We could calculate the others unit vectors also, but there is no need. We only need one unit vector to compute the angle of the reference frame with respect to the inertial reference frame. We need to express everything in terms of an inertial reference frame, because that is the reference frame in which the angular velocity of the satellite is defined. If we define the inertial frame with \(\hat{x}_1\) pointing to the right and \(\hat{x}_2\) pointing up, then we can write </div><div>\[\hat{e}_r = c_\alpha \hat{x}_1 + s_\alpha \hat{x}_2 \]</div><div>\[\hat{e}_\theta = -s_\alpha \hat{x}_1 + c_\alpha \hat{x}_2 \]</div><div>Substitution into the equation above (and a little bit of algebra) gives</div><div>\[\hat{x}_{mz} = \left(\frac{-3s_\alpha c_\alpha}{\sqrt{1+3s_\alpha^2}}\right)\hat{x}_1 + \left(\frac{-2s_\alpha^2+c_\alpha^2}{\sqrt{1+3s_\alpha^2}}\right)\hat{x}_2\]</div><div>Ultimately, we need an expression of the angle between the magnetic reference frame and the inertial reference frame so we can take a derivative. To get this angle, we simply take a dot product. Defining \(\beta\) as the positive angle between \(\hat{x}_1\) and \(\hat{x}_{mz}\),</div><div>\[ \cos(\beta) = \hat{x}_{mz} \cdot \hat{x}_1 = \frac{-3s_\alpha c_\alpha}{\sqrt{1+3s_\alpha^2}} \]</div><div>Take the inverse cosine, and we get</div><div>\[ \beta = \cos^{-1}\left(\frac{-3s_\alpha c_\alpha}{\sqrt{1+3s_\alpha^2}}\right) \]</div><div>Now we have our expression of the angle of the magnetic reference frame. All that's left to do is take a derivative. I used Wolfram Alpha for help with this, although it can be done by hand with a fair amount of patience.</div><div>\[\dot{\beta} = \frac{-3\dot{\alpha}\left(3 s_\alpha^4 + s_\alpha^2 - c_\alpha^2 \right)}{\left(3s_\alpha^2+1\right)\sqrt{3s_\alpha^2-9s_\alpha^2 c_\alpha^2 + 1}}\]</div><div>Note that \(\dot{\beta}=\omega_{mz}\) and \(\dot{\alpha} = \omega_0\), so what we have is an expression for \(\omega_{mz}\) in terms of \(\omega_0\). With some trigonometric identities and simplification, we get a slightly nicer expression.</div><div>\[\omega_{mz} = \frac{-3\left(\cos(2\alpha)-3\right)}{2\left(\sin^2(\alpha)+1\right)}\omega_0\]</div><div>The question asks for an expression for when the angular rate of rotation of the magnetic reference frame is exactly twice the orbit rate. In other words, when does that coefficient equal 2? There are a few ways to do this, but Wolfram Alpha again makes short work of the problem.</div><div>\[\frac{-3\left(\cos(2\alpha)-3\right)}{2\left(\sin^2(\alpha)+1\right)} = 2 \quad\textrm{when}\quad \alpha=\pi n \pm \frac{1}{2}\cos^{-1}\left(\frac{1}{3}\right) \quad\textrm{for}\quad\ n\in\mathbb{Z}\].<br />If you just need the first instance of the angle (which in this case also happens to be the magnitude angle from the horizon in all four quadrants) just take \(n=0\) and you get<br />\[\alpha = \frac{1}{2}\cos^{-1}\left(\frac{1}{3}\right)\]<br />That's it. It takes no more than a page or so of work, and it can be done in 30 minutes if you know what you are doing.</div><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/D_kjPO1lbaQ" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com0http://www.cdsimpson.net/2019/09/opening-space-box.htmltag:blogger.com,1999:blog-1308139960516080412.post-28525639409830973532019-01-09T23:50:00.001-06:002019-01-09T23:50:50.730-06:00Teensy TroubleIn my Ph.D. studies, I have been doing a fair amount of embedded programming. My favorite microcontroller by far has been the aptly-named <a href="http://www.pjrc.com/teensy">Teensy</a>. These microcontrollers have some serious hardware at a very reasonable price. They are also incredibly easy to program due to the fact that they can be programmed with arduino code. I could go on and on about why I like them, but today I found a relatively minor issue that caused me a lot of trouble.<br /><br />The Teensy 3.6 has many pins capable of reading analog voltages. The two pins I decided to use for my project are pins 32 (A13) and 35 (A16). I implemented the functionality using <code>analogRead(13)</code> and <code>analogRead(16)</code>, respectively. Interestingly, A13 worked as expected, but A16 did not. Instead, it behaved like a floating pin, and even shorting that pin to ground had no impact on the reading (this should have been a clue, actually). It wasn't just A16, either. All analog pins numbered A14 or higher were not working, while all numbered A13 or lower were working as expected.<br /><br />After much frustration and much investigation, I found the answer. Since <code>analogRead()</code> can accept either analog pin numbers or actual pin numbers as its argument, numbers over 14 are ambiguous to the code, so it assumes the argument is the actual pin number. To avoid this, you can either use the actual pin number or use the analog pin number preceded by "A" (i.e. A14, A15, etc.). This works because those are just macros to the actual pin numbers. In my case, using <code>analogRead(A16)</code> solved the problem.<img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/KjUT0JciK2Y" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com0http://www.cdsimpson.net/2019/01/teensy-trouble.htmltag:blogger.com,1999:blog-1308139960516080412.post-38019262589725277192016-09-27T09:26:00.000-05:002016-09-27T09:26:07.177-05:00Predicting Control Surface Hinge Moments using Computational Fluid DynamicsMy Master's thesis has finally been published! You can view my thesis <a href="https://drive.google.com/open?id=0B_OtEk0gkL6XMFZ5TjlxeHlCcTg">here</a> (faster) or on <a href="http://purl.lib.ua.edu/149910">The University of Alabama Electronic Theses and Dissertations database</a>. I started this research with absolutely no CFD experience. I am quite proud of the fact that I was able to learn enough that I am now considered by many to be an expert. When I look at how far I have come in the last three years, it truly amazes me. I never thought I would learn so much.<br /><br />Thanks to everyone who supported me while I worked on this research, and thanks for your continued support as I continue in my studies.<img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/tWWNJqyUUWA" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com0http://www.cdsimpson.net/2016/09/predicting-control-surface-hinge.htmltag:blogger.com,1999:blog-1308139960516080412.post-36880026957705343092016-09-03T18:03:00.000-05:002016-09-03T18:03:18.176-05:00Rebuilding a Three-Channel Aircraft (Dionysus)When I was a kid, I had a 3-channel Cub-like RC airplane. While in storage, the aircraft suffered some rather severe hangar rash. The damage was the result of the unexpected compression load applied to the fuselage due to poor visibility conditions. In other words, I stepped on the plane because it was dark.<br /><br /><h2>Bringing It Back From The Dead</h2>After I recently "retired" another aircraft when the wing incidence started varying mid-flight, I decided it was time to rebuild the old . The fuselage had been collapsed into itself, so I still had most of the broken pieces after all the years. This made reconstruction a great deal easier. Using copious amounts of wood glue and CA, and just a couple pieces of new balsa, the aircraft was ready for a new skin.<br /><br />I opted to replace all of the skin because the old skin would no longer adhere to the balsa. I removed the existing skin and gave the airframe a new skin of Ultracote Lite (transparent yellow). I found that it is easy to work with, it adheres well, and it gives a very nice finish. The "Lite" version is substantially thinner than the normal stuff.<br /><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="https://1.bp.blogspot.com/-Srzs8d2bZU8/V8st0dMV7yI/AAAAAAAABXc/C_DXcNex_Lk8pEcPgkvOINUtsUdCj_1SACLcB/s1600/20160903_113851.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="360" src="https://1.bp.blogspot.com/-Srzs8d2bZU8/V8st0dMV7yI/AAAAAAAABXc/C_DXcNex_Lk8pEcPgkvOINUtsUdCj_1SACLcB/s640/20160903_113851.jpg" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">The completed repair of the fuselage damage.</td></tr></tbody></table><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="https://1.bp.blogspot.com/-mUoyi2AwxgQ/V8st4Sm-FwI/AAAAAAAABX0/cDvUm6EqkkY-W9l6lbU7GfEEtSeWlccBwCEw/s1600/20160903_114248.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="360" src="https://1.bp.blogspot.com/-mUoyi2AwxgQ/V8st4Sm-FwI/AAAAAAAABX0/cDvUm6EqkkY-W9l6lbU7GfEEtSeWlccBwCEw/s640/20160903_114248.jpg" width="640" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">This spot on the wing had also been damaged. The broken piece was glued back in place.</td></tr></tbody></table><br /><br />In addition to structural repairs, I decided that it was a good time to modernize the components. Gone are the 72 MHz receiver, brushed motor, gearbox, ESC, and NiMh (NiCd?) battery. I now have a Spektrum 2.4 GHz receiver, brushless motor and ESC, and a 3s LiPo battery.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://3.bp.blogspot.com/-B2BWo5bHogU/V8st0jDrg2I/AAAAAAAABYA/AdJsuorvTO8oeUzGHH9x0RE9HLuwgviJwCEw/s1600/20160903_113930.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://3.bp.blogspot.com/-B2BWo5bHogU/V8st0jDrg2I/AAAAAAAABYA/AdJsuorvTO8oeUzGHH9x0RE9HLuwgviJwCEw/s400/20160903_113930.jpg" width="225" /></a></div><br /><br />I fabricated a motor mount to convert the aircraft from a brushed inrunner motor to a brushless outrunner. The entire motor mount assembly is attached to the airframe using a couple of cable ties through existing holes, so I can easily remove it if needed. With this motor, the aircraft has much more power than it needs. The thrust-to-weight ratio is about 2, so it can climb vertically.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://2.bp.blogspot.com/-lSJfPQGSWX0/V8stenipMoI/AAAAAAAABYA/XzOm2EGYFcYxdxKV-QcKmOrPPv_locLwQCEw/s1600/20160903_113955.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="320" src="https://2.bp.blogspot.com/-lSJfPQGSWX0/V8stenipMoI/AAAAAAAABYA/XzOm2EGYFcYxdxKV-QcKmOrPPv_locLwQCEw/s320/20160903_113955.jpg" width="180" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://4.bp.blogspot.com/-lpnl-Vjl0Ww/V8st3D4f2tI/AAAAAAAABYA/l4TNvGDe3YwpPbNHN1yWCwYbNgZFEsAmQCEw/s1600/20160903_114007.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://4.bp.blogspot.com/-lpnl-Vjl0Ww/V8st3D4f2tI/AAAAAAAABYA/l4TNvGDe3YwpPbNHN1yWCwYbNgZFEsAmQCEw/s320/20160903_114007.jpg" width="320" /></a></div><div class="separator" style="clear: both; text-align: center;"></div><br /><br /><br />I also added a Mobius ActionCam to record video of my flights. I mounted the camera on the bottom of the fuselage, directly at the center of gravity. A piece of 1/8-inch birch plywood on the inside of the fuselage acts as a reinforcement plate, distributing the load to the balsa fuselage. The camera mount simply bolts onto the fuselage. You can see some footage recorded with this camera below and on my YouTube channel.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://3.bp.blogspot.com/-fyk-FxOZtEA/V8st3F6JX_I/AAAAAAAABYA/OfSI7DdoDdQ7DCJlFAW-0QDV2G1Uk7MWQCEw/s1600/20160903_114052.jpg" imageanchor="1"><img border="0" height="320" src="https://3.bp.blogspot.com/-fyk-FxOZtEA/V8st3F6JX_I/AAAAAAAABYA/OfSI7DdoDdQ7DCJlFAW-0QDV2G1Uk7MWQCEw/s320/20160903_114052.jpg" width="180" /></a></div><br /><table align="center" cellpadding="0" cellspacing="0" class="tr-caption-container" style="margin-left: auto; margin-right: auto; text-align: center;"><tbody><tr><td style="text-align: center;"><a href="https://2.bp.blogspot.com/-9gIaoDeTLOM/V8st3bqeGRI/AAAAAAAABYA/ahT8CuVRVV87NDg715b06xIB-Bm78kESwCEw/s1600/20160903_114129.jpg" imageanchor="1" style="margin-left: auto; margin-right: auto;"><img border="0" height="180" src="https://2.bp.blogspot.com/-9gIaoDeTLOM/V8st3bqeGRI/AAAAAAAABYA/ahT8CuVRVV87NDg715b06xIB-Bm78kESwCEw/s320/20160903_114129.jpg" width="320" /></a></td></tr><tr><td class="tr-caption" style="text-align: center;">You can see the reinforcement plate for the camera at the bottom of the fuselage.</td></tr></tbody></table><br /><br />The total weight of the aircraft with battery and camera installed is 1 lb 15 oz. Because I used parts that I already had lying around, the total cost of the rebuild was about $5 plus my time. Overall, I am pleased with the quality of the rebuild.<br /><br /><h2>Test Flight Footage and Observations</h2><div>Of course, the point of undertaking this repair was to get the aircraft flying again. So early one morning, I took it to the soccer field for a test flight. It turns out that 8:00 AM is a good time to fly for a number of reasons. Not only is it the most convenient time for me, but also the soccer fields are not in use, the weather is usually calm, and the lighting is perfect for filming to the west.</div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><iframe width="320" height="266" class="YOUTUBE-iframe-video" data-thumbnail-src="https://i.ytimg.com/vi/b0bSux1TyWI/0.jpg" src="https://www.youtube.com/embed/b0bSux1TyWI?feature=player_embedded" frameborder="0" allowfullscreen></iframe></div><div><br /></div><div><br /></div><div>The aircraft itself performs very well. Since the aircraft has no ailerons, it relies on the dihedral of the wing to produce a rolling moment with rudder deflection. I mapped the rudder to the aileron channel on my transmitter. The dynamic roll response is surprisingly fast. Truth be told, I only notice the lack of ailerons during slow flight.</div><div><br /></div><div>With the salvaged parts I had on hand, the aircraft has much more power than it needs. I estimate the thrust-to-weight ratio to be about 2, so it can climb vertically with no problem. The recorded flight lasted 10 minutes, and the 3S 2200 mAh battery was still at 11.8 V at the end of the flight.</div><div><br /></div><div>Having successfully resurrected the aircraft, it is only appropriate to give it a suitable name. I like to name things after mythological entities, and I think Dionysus fits pretty well given the story. I just hope I don't have to change the name to Osiris down the road.</div><div><br /></div><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/9J9tq-hZyUA" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com4http://www.cdsimpson.net/2016/09/rebuilding-three-channel-aircraft.htmltag:blogger.com,1999:blog-1308139960516080412.post-57388812852014564672015-04-30T00:17:00.001-05:002015-07-26T20:07:26.810-05:00Temporal and Spatial Stability Analysis of the Orr-Sommerfeld EquationThis is the second and final part of the stability analysis of the Orr-Sommerfeld Equation. In my <a href="http://www.cdsimpson.net/2015/04/derivation-of-orr-sommerfeld-equation.html">previous post</a>, I went through the derivation and nondimensionalization of the Orr-Sommerfeld Equation. In this part, I show how to perform the temporal and spatial stability analyses.<br /><br /><h2>Blasius Velocity Profile</h2>A 2-D Blasius boundary layer can be expressed as<br />$$\bar{U}=f'\left(\eta\right)\tag{1}\label{blasius}$$<br />where<br />$$\begin{align}\eta=y\sqrt{\frac{U_{\infty}}{2\nu x}}\\f'''+ff''&=0\\f(0)=f'(0)&=0\\f'(\infty)&=1\end{align}$$<br />This is a nonlinear ordinary differential equation that is most easily solved using a shooting method. In a shooting method, the boundary value problem is made into an initial value problem, whereby an initial guess of one parameter is iteratively updated until the opposite boundary conditions are met. In the present case, the unknown boundary condition is \(f''(0)\) , so it is iteratively adjusted until \(|f'(10)-1|\leq10^{-6}\). The present case uses a fourth-order Runge-Kutta method of integration. The initial condition The transformation of the Blasius velocity profile from \(\eta\) to \(\xi\) is performed by determining the distance from the wall, \(\delta\), where the stream-wise velocity is \(99.9\%\) of the free-stream velocity.<br />$$\begin{align}\bar{U}\left(\xi\right)&=\bar{U}\left(\frac{\eta}{\delta}\right)\\\bar{U}''\left(\xi\right)&=\frac{1}{\delta^{2}}\bar{U}''\left(\frac{\eta}{\delta}\right)\end{align}\tag{2}\label{transform}$$<br />The figure below shows the Blasius velocity profile and its derivatives.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-Jak55SmZBnI/VUGajd2rqmI/AAAAAAAAAX4/-wFfZ9i06kc/s1600/blasius.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="351" src="http://2.bp.blogspot.com/-Jak55SmZBnI/VUGajd2rqmI/AAAAAAAAAX4/-wFfZ9i06kc/s1600/blasius.png" width="520" /></a></div><br /><h2>Temporal Stability</h2><div><div>The Orr-Sommerfeld equation can be analyzed for temporal stability by assuming that \(\bar{\alpha}\) is real. This enables us to rearrange the non-dimensional Orr-Sommerfeld equation as follows.</div></div><div>$$\left(-\bar{U}\bar{\alpha}^{2}-U''+\frac{i\bar{\alpha}^{3}}{Re_{\delta}}\right)\bar{\phi}+\left(\bar{U}-\frac{2i\bar{\alpha}}{Re_{\delta}}\right)\bar{\phi}''+\left(\frac{i}{\bar{\alpha}Re_{\delta}}\right)\bar{\phi}''''=\bar{c}\left(\bar{\phi}''-\bar{\alpha}^{2}\bar{\phi}\right)\tag{3}\label{realalpha}$$</div><div><div>If we discretize this equation using a finite difference approximation for \(\bar{\phi}''\) and \(\bar{\phi}''''\), we get an equation of the form</div></div><div>$$A\Phi=\tilde{c}B\Phi\tag{4}\label{eigmatrixform}$$</div><div><div>where \(\Phi\) is a vector containing values \(\bar{\phi}_{i}\) at discrete locations. This is nothing more than an eigenvalue problem, where \(\tilde{c}\) is a diagonal matrix of eigenvalues for the discrete system. The present calculations use a second order central differencing scheme to approximate the second and fourth derivative terms.</div></div><div>$$\begin{align}\bar{\phi}_{i}''&\approx\frac{\bar{\phi}_{i-1}-2\bar{\phi}_{i}+\bar{\phi}_{i+1}}{h^{2}}\\\bar{\phi}_{i}''''&\approx\frac{\bar{\phi}_{i-2}-4\bar{\phi}_{i-1}+6\bar{\phi}_{i}-4\bar{\phi}_{i+1}+\bar{\phi}_{i+2}}{h^{4}}\end{align}\tag{5}\label{eqfindiff}$$</div><div><div>At the boundaries, \(\bar{\phi}=0\), so the first and last columns of \(A\) and \(B\) are removed. However, \(\bar{\phi}'=0\) also at the boundaries. Fortunately, using a backward difference method at the boundary gives \(\bar{\phi}_{-1}\approx\bar{\phi}_{0}\), so the \(\bar{\phi}_{i-2}\) term can be omitted from the \(\bar{\phi}''''\) approximation just inside the wall boundary. Similarly, the \(\bar{\phi}_{i+2}\) term can be omitted just inside the free-stream boundary.</div></div><div><div>The finite difference approximations in Equation \ref{eqfindiff} gives a pentadiagonal \(A\) and tridiagonal \(B\).</div></div><div><div>$$\begin{align}A&=\frac{1}{h^{4}}\left[\begin{array}{ccccc}a_{11} & a_{21} & a_{3} & & 0\\a_{21} & \ddots & \ddots & \ddots\\a_{3} & \ddots & \ddots & \ddots & a_{3}\\ & \ddots & \ddots & \ddots & a_{2n}\\0 & & a_{3} & a_{2n} & a_{1n}\end{array}\right]\\B&=\frac{1}{h^{2}}\left[\begin{array}{cccc}b_{1} & b_{2} & & 0\\b_{2} & \ddots & \ddots\\ & \ddots & \ddots & b_{2}\\0 & & b_{2} & b_{1}\end{array}\right]\end{align}\tag{6}\label{fdmatrix}$$</div></div><div>where</div><div>$$\begin{align}a_{1i}&=h^{4}\left(-\bar{U}_{i}\bar{\alpha}^{3}-\bar{U}_{i}''\bar{\alpha}+\frac{i\bar{\alpha}}{Re_{\delta}}\right)-2h^{2}\left(\bar{U}_{i}\bar{\alpha}-\frac{2i\bar{\alpha}}{Re_{\delta}}\right)+6\left(\frac{i}{Re_{\delta}}\right)\\a_{2i}&=h^{2}\left(\bar{U}_{i}\bar{\alpha}-\frac{2i\bar{\alpha}}{Re_{\delta}}\right)-4\left(\frac{i}{Re_{\delta}}\right)\\a_{3}&=\frac{i}{Re_{\delta}}\\b_{1}&=-h^{2}\bar{\alpha}^{3}-\bar{\alpha}\\b_{2}&=\bar{\alpha}\end{align}\tag{7}\label{fdmatrixcoef}$$</div><div><div>The coefficients in \(A\) depend on the velocity profile, which varies with the distance from the wall. However, \(B\) is constant for a given \(\bar{\alpha}\). In addition, \(B\) is guaranteed to be real and symmetric, and is thus guaranteed to be invertible. We can thus left multiply Equation \ref{eigmatrixform} by \(B^{-1}\) and rewrite the equation as follows.</div></div><div>$$B^{-1}A\Phi=B^{-1}\tilde{c}B\Phi\tag{8}\label{newmatform}$$</div><div><div>It is clear that the eigenvalues of \(B^{-1}A\) are the eigenvalues in \(\tilde{c}\). Matlab's <span style="font-family: Courier New, Courier, monospace;">eig</span> function was used to solve for the eigenvalues for each combination of \(\bar{\alpha}\) and \(Re\). The following figure shows the eigenvalues of the discrete system using one value of \(\bar{\alpha}\) and \(Re\). Each combination of \(\bar{\alpha}\) and \(Re\) gives a unique set of eigenvalues similar to these.<br /><br /></div></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-EP0poIpgL_4/VUGqkXIYHAI/AAAAAAAAAYM/EH6VATkaAw8/s1600/temporal_eigenvalues_example.png" imageanchor="1"><img border="0" height="351" src="http://1.bp.blogspot.com/-EP0poIpgL_4/VUGqkXIYHAI/AAAAAAAAAYM/EH6VATkaAw8/s1600/temporal_eigenvalues_example.png" width="520" /></a></div><div><br />Because this is a stability analysis, we only care about the most unstable eigenvalue. The governing equation is related to \(e^{i\alpha(x-\left(c_{r}+ic_{i}\right)t)}\). If \(\alpha\) is real, then any instabilities must come from an eigenvalue with a positive imaginary component. For this reason, the interesting eigenvalue is the eigenvalue with the maximum imaginary component. This analysis was done over a range of \(\bar{\alpha}\) and \(Re\), and the most unstable eigenvalue was stored for each combination. Matlab's <span style="font-family: Courier New, Courier, monospace;">contour</span> function was then used to plot the contours shown in the following figure.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://3.bp.blogspot.com/-PBP5cxsygw0/VUGs1w2efQI/AAAAAAAAAYc/u0m60oZ-JxY/s1600/temporal_stability.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="351" src="http://3.bp.blogspot.com/-PBP5cxsygw0/VUGs1w2efQI/AAAAAAAAAYc/u0m60oZ-JxY/s1600/temporal_stability.png" width="520" /></a></div><br />It is important to note that the contours shown here differ from those found by Wazzan in [7] because the Orr-Sommerfeld equation has been nondimensionalized in a different manner. The present case used the boundary layer thickness \(\delta\), while Wazzan used displacement thickness \(\delta^{*}\).<br /><br /><h2>Spatial Stability</h2></div><div><div>The temporal stability analysis is performed assuming that \(\bar{\alpha}\) is real. In the spatial stability analysis, we assume that \(\omega\) is real. It is possible to derive a polynomial eigenvalue problem for the eigenvalues \(\bar{\alpha}\), but solving the polynomial eigenvalue problem presented some challenges. Namely, at least one of the coefficient matrices was singular.</div></div><div><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-qpxXjRGwBwo/VUGtuukLK2I/AAAAAAAAAYk/IuERClemElQ/s1600/spatial_map.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="228" src="http://1.bp.blogspot.com/-qpxXjRGwBwo/VUGtuukLK2I/AAAAAAAAAYk/IuERClemElQ/s1600/spatial_map.png" width="520" /></a></div><div><br />To circumvent this challenge, a mapping method has been used to give a relation between a complex \(\bar{\alpha}\) and a complex \(\omega\). In the present case, we hold \(\bar{\alpha}_{i}\) at a specific value and vary \(\bar{\alpha}_{r}\). For each \(\bar{\alpha}_{r}\), the discrete eigenvalue problem is solved for \(\omega\), and the most unstable eigenvalue is kept as before. By varying \(\bar{\alpha}_{r}\) with \(\bar{\alpha}_{i}\), a curve of \(\omega\) values can be plotted in the complex plane, as shown in the figure above. We then determine if there is an \(\bar{\alpha}_{r}\) that gives \(\omega_{i}=0\) and store any locations that satisfy this condition. This process is repeated for a range of \(Re\) for all interesting \(\bar{\alpha}_{i}\) values. The results of this process are shown in the figure below.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://1.bp.blogspot.com/-16e5jaNRuPs/VUGujbnrmLI/AAAAAAAAAY0/PyGt6anvPz0/s1600/spatial_stability.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="351" src="http://1.bp.blogspot.com/-16e5jaNRuPs/VUGujbnrmLI/AAAAAAAAAY0/PyGt6anvPz0/s1600/spatial_stability.png" width="520" /></a></div><br /><h2>Discussion</h2></div><div><div>This report is not the first on this topic, and the results do match fairly well with the references. Once the discretization was properly arranged, the eigenvalue computation was trivial using Matlab's built-in functions. Some factors added complexity to the project. The temporal stability curves are sensitive to the discretization step size and the number of \(Re\) and \(\bar{\alpha}\) values used. However, the spatial stability curves were far more sensitive. </div><div><br /></div><div>It turns out that an insufficient number of \(\bar{\alpha}_{r}\) values will give bad resolution in the complex \(\omega\) plane, periodically overestimating the zero locations. This caused some oscillations at larger values of \(Re\). Additionally, a large number of Reynold's Number steps were required to ensure that the each of the desired contours were visible.Compounding the issue was the fact that an extra variable effectively raised computation time by nearly an order of magnitude.<br /><br /></div></div><h2>References</h2><div><ol><li>J. D. Anderson. <i>Fundamentals of Aerodynamics</i>. Aeronautical and Aerospace Engineering Series. McGraw Hill, 4th edition, 2007.</li><li>M. J. Maghrebi. Orr Sommerfeld Solver using Mapped Finite Di fference scheme for plane wake fl ow. <i>Journal of Aerospace Science and Technology</i>, 2(4):55-63, December 2005</li><li>P. J. J. Moeleker. Linear Temporal Stability Analysis. Technical report, Delft University of Technology, 1998. Series 01: Aerodynamics 07.</li><li>B. S. Ng and W. H. Reid. Simple Asymptotics for the Temporal Spectrum of an Orr-Sommerfeld Problem. <i>Applied Mathematics Letters</i>, 13:51-55, 2000.</li><li>S. A. Orszag. Accurate Solution of the Orr-Sommerfeld Stability Equation. Technical report, Massachusetts Institute of Technology, May 1971.</li><li>T. Patel. Stability Properties of Self-propelled Wakes. Master's Thesis, University of California, San Diego, 2012.</li><li>A. R. Wazzan, T. T. Okamura, and A. M. O. Smith. Spatial and Temporal Stability Charts for the Falkner-Skan Boundary-layer Profiles. Technical report, Douglas Aircraft Company, 1968.</li><li>F. M. White. <i>Viscous Fluid Flow</i>. McGraw-Hill Series in Mechanical Engineering. McGraw Hill, 2nd edition, 1991.</li></ol><div><br /></div><h2>Matlab Codes</h2></div><div>These are the Matlab files needed to perform this analysis. Please feel free to use them as you will, but please give me credit. A link to my blog is always appreciated.</div><div>Main Stability Analysis Code: <a href="https://drive.google.com/file/d/0B_OtEk0gkL6XS1N6NU9XNks2OFU/view?usp=sharing">cs_orr_sommerfeld_stability.m</a><br />Fourth Order Runge-Kutta Routine: <a href="https://drive.google.com/file/d/0B_OtEk0gkL6XYkthSjRSX1NFdFE/view?usp=sharing">cs_rk4.m</a></div><div class="separator" style="clear: both; text-align: center;"></div><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/wWWmzFj9O0Y" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com10http://www.cdsimpson.net/2015/04/temporal-and-spatial-stability-analysis.htmltag:blogger.com,1999:blog-1308139960516080412.post-20861992575354901922015-04-03T22:19:00.000-05:002015-04-30T20:13:03.429-05:00Derivation and Nondimensionalization of the Orr-Sommerfeld EquationThe Orr-Sommerfeld equation is a famous equation that can give some insight into the stability of the velocity profile of a fluid flow. This is one of two parts on the derivation and stability analysis of the Orr-Sommerfeld equation. In this post, I derive the Orr-Sommerfeld equation starting from the 2-D Navier-Stokes equations. I then show how it can be nondimensionalized. It may look like a lot of math at first glance, but it is all relatively simple.<br /><br /><h2>Derivation</h2><div>The 2-D Navier-Stokes equations are given as follows:<br />$$\begin{align}<br />\nabla\vec{V}&=0\\<br />\rho\frac{D\vec{V}}{Dt}&=-\nabla p+\mu\nabla^{2}\vec{V}<br />\end{align}$$</div><div>Letting \(V_{x}=U+u'\), \(V_{y}=V+v'\), and \(p=P+p'\) and performing a small disturbance analysis gives the small perturbation version of the Navier-Stokes Equations<br />$$\begin{align}\frac{\partial u'}{\partial x}+\frac{\partial v'}{\partial y}&=0\\<br />\frac{\partial u'}{\partial t}+U\frac{\partial u'}{\partial x}+V\frac{\partial u'}{\partial y}+u'\frac{\partial U}{\partial x}+v'\frac{\partial U}{\partial y}&=-\frac{1}{\rho}\frac{\partial p'}{\partial x}+\frac{\mu}{\rho}\left(\frac{\partial^{2}u'}{\partial x^{2}}+\frac{\partial^{2}u'}{\partial y^{2}}\right)\\<br />\frac{\partial v'}{\partial t}+U\frac{\partial v'}{\partial x}+V\frac{\partial v'}{\partial y}+u'\frac{\partial V}{\partial x}+v'\frac{\partial V}{\partial y}&=-\frac{1}{\rho}\frac{\partial p'}{\partial y}+\frac{\mu}{\rho}\left(\frac{\partial^{2}v'}{\partial x^{2}}+\frac{\partial^{2}v'}{\partial y^{2}}\right)<br />\end{align}\tag{1}$$</div><div>Assuming parallel flow, where \(U\approx U(y)\) and \(V\approx0\), we can simplify this to the the following form of the Navier-Stokes equations.</div><div>$$\begin{align}\frac{\partial u'}{\partial x}+\frac{\partial v'}{\partial y}&=0\\\frac{\partial u'}{\partial t}+U\frac{\partial u'}{\partial x}+v'\frac{\partial U}{\partial y}&=-\frac{1}{\rho}\frac{\partial p'}{\partial x}+\frac{\mu}{\rho}\left(\frac{\partial^{2}u'}{\partial x^{2}}+\frac{\partial^{2}u'}{\partial y^{2}}\right)\\\frac{\partial v'}{\partial t}+U\frac{\partial v'}{\partial x}&=-\frac{1}{\rho}\frac{\partial p'}{\partial y}+\frac{\mu}{\rho}\left(\frac{\partial^{2}v'}{\partial x^{2}}+\frac{\partial^{2}v'}{\partial y^{2}}\right)\end{align}\label{simp_NS}\tag{2}$$</div><div>In this analysis, disturbances are assumed to be Tollmien-Schlichting waves, with the general form as follows.</div><div>$$\begin{align}\psi&=\phi(y)e^{i(\alpha x-\omega t)}\\u'&=\frac{\partial\psi}{\partial y}=\frac{\partial\phi}{\partial y}e^{i(\alpha x-\omega t)}\\v'&=-\frac{\partial\psi}{\partial x}=-i\alpha\phi e^{i(\alpha x-\omega t)}\end{align}\tag{3}$$</div><div>The temporal and spatial derivatives are then calculated as follows.</div><div>$$\begin{align}\frac{\partial u'}{\partial t}&=-i\omega\frac{\partial\phi}{\partial y}e^{i(\alpha x-\omega t)}\\\frac{\partial u'}{\partial x}&=i\alpha\frac{\partial\phi}{\partial y}e^{i(\alpha x-\omega t)}\\\frac{\partial^{2}u'}{\partial x^{2}}&=-\alpha^{2}\frac{\partial\phi}{\partial y}e^{i(\alpha x-\omega t)}\\\frac{\partial u'}{\partial y}&=\frac{\partial^{2}\phi}{\partial y^{2}}e^{i(\alpha x-\omega t)}\\\frac{\partial^{2}u'}{\partial y^{2}}&=\frac{\partial^{3}\phi}{\partial y^{3}}e^{i(\alpha x-\omega t)}\\\frac{\partial v'}{\partial t}&=-\alpha\omega\phi e^{i(\alpha x-\omega t)}\\\frac{\partial v'}{\partial x}&=\alpha^{2}\phi e^{i(\alpha x-\omega t)}\\\frac{\partial^{2}v'}{\partial x^{2}}&=i\alpha^{3}\phi e^{i(\alpha x-\omega t)}\\\frac{\partial v'}{\partial y}&=-i\alpha\frac{\partial\phi}{\partial y}e^{i(\alpha x-\omega t)}\\\frac{\partial^{2}v'}{\partial y^{2}}&=-i\alpha\frac{\partial^{2}\phi}{\partial y^{2}}e^{i(\alpha x-\omega t)}\end{align}\tag{4}$$</div><div>We can then substitute each of these derivatives into Equation \(\ref{simp_NS}\) and we get the following relations.</div><div>$$\begin{align}e^{i(\alpha x-\omega t)}\left[i\alpha\frac{\partial\phi}{\partial y}-i\alpha\frac{\partial\phi}{\partial y}\right]&=0\\-\rho e^{i(\alpha x-\omega t)}\left[-i\omega\frac{\partial\phi}{\partial y}+i\alpha U\frac{\partial\phi}{\partial y}+-i\alpha\phi\frac{\partial U}{\partial y}-\frac{\mu}{\rho}\left(-\alpha^{2}\frac{\partial\phi}{\partial y}+\frac{\partial^{3}\phi}{\partial y^{3}}\right)\right]&=\frac{\partial p'}{\partial x}\\-\rho e^{i(\alpha x-\omega t)}\left[-\alpha\omega\phi+U\alpha^{2}\phi-\frac{\mu}{\rho}\left(i\alpha^{3}\phi-i\alpha\frac{\partial^{2}\phi}{\partial y^{2}}\right)\right]&=\frac{\partial p'}{\partial y}\end{align}$$</div><div>To eliminate the pressure fluctuation term, differentiate the x- and y-momentum equations by \(y\) and \(x\), respectively.<br />$$\begin{align}\frac{1}{-\rho e^{i(\alpha x-\omega t)}}\frac{\partial^{2}p'}{\partial x\partial y}&=-i\omega\frac{\partial^{2}\phi}{\partial y^{2}}+i\alpha\frac{\partial U}{\partial y}\frac{\partial\phi}{\partial y}+i\alpha U\frac{\partial^{2}\phi}{\partial y^{2}}-i\alpha\frac{\partial\phi}{\partial y}\frac{\partial U}{\partial y}\\&\quad-i\alpha\phi\frac{\partial^{2}U}{\partial y^{2}}+\frac{\mu}{\rho}\left(\alpha^{2}\frac{\partial^{2}\phi}{\partial y^{2}}-\frac{\partial^{4}\phi}{\partial y^{4}}\right)\\\frac{1}{-i\alpha\rho e^{i(\alpha x-\omega t)}}\frac{\partial^{2}p'}{\partial x\partial y}&=-\alpha\omega\phi+U\alpha^{2}\phi+\frac{\mu}{\rho}\left(-i\alpha^{3}\phi+i\alpha\frac{\partial^{2}\phi}{\partial y^{2}}\right)\end{align}\tag{5}$$<br />Equating the two momentum equations gives<br />$$-i\omega\frac{\partial^{2}\phi}{\partial y^{2}}+i\alpha U\frac{\partial^{2}\phi}{\partial y^{2}}-i\alpha\phi\frac{\partial^{2}U}{\partial y^{2}}+\frac{\mu}{\rho}\left(2\alpha^{2}\frac{\partial^{2}\phi}{\partial y^{2}}-\frac{\partial^{4}\phi}{\partial y^{4}}-\alpha^{4}\phi\right)+i\alpha^{2}\omega\phi-iU\alpha^{3}\phi=0$$<br />This simplifies to the Orr-Sommerfeld Equation.<br />$$\left(U-\frac{\omega}{\alpha}\right)\left(\frac{\partial^{2}\phi}{\partial y^{2}}-\alpha^{2}\phi\right)-\phi\frac{\partial^{2}U}{\partial y^{2}}+\frac{i\nu}{\alpha}\left(\frac{\partial^{4}\phi}{\partial y^{4}}-2\alpha^{2}\frac{\partial^{2}\phi}{\partial y^{2}}+\alpha^{4}\phi\right)=0\label{orrsommerfeld}\tag{6}$$<br /><h2>Nondimensionalization</h2></div><div>The Orr-Sommerfeld equation is nondimensionalized using the following nondimensional parameters,</div><div>$$\bar{U}=\frac{U}{U_{\infty}}\quad\xi=\frac{y}{\delta}\quad\bar{\phi}=\frac{\phi}{U_{\infty}\delta}\quad\bar{c}=\frac{c}{U_{\infty}}\quad\bar{\alpha}=\alpha\delta\quad Re_{\delta}=\frac{U_{\infty}\delta}{\nu}\tag{7}$$</div><div>where \(c=\frac{\omega}{\alpha}\) and \(\delta\) is the boundary layer thickness. Substituting these into Equation \(\ref{orrsommerfeld}\) gives</div><div>$$\begin{align}0&=\left(\bar{U}U_{\infty}-\bar{c}U_{\infty}\right)\left(\frac{1}{\delta^{2}}\frac{\partial^{2}}{\partial\xi^{2}}\left(\bar{\phi}U_{\infty}\delta\right)-\left(\frac{\bar{\alpha}}{\delta}\right)^{2}\left(\bar{\phi}U_{\infty}\delta\right)\right)-\frac{\bar{\phi}U_{\infty}\delta}{\delta^{2}}\frac{\partial^{2}}{\partial\xi^{2}}\left(\bar{U}U_{\infty}\right)\\&\quad+\frac{i\nu\delta}{\bar{\alpha}}\left(\frac{1}{\delta^{4}}\frac{\partial^{4}}{\partial\xi^{4}}\left(\bar{\phi}U_{\infty}\delta\right)-\frac{2\bar{\alpha}^{2}}{\delta^{4}}\frac{\partial^{2}}{\partial\xi^{2}}\left(\bar{\phi}U_{\infty}\delta\right)+\frac{\bar{\alpha}^{4}}{\delta^{4}}\left(\bar{\phi}U_{\infty}\delta\right)\right)\end{align}\tag{8}$$</div><div>Applying the chain rule for each partial derivative gives</div><div>$$\begin{align}0&=U_{\infty}\left(\bar{U}-\bar{c}\right)\left(\frac{U_{\infty}}{\delta}\frac{\partial^{2}\bar{\phi}}{\partial\xi^{2}}-\frac{\bar{\alpha}^{2}U_{\infty}}{\delta}\bar{\phi}\right)-\frac{U_{\infty}^{2}}{\delta}\frac{\partial^{2}\bar{U}}{\partial\xi^{2}}\bar{\phi}\\&\quad+\frac{i\nu\delta}{\bar{\alpha}}\left(\frac{U_{\infty}}{\delta^{3}}\frac{\partial^{4}\bar{\phi}}{\partial\xi^{4}}-\frac{2\bar{\alpha}^{2}U_{\infty}}{\delta^{3}}\frac{\partial^{2}\bar{\phi}}{\partial\xi^{2}}+\frac{\bar{\alpha}^{4}U_{\infty}}{\delta^{3}}\bar{\phi}\right)\end{align}\tag{9}$$</div><div><div>Finally, factor out \(\frac{U_{\infty}^{2}}{\delta}\) and substitute for \(Re_{\delta}\) to get</div></div><div>$$\left(\bar{U}-\bar{c}\right)\left(\frac{\partial^{2}\bar{\phi}}{\partial\xi^{2}}-\bar{\alpha}\bar{\phi}\right)-\frac{\partial^{2}\bar{U}}{\partial\xi^{2}}\bar{\phi}+\frac{i}{\bar{\alpha}Re_{\delta}}\left(\frac{\partial^{4}\bar{\phi}}{\partial\xi^{4}}-2\bar{\alpha}^{2}\frac{\partial^{2}\bar{\phi}}{\partial\xi^{2}}+\bar{\alpha}^{4}\bar{\phi}\right)=0\tag{10}$$</div><div><div>For convenience, derivatives with respect to the station coordinate \(\xi\) are hereafter denoted with prime notation. This gives the final nondimensional form of the Orr-Sommerfeld equation:</div></div><div>$$\left(\bar{U}-\bar{c}\right)\left(\bar{\phi}''-\bar{\alpha}\bar{\phi}\right)-\bar{U}''\bar{\phi}+\frac{i}{\bar{\alpha}Re_{\delta}}\left(\bar{\phi}''''-2\bar{\alpha}^{2}\bar{\phi}''+\bar{\alpha}^{4}\bar{\phi}\right)=0\tag{11}$$</div><div><br /></div><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/0GBOT693lCw" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com7http://www.cdsimpson.net/2015/04/derivation-of-orr-sommerfeld-equation.htmltag:blogger.com,1999:blog-1308139960516080412.post-72827828300182770992015-02-05T17:33:00.000-06:002015-02-06T10:06:21.123-06:00Delaunay TriangulationGeneration of unstructured grids occurs frequently when solving numerical problems in engineering. While structured grids typically use quadrilaterals (and hexahedra in 3D) which are very simple to generate, unstructured grids use triangles (and tetrahedra in 3D), which are much more difficult to generate. Fortunately, there are many programs available to help generate these triangulations.<br /><br />It is important to know how these triangulations are generated. While not the only triangulation method, one method that is almost universally included is Delaunay triangulation. In this post, I explain what Delaunay triangulation is, why it is so popular, and the process that is followed by many programs.<br /><br /><h2>What is a Delaunay triangulation?</h2>By definition, a Delaunay triangulation is a triangulation in which no vertex lies within the circumcircle of any triangle. The circumcircle of a triangle is simply the circle that passes through the three vertices of that triangle. Here is one very simple example. Each domain in the image below contains four points. There are exactly two triangulations for this domain. On the left, we see that one corner of each triangle is inside the circumcircle of the other. On the right, neither circumcircle encompasses the other triangle's opposing vertex. The case on the right is a Delaunay triangulation.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-5Cox0Eq6Odc/VNG6Uj4hDtI/AAAAAAAAALE/KnAX1pr1xeQ/s1600/delaunay_triangulation.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://4.bp.blogspot.com/-5Cox0Eq6Odc/VNG6Uj4hDtI/AAAAAAAAALE/KnAX1pr1xeQ/s1600/delaunay_triangulation.png" /></a></div><br />Here you see the basic concept behind an implementation. First generate a triangulation, then check to see if any vertices fall within the cirumcircle of another triangle. If they do, simply swap the diagonal separating the two triangles. This will always give a Delaunay triangulation.<br /><br /><h2>Why is it so popular?</h2>Delaunay triangulation is fairly simple conceptually, but why is it so popular? The primary reason for its popularity is that the resulting mesh is inherently good quality. For a two-dimensional Delaunay triangulation, it can be shown that the minimum interior angle of each triangle is maximized, and that the maximum interior angle is minimized. The resulting triangles are as equiangular as possible. Another benefit is that the triangulation is independent of the order in which nodes are placed. It is also a relatively simple triangulation method to implement for convex hulls (think of a convex hull as a domain with no "dents" on the boundary). Implementation for non-convex hulls are a bit more challenging, but not impossible.<br /><br /><h2>The Process</h2>The algorithm presented here is as described by S. W. Sloan in Ref. [1]. This routine generates a Delaunay triangulation for a set of predetermined coordinates.<br /><br />First, generate a triangle that encloses the entire domain to be triangulated. This triangle is called the supertriangle. The size and shape of the supertriangle is irrelevant, as long as all points in the domain are contained within.<br /><br />The next portion of the algorithm is an iterative process. For each point in the domain, add the point to the triangulation by subdividing the triangle containing that point. You will then have a new node surrounded by three triangles, seen in green in the image below, and across the far edge of each new triangle is an opposing node, shown in red.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-F5AcoI2OGyQ/VNPmK3DIunI/AAAAAAAAASg/dbpJ6QUki5M/s1600/opposing_nodes.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-F5AcoI2OGyQ/VNPmK3DIunI/AAAAAAAAASg/dbpJ6QUki5M/s1600/opposing_nodes.png" /></a></div><br />For each new triangle created, if the opposing node for that triangle is not part of the supertriangle, check whether the opposing node is within the circumcircle of the new triangle. If it is, swap the diagonal separating the newly inserted node and the opposing node. In the example, two diagonals need to be swapped.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-jK6-lT15_rQ/VNPpSj8vnNI/AAAAAAAAASs/q7lpCfkOs64/s1600/cctest_edge_swap.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-jK6-lT15_rQ/VNPpSj8vnNI/AAAAAAAAASs/q7lpCfkOs64/s1600/cctest_edge_swap.png" /></a></div><br />Swapping the diagonal between two triangles gives two new triangles. Again, a circumcircle test is required for the new opposing nodes, and the process continues. until the circumcircle test is satisfied for all triangles connected to the recently inserted node.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-813w7fYIqBQ/VNPrbLg-9OI/AAAAAAAAAS4/WTtUumZFZfQ/s1600/neighbortest.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://4.bp.blogspot.com/-813w7fYIqBQ/VNPrbLg-9OI/AAAAAAAAAS4/WTtUumZFZfQ/s1600/neighbortest.png" /></a></div><br />The entire process then repeats until all nodes in the domain have been inserted and triangulated. The final step is to remove any triangles connected to the supertriangle vertices, leaving only the domain of interest. Below is a simple example to demonstrate the process one step at a time.<br /><br /><iframe allowfullscreen="true" frameborder="0" height="389" mozallowfullscreen="true" src="https://docs.google.com/presentation/d/1yuAyOCqIwSlo_CxXpsBh8hSm7pEKDR_au0QC1coHYDA/embed?start=false&loop=false&delayms=2000" webkitallowfullscreen="true" width="528"></iframe> <br /><br /><h2>Extension into Three Dimensions</h2>It is not difficult to see that Delaunay triangulation is possible in three dimensions. The obvious difference for a tetrahedral mesh is that a circumsphere is used instead of a circumcircle. One difference that is not so obvious is that instead of swapping the edge dividing two triangles, you need to swap a dividing face. However, when swapping a face, there are two alternative positions, so some care needs to be taken when choosing which option to use.<br /><h2>Potential Uses</h2>Delaunay triangulation, or any triangulation scheme for that matter, is great for connecting a known set of data points. I have used this in conjunction with <a href="http://www.cdsimpson.net/2014/10/barycentric-coordinates.html">barycentric interpolation</a> to create a program that quickly interpolates to find values between known data points. It can also be used to generate a mesh for finite element and finite volume programs. Because the nodes can be inserted in an arbitrary order, it could even be used as a strategy to adapt a mesh in real-time.<br /><div><br />If this post was helpful, let me know in the comments. Ask questions if I have left something unclear, and I'll try to elaborate on the subject.</div><br /><h2>References</h2><ol><li>S. W. Sloan, "A fast algorithm for constructing Delaunay triangulations in the plane", Adv. Eng. Software, Vol. 9, No. 1, pp. 34-55 1987</li></ol><br /><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/elGqTx3yy1M" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com3http://www.cdsimpson.net/2015/02/delaunay-triangulation.htmltag:blogger.com,1999:blog-1308139960516080412.post-45858901253809964912014-10-20T17:25:00.000-05:002015-04-30T20:12:10.186-05:00Barycentric CoordinatesIn the realm of three-dimensional CFD simulations, we tend to use unstructured meshes based on triangles and tetrahedra. If you want to write a program to modify these meshes or even interpolate values, barycentric coordinates can make life much easier.<br /><br />For simplicity, I'll start with the two-dimensional case of a triangle, but this can be extended easily into three dimensions. For a triangle, the barycentric coordinate system is comprised of three dependent basis vectors. Each basis vector extends from the midpoint of one side of the triangle to the opposite corner as shown below.<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-bUBB7qGCmoc/VEVWRj8-xmI/AAAAAAAAAHE/aR2zPQV55lE/s1600/barycentric.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-bUBB7qGCmoc/VEVWRj8-xmI/AAAAAAAAAHE/aR2zPQV55lE/s1600/barycentric.png" /></a></div><br />Notice that each corner of the triangle is now represented by a coordinate triple. In fact, any point \((x,y)\) can be represented by a coordinate triple \((\xi_1,\xi_2,\xi_3)\). We also notice that any point outside the triangle will have at least one negative coordinate. This makes barycentric coordinates extremely useful when determining whether a point is inside a triangle. In addition, the coordinates of a point must always add up to 1.<br /><br />Barycentric coordinates of a point are quite simple to calculate. If we take a point inside a triangle and connect that point to each vertex of the triangle, the result is three subtriangles with areas \(A_1, A_2,\) and \(A_3\) as shown below. If the total area of the triangle is \(A\), then the barycentric coordinates of point \(P\) are simply \(\left(\frac{A_1}{A},\frac{A_2}{A},\frac{A_3}{A}\right)\). <br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://4.bp.blogspot.com/-JcWH4bvsOHQ/VEVt5yF0hFI/AAAAAAAAAHk/wfAajOe7Ids/s1600/barycentric_subtriangles.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://4.bp.blogspot.com/-JcWH4bvsOHQ/VEVt5yF0hFI/AAAAAAAAAHk/wfAajOe7Ids/s1600/barycentric_subtriangles.png" /></a></div><div class="separator" style="clear: both; text-align: center;"></div><br />The barycentric coordinates are simply a ratio of the area of each subtriangle to the total area of the triangle. Now you can see why the coordinates must add up to one. The overall area of the triangle is the sum of the areas of the three subtriangles. "But wait," you may say, "this won't work for points outside the triangle." Actually, it does work. You may recall that area is a vector quantity, which means there is a direction associated with it. This means that a point outside the triangle will give a negative area for one or more triangles. Let's take a look at this in more detail.<br /><br />We will define some vectors to help us out. I'll use the notation \(\vec{\mathbf{V}}_{ab}\) to represent a vector from point \(a\) to point \(b\) in the figure above, and similarly for other vectors. Now, let us define some vectors normal to this triangle:\begin{aligned}<br />\vec{\mathbf{n}}=\vec{\mathbf{V}}_{ab} \times \vec{\mathbf{V}}_{ac} \\<br />\vec{\mathbf{n}}_1=\vec{\mathbf{V}}_{bc} \times \vec{\mathbf{V}}_{bp} \\<br />\vec{\mathbf{n}}_2=\vec{\mathbf{V}}_{ca} \times \vec{\mathbf{V}}_{cp} \\<br />\vec{\mathbf{n}}_3=\vec{\mathbf{V}}_{ab} \times \vec{\mathbf{V}}_{ap} \end{aligned}<br />Now, you may recall that the area of a triangle can be calculated as half of the cross product of the vectors forming two sides of the triangle. We can thus write the areas as follows:<br />\begin{aligned}<br />A=\frac{1}{2}|\vec{\mathbf{n}}| \\<br />A_1=\frac{1}{2}|\vec{\mathbf{n}}_1| \\<br />A_2=\frac{1}{2}|\vec{\mathbf{n}}_2| \\<br />A_3=\frac{1}{2}|\vec{\mathbf{n}}_3| \\ <br />\end{aligned} <br />These always give positive values for area, but we can check whether the area is really negative by comparing the normal vector of the subtriangle with the normal vector of the overall triangle. We do this with a dot product. If the normal vectors are in the same direction, the result will be positive. Likewise, if the normal vectors are in opposite directions, the dot product will be negative. The way we have defined our normal vectors ensures the correct orientation. For example,<br />\begin{equation}\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_1}{|\vec{\mathbf{n}}||\vec{\mathbf{n}}_1|}=%<br />\begin{cases}<br />1 & \text{if P is inside the triangle} \\<br /> -1 & \text{if P is outside the triangle, across edge b-c} \end{cases}<br />\end{equation}<br />It follows that the barycentric coordinates can be calculated by multiplying the area ratio by the corresponding dot-product term.<br />\begin{aligned}<br />\xi_1=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_1}{|\vec{\mathbf{n}}|^2} \\<br />\xi_2=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_2}{|\vec{\mathbf{n}}|^2} \\<br />\xi_3=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_3}{|\vec{\mathbf{n}}|^2} <br />\end{aligned}<br />This works with a triangle in two or three dimensions. We can do something very similar for a tetrahedron. If the barycentric coordinates of a triangle are based on the area of subtriangles, you may have guessed that the barycentric coordinates of a tetrahedron are based on the volume of subtetrahedra, and the coordinates will have four components. The tetrahedron below, formed by vertices \(a,b,c,d\), can be subdivided by inserting a point \(p\).<br /><br /><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-hoCqxpGKjiI/VEWGST2189I/AAAAAAAAAIA/Q0lectt3OW4/s1600/barycentric_tetrahedron.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="http://2.bp.blogspot.com/-hoCqxpGKjiI/VEWGST2189I/AAAAAAAAAIA/Q0lectt3OW4/s1600/barycentric_tetrahedron.png" /></a></div><div class="separator" style="clear: both; text-align: center;"><a href="http://2.bp.blogspot.com/-Bd98NDj959Q/VEWCjSPE-6I/AAAAAAAAAH0/Fxk7zgJ2B7w/s1600/barycentric_tetrahedron.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><br /></a></div>The volume of a tetrahedron can be calculated with a scalar triple product. Let \(V\) be the total volume of the tetrahedron, \(V_a\) be the volume of the subtetrahedron opposite vertex \(a\), and similarly for the others. Then,<br />\begin{aligned}<br />V_a=\frac{1}{6}(\vec{\mathbf{V}}_{bp}\cdot(\vec{\mathbf{V}}_{bd}\times \vec{\mathbf{V}}_{bc})) \\<br />V_b=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ac}\times \vec{\mathbf{V}}_{ad})) \\<br />V_c=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ad}\times \vec{\mathbf{V}}_{ab})) \\ <br />V_d=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ab}\times \vec{\mathbf{V}}_{ac})) \\ <br />V=\frac{1}{6}(\vec{\mathbf{V}}_{ab}\cdot(\vec{\mathbf{V}}_{ac}\times \vec{\mathbf{V}}_{ad})) <br />\end{aligned}<br />The barycentric coordinates can then be calculated as<br />\begin{aligned}<br />\xi_1=\frac{V_a}{|V|} \\<br />\xi_2=\frac{V_b}{|V|} \\<br />\xi_3=\frac{V_c}{|V|} \\<br />\xi_4=\frac{V_d}{|V|} <br />\end{aligned}<br /><br />Both of these algorithms are easily implemented in a computer program, and I am constantly finding uses. Here are a few examples of how you can use barycentric coordinates in your programs.<br /><h2>Barycentric Interpolation</h2>Especially if you work with CFD or finite element codes, you will run into triangular and tetrahedron-based meshes. Sometimes, we want to know the value of some property between nodes on our mesh. We can do this very quickly with barycentric interpolation. If the value at node \(i\) on your element is \(Q_i\), then the value at a point inside that element is given by<br />\[Q=\sum_{i=1}^n Q_i\xi_i\]<br /><h2>Locating the Element That Contains a Point</h2>As I mentioned earlier, if any of the barycentric coordinate components are negative, then the point lies outside of your element. Even cooler is the fact that the point will lie on the side of the element with the negative component. This enables us to create a semi-intelligent search algorithm. We can take an initial guess at which element the point is in. If it's not the right element, we simply check the element across the side with the negative barycentric component. If we continue this, we will eventually find an element that gives us all positive barycentric components. That element contains the point.<br /><br />There are many other uses for barycentric coordinates. If you have used barycentric coordinates for something interesting, let us know in the comments.<br /><h2>Source Code</h2><div>I have written a simple c library that you can use as you desire. Feel free to use, modify, or distribute it as you see fit. All I ask is you give me credit.</div><div><a href="https://drive.google.com/file/d/0B_OtEk0gkL6XSkVSMi0tendDaG8/view?usp=sharing">bary.c</a></div><div><a href="https://drive.google.com/file/d/0B_OtEk0gkL6XaEpzR3FSbDl6S1k/view?usp=sharing">bary.h</a></div><img src="http://feeds.feedburner.com/~r/ChristopherSimpson/~4/JDYkfeagCfs" height="1" width="1" alt=""/>Christopherhttp://www.blogger.com/profile/14768438682715700162noreply@blogger.com0http://www.cdsimpson.net/2014/10/barycentric-coordinates.html