Skip to content
Snippets Groups Projects
Commit dbed6fd2 authored by Jenni Rinker's avatar Jenni Rinker
Browse files

fixing download links again

parent fd2f06bc
No related branches found
No related tags found
No related merge requests found
Pipeline #25072 passed
%% Cell type:markdown id:fa60cf2d tags:
# Forced response and save to file
%% Cell type:raw id:7b072291 tags:
To do the exercises in this notebook, please `download:`download the wind file <wind_12_ms_TI_0.1.txt>` and `download:`download the results file <res_12_ms_TI_0.1.txt>`
To do the exercises in this notebook, please download the wind file and save it as ``wind_12_ms_TI_0.1.txt``: :download:`download file <wind_12_ms_TI_0.1.txt>`. You will also need the results file from the loading and plotting exercises.
%% Cell type:markdown id:fd77dbf1 tags:
It is now time to add the wind to our Turbie simulations.
Remember from the definition of Turbie that
\begin{equation}
f_{aero}(t) = \frac12\, \rho\, C_T\, A\, (u(t) - \dot{x}_1)\,|u(t) - \dot{x}_1|,
\end{equation}
where $C_T$ is constant from a simulation and calculated from a look-up table using the mean wind speed of a simulation, $U=\overline{u(t)}$. Furthermore, the forcing vector for our 2-DOF system is
\begin{equation}
\overline{F}(t) = \left[ \begin{array}{c} f_{aero}(t) \\ 0\end{array} \right].
\end{equation}
### Exercises for the reader
1. Create a function that takes as input the path to a wind file and returns the time array and the wind array. Try it on the wind file. What is the wind speed at $t=1$ s?
2. Create a function that takes as input a wind array and returns the value for $C_T$ to be used in a simulation. Test it on the wind file.
3. Update your `dydt` function from the homogeneous response to include the aerodynamic forcing as defined above. What is the value of your new `dydt` function when $t=1$ and $y=[1, 2, 3, 4]$?
4. Use your numerical integrator, just like in the homogeneous response, to simulate Turbie's response with the following conditions.
* Time vector spans from 0 to 660 in spaces of 0.01 s.
* Turbie initial conditions are `y = [0, 0, 0, ]`.
5. Load the tower deflection from the results file. Plot the tower deflection in that file versus the tower deflection you just simulated. Do they match?
* Don't worry about the blade deflections matching -- the res file has a different coordinate system.
6. Save the results to a text file in a format similar to the results file we give you.
* **Python hint**: Use `np.savetxt` with `delimiter='\t'` and `fmt='%.3f'` to get a nice file. You should also pass in a header string that explains the columns.
### Answers to the exercises
#### Exercise 1
%% Cell type:code id:305fdb38 tags:
``` python
import numpy as np
def load_wind(path):
"""Load the time and wind speed from a file.
Returns t, u.
"""
t_wind, wind = np.loadtxt(path, skiprows=1).T
return t_wind, wind
# test the function
path = 'wind_12_ms_TI_0.1.txt'
t_test = 1 # time we want the wind speed
t_wind, wind = load_wind(path) # load the wind
idx_test = np.argmin(np.abs(t_wind - t_test)) # index closest to t_test
print(wind[idx_test]) # wind speed at that index
```
%% Cell type:markdown id:047090f9 tags:
#### Exercise 2
%% Cell type:code id:131be2c3 tags:
``` python
def calculate_ct(u, ct_path='CT.txt'):
"""Calculate Ct using lookup table and wind
time series. ct = CT(U).
"""
u_mean = np.mean(u) # calculate mean of wind time series
u_lut, ct_lut = np.loadtxt(ct_path, skiprows=1).T # load lookup table
ct = np.interp(u_mean, u_lut, ct_lut) # interp Ct curve
return ct
print('The Ct for this wind series is: ', calculate_ct(wind))
```
%% Cell type:markdown id:6c826bae tags:
#### Exercise 3
%% Cell type:code id:c7fd4f63 tags:
``` python
def dydt(t, y, M, C, K, ct, rho, area, t_wind, wind):
"""Differential function for Turbie (with forcing)"""
Minv = np.linalg.inv(M) # save inverse
ndof = M.shape[0] # get number of degrees of freedon
# ---- assemble matrix A ----
I, O = np.eye(ndof), np.zeros((ndof, ndof))
A = np.block([[O, I], [-Minv @ K, -Minv @ C]])
# ---- define forcing vector ----
F = np.zeros(ndof) # initialize forcing vector
# calculate aerodynamic force
x1dot = y[2] # third element is blade velocity
u = np.interp(t, t_wind, wind) # interpolate current wind
faero = 0.5 * rho * area * ct * (u - x1dot) * np.abs(u - x1dot)
F[0] = faero
# assemble array B
B = np.zeros(2*ndof) # initialize the array
B[ndof:] = Minv @ F
return A @ y + B
from turbie import get_turbie_system_matrices
# load/calculate necessary parameters
Dr, rho = np.loadtxt('turbie_parameters.txt', comments='%')[-2:]
area = np.pi * (Dr / 2) **2 # rotor area
M, C, K = get_turbie_system_matrices()
ct = calculate_ct(wind)
# test the function
t_test, y_test = 1, [1, 2, 3, 4]
dydt(t_test, y_test, M, C, K, ct, rho, area, t_wind, wind)
```
%% Cell type:markdown id:1af79848 tags:
#### Exercise 4
%% Cell type:code id:5d118383 tags:
``` python
from scipy.integrate import solve_ivp
# define our time vector
t0, tf, dt = 0, 660, 0.01
# inputs to solve ivp
tspan = [t0, tf] # 2-element list of start, stop
y0 = [0, 0, 0, 0] # initial condition
t_eval = np.arange(t0, tf, dt) # times at which we want output
args = (M, C, K, ct, rho, area, t_wind, wind) # extra arguments to dydt besides t, y
# run the numerical solver
res = solve_ivp(dydt, tspan, y0, t_eval=t_eval, args=args)
# extract the output
t, y = res['t'], res['y']
```
%% Cell type:markdown id:2cd5b73e tags:
#### Exercise 5
%% Cell type:code id:aeb17bc5 tags:
``` python
import matplotlib.pyplot as plt
from turbie import load_results
# load values from given results file
path = 'res_12_ms_TI_0.1.txt'
t_wind, _, xb_txt, xt_txt = load_results(path) # load results file (skip wind)
# get out tower deflections
xb_sim = y[0, :] # blade deflection is first element in state vector
xt_sim = y[1, :] # tower deflection is second element in state vector
# plot them
plt.plot(t_wind, xt_txt, label='In file')
plt.plot(t, xt_sim, label='Simulated')
plt.legend()
plt.xlim([0, 200])
```
%% Cell type:markdown id:756cf057 tags:
They match, yay!
%% Cell type:markdown id:c70b6f82 tags:
#### Exercise 6
%% Cell type:code id:a1144c1b tags:
``` python
# create the 2D array we want to save
nt = t_eval.size # number of time steps
ncol = 4 # time, wind speed, xb, xt
out_arr = np.empty((nt, ncol))
out_arr[:, 0] = t_eval
out_arr[:, 1] = np.interp(t_eval, t_wind, wind)
out_arr[:, 2] = xb_sim
out_arr[:, 3] = xt_sim
# save the file
out_path = 'myres.txt'
header = 'Time(s)\tV(m/s)\txb(m)\txt(m)'
np.savetxt(out_path, out_arr, header=header, fmt='%.3f', delimiter='\t', comments='')
```
%% Cell type:markdown id:d75df559 tags:
......
%% Cell type:markdown id:fa60cf2d tags:
# Bring it all together
This is the final step before the final project!
### Exercises for the reader
%% Cell type:raw id:3af72a2b tags:
To do the exercises in this notebook, you will need the wind file from the forced response exercise and the results file from the loading and plotting exercise.
To do the exercises in this notebook, please `download:`download the wind file <wind_12_ms_TI_0.1.txt>`
### Exercises for the reader
%% Cell type:markdown id:d272639b tags:
1. Write a function that:
* Takes as input:
* a path to a wind file,
* time to start simulating,
* time to stop simulating,
* time step, and
* a path to an output file to write,
* Simulates turbie's response to that wind, and
* Saves the simulation to a Turbie result file with the name given in the input.
2. Test this function on the wind file you can download at the top of this page and compare to your results for the forced response.
You can assume that the initial conditions are all zero.
That's it for the exercises! When you're done, you can take a shot at the final project. Have fun!
### Answers to the exercises
%% Cell type:code id:f8b04af3 tags:
``` python
import numpy as np
from scipy.integrate import solve_ivp
from turbie import load_wind, calculate_ct, get_turbie_system_matrices, dydt
def simulate_turbie(wind_path, t0, tf, dt, save_path):
"""Run Turbie on a wind file and save the result
"""
print('Loading values...')
# load/calculate necessary parameters
t_wind, wind = load_wind(wind_path) # wind for aero forcing
Dr, rho = np.loadtxt('turbie_parameters.txt', comments='%')[-2:]
area = np.pi * (Dr / 2) **2 # rotor area
M, C, K = get_turbie_system_matrices()
ct = calculate_ct(wind)
# inputs to solve_ivp
tspan = [t0, tf] # 2-element list of start, stop
y0 = [0, 0, 0, 0] # initial conditions all zero
t_eval = np.arange(t0, tf, dt) # times at which we want output
args = (M, C, K, ct, rho, area, t_wind, wind) # extra arguments to dydt besides t, y
# run the numerical solver
print('Simulating...')
res = solve_ivp(dydt, tspan, y0, t_eval=t_eval, args=args)
# extract the output
t, y = res['t'], res['y']
# create the 2D array we want to save
nt = t_eval.size # number of time steps
ncol = 4 # time, wind speed, xb, xt
out_arr = np.empty((nt, ncol))
out_arr[:, 0] = t_eval
out_arr[:, 1] = np.interp(t_eval, t_wind, wind)
out_arr[:, 2] = y[0, :] # blade deflection
out_arr[:, 3] = y[1, :] # tower deflection
# save the file
header = 'Time(s)\tV(m/s)\txb(m)\txt(m)'
np.savetxt(save_path, out_arr, header=header, fmt='%.3f', delimiter='\t', comments='')
print(f'File saved to {save_path}.')
return
# test the function
wind_path = 'wind_12_ms_TI_0.1.txt'
save_path = 'myres.txt'
t0, tf, dt = 0, 660, 0.01
simulate_turbie(wind_path, t0, tf, dt, save_path)
```
%% Cell type:markdown id:bec79f56 tags:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment