Kinematics is mathematical description of motion without considering the physical forces needed to perform the movement. Kinematics is used in this project to change the posture of the cyclist, by the principle of inverse kinematics, which is the same principle as for inverse mathematics, you know the answer but you do not know how to calculate it. Consider the robot arm consisting of two revolute joints presented in figure 2.5. Let us assume that the tool-tip position, i.e. the red circle in Figure 2.5, of the robot is a known position (x0,y0). However, we do not know the joint anglesθa andθb. By the principle of inverse kinematics, these angles can be calculated for any tool-tip position.
The tooltip position of the robot is given by
x=l1cos(θ1) +l2cos(θ1+θ2) (2.38) y=l1sin(θ1) +l2sin(θ1+θ2) (2.39)
Figure 2.6: Figure 2.5 with additional definitions needed to find the equation forθ1
which is known as the forward kinematics. Now we start our work on deriving the inverse kinematics, i.e. explicit equations forθ1 andθ2. First, we square and sum Equation 2.38 and 2.39, which gives us x2+y2=l21cos2(θ1) +l22cos2(θ1+θ2) + 2l1l2cos(θ1)cos(θ1+θ2) +l12sin2(θ1) (2.40)
+l22sin2(θ1+θ2) + 2l1l2sin(θ1)sin(θ1+θ2)
=l21+l22+ 2l1l2[cos(θ1)cos(θ1+θ2+ sin(θ1)sin(θ1+θ2)]
Next, consider the trigonometric identities
sin(x±y) = sin(x)cos(y)±cos(x)sin(y) (2.41) cos(x±y) = cos(x)cos(y)±sin(x)sin(y) (2.42) By applying these on Equation 2.40, it can be rewritten as
x2+y2=l21+l22+ 2l1l2[cos(θ1)(cos(θ1)cos(θ2)−sin(θ1)sin(θ2)) (2.43) + sin(θ1)(sin(θ1)cos(θ2)−cos(θ1)sin(θ2))]
=l21+l22+ 2l1l2[cos2(θ1)cos(θ2) + sin2(θ2)cos(θ2)]
=l21+l22+ 2l1l2cos(θ2) Next, we can use Equation 2.44 to expressθ2 explicit
cos(θ2) =x2+y2−l21−l22
Our first step in finding the equation forθ1is to rewrite the forward kinematics equations with the basis of the definitions presented in Figure 2.6.
x=k1cos(θ1)−k2sin(θ1) (2.46)
y=k1sin(θ1) +k2cos(θ1)where (2.47) where
k1=l1+l2cos(θ2) (2.48)
k2=l2sin(θ2) (2.49)
whereatan2 is the four-quadrant inverse tangent of a coordinate. Based on the coordinate, the angle αbetween the coordinate and the positive defined horizontal axis within in the interval−π < α < πis returned. As opposed to the traditionalatan, which returns an angleαin the interval−π/2< α < π/2.
Inserting Equation 2.52 and 2.53 into Equation 2.46 and 2.47 yields
x=rcos(γ)cos(θ1)−rsin(γ)sin(θ1) (2.54)
Finally, applying the atan2 function γ+θ1atan2(y
The software used in this project for adjusting the cyclist posture, which is Blender, uses inverse kinematics to calculate corresponding joint angles between limbs when for example a hand is moved.
The calculation is embedded in the software.
Chapter 3
Simulation of flow over a sphere
As part of the process of developing the computational setup, it was tested on a simpler geometry, for which information about physical experiments is available. Simulation of flow around a sphere challenges the computational setup, as well as being modelling-wise a simple geometry, with several publications on both physical experiments and numerical studies available. In this thesis, the sphere is simulated for Reynolds number equal to 104 and 106. By Equation 2.2, these Reynolds numbers corresponds to velocities of 0.15 and 15 m/s, where the latter is roughly the typical velocity for a professional cyclist by assuming the properties of air to beρ=1 kg/m3andµ=1.5·10−5 Pa·s and L=1 m.
Research has been carried out on the flow around a sphere for almost a 100 years, with physical experiments of Wieselsberger dating back to 1922. Wieselsberger (1922) measured the drag force from Re= 10 through the critical flow regime, i.e. up to Re = 106, and observed the characteristic drop of CD in the critical flow regime. The reported results of the drag coefficient at high Reynolds numbers, together with results from the researchers presented in the following, are plotted in Figure 3.1. Milikan & Klein (1933) measured in the range fromRe= 2·105 to Re= 8·105, that is into the supercritical regime, and observed the characteristic drop ofCd, similar, although far from identical, to Wieselsberger. Achenbach (1972) extended the measurement range, and measured fromRe= 5·104to Re= 6·106, and obtained results quite similar to Milikan. The results obtained Achenbach are widely cited by later work performed in the recent years, for which the trends have shifted from physical experiments to computational experiments using CFD.
Constantinescu, Chapelet & Squires (2003) studied the turbulence modelling applied to flow over a sphere, and presented numerical simulations of the subcritical flow regime, with the aim of comparing the results from three different approaches to model turbulence; unsteady Reynolds-Averaged Navier-Stokes (URANS), detached eddy simulation(DES) and large eddy simulation (LES). They concluded that the URANS predictions, withk−ωturbulence modelling, of the pressure coefficient, skin friction, and (by association) the streamwise drag, were in reasonable agreement with measurements. Constan-tinescu & Squires (2004) performed a numerical investigation of flow over a sphere in the subcritical and supercritical regimes using DES simulations and were able to capture many of the features that characterise the subcritical and supercritical regimes, as revealed by experimental investigations. Jones
& Clarke (2008) simulated the flow past a sphere using the commercial CFD Fluent code in multi-ple flow regimes and confirmed that the capabilities of Fluent to accurately reproduce typical flow structures observed for both time-independent/time-dependent and laminar/turbulent flow regimes.
3.1 Computational setup
The sphere, which has a diameter of 1m, is modelled in a computer-assisted drawing (CAD) program, saved in STL file format, and imported and meshed in OpenFOAM using the two utilities; blockMesh
104 105 106 107
Figure 3.1: A comparison of the literature on the (experimental) drag coefficient of the sphere as a function of the Reynolds number
and snappyHexMesh. blockMesh is used to create what is known as the background mesh in the computational domain, while snappyHexMesh refines the mesh around the sphere according to the user’s specifications.
The grid is generated based on the recommendations presented in Chapter 2.3, withy+set to 1,y is calculated by Equation 2.36, which is
y= 1
Inserting the Re corresponding to each simulation in the latter equation and the constants for air yields Further,nce is calculated by Equation 2.37, which is
nce=ln ycl−ln y lnera
=ln ycl−lny
ln (1.2) (3.4)
ycl is set equal to the uniform grid which is surrounding the sphere, which is treated as a parameter to determine a sufficient grid size. As we will see in the results, a grid consisting of 5.13·105 and 6.55·106 cells are sufficient to reach convergence withRe equal to 104 and 106, respectively. These meshes are constructed by setting the sides of the cubic grid cells to be 2.27·10−2mand 6.25·10−3m, respectively. Thus, by settingycl equal to these values,nce can be calculated for each mesh
Re= 104, nce= ln(2.27·10−2)−ln(1.7·10−3)
ln(1.2) ≈14 (3.5)
Re= 106, nce= ln(6.25·10−3)−ln(2.4·10−5)
ln(1.2) ≈30 (3.6)
An overview of the meshes used withRe = 104 andRe = 106, which had sufficient number of cells, are presented in Figure 3.2 and 3.4, with additional zoom on the viscous mesh in Figure 3.3 and 3.5, respectively.
The fluid flow surrounding the sphere is simulated in steady-state as an incompressible Newtonian turbulent fluid, which by Chapter 2.2 leave us with the obvious solver choice of simpleFoam. The governing equations presented in Chapter 2.1 applies to these simulations together with the with the turbulence model SST k−ω. The convergence criteria are based on the residuals for pressure and velocity in the x-, y- and z-direction, which are to be smaller than 10−5 before the simulation converges. A uniform constant horizontal velocity of 0.15 m/s and 15 m/s and zero gradient for pressure are imposed at the inlet of the fluid domain, which corresponds to the left vertical edge of the meshes presented in Figure 3.2 and 3.4, respectively. At the outlet, i.e. the right vertical edge of the mesh in the two latter figures, a pressure outlet condition with ambient static pressure is imposed. A zero gradient for pressure is imposed on the other sides of the domain, together with slip for velocity.
Finally, for the sphere, a no-slip for velocity- and a zero gradient for pressure boundary condition is imposed. Additionally, the turbulent boundary conditions described in Chapter 2.3 are embedded in the software, except for k and ω at the inlet. Their numerical value are approximated by Equation (2.23) and (2.24), which are
k= 3
Since we are to compare the simulation results with wind tunnel experiments performed in the lit-erature, a lowTi is expected. Hence, it’s set equal to 0.5%. As opposed to engineering flows which typically have aTi between 2 and 5%(Versteeg & Malalasekera, 2007). By inserting the velocity, the correspondingkandω can be calculated for both simulations
k(U∞= 0.15 m/s) = 3
Figure 3.2: Overview of the converging mesh atRe= 104 in the computational domain, its is refined in areas of particular interest
Figure 3.3: Figure 3.2 with additional zoom on the viscous mesh on the sphere
Figure 3.4: Overview of the converging mesh atRe= 106 in the computational domain, its is refined in areas of particular interest, similar to the converging mesh atRe= 104
Figure 3.5: Figure 3.4 with additional zoom on the viscous mesh on the sphere