SEAMRotor.py 11.6 KB
Newer Older
Frederik Zahle's avatar
Frederik Zahle committed
1
__author__ = 's127504'
Frederik Zahle's avatar
Frederik Zahle committed
2

Frederik Zahle's avatar
Frederik Zahle committed
3
4
5
6
7
8
from decimal import *

import matplotlib.pyplot as plt
from math import pow
import numpy as np
from scipy.optimize import minimize
Frederik Zahle's avatar
Frederik Zahle committed
9

Frederik Zahle's avatar
Frederik Zahle committed
10
from openmdao.main.api import Component
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
11
from openmdao.lib.datatypes.api import Float, Array, Bool, Int
Frederik Zahle's avatar
Frederik Zahle committed
12

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
13
class SEAMBladeStructure(Component):
Frederik Zahle's avatar
Frederik Zahle committed
14

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
15
    Nsections = Int(iotype='in', desc='number of sections')
Frederik Zahle's avatar
Frederik Zahle committed
16
17
18
19
20
    Neq = Float(1.e7, iotype='in', desc='')

    WohlerExpFlap = Float(iotype='in', desc='')
    PMtarget = Float(iotype='in', desc='')

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
21
    rotor_diameter = Float(iotype='in', units='m', desc='') #[m]
Frederik Zahle's avatar
Frederik Zahle committed
22
23
24
25
#    RootChord = Float(iotype='in', units='m', desc='') #[m]
#    MaxChord = Float(iotype='in', units='m', desc='') #[m]
    MaxChordrR = Float(iotype='in', units='m', desc='') #[m]

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
26
27
    overallMaxFlap = Float(iotype='in', desc='')
    overallMaxEdge = Float(iotype='in', desc='')
Frederik Zahle's avatar
Frederik Zahle committed
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
    TIF_FLext = Float(iotype='in', desc='') # Tech Impr Factor _ flap extreme
    TIF_EDext = Float(iotype='in', desc='')

    FlapLEQ = Float(iotype='in', desc='')
    EdgeLEQ = Float(iotype='in', desc='')
    TIF_FLfat = Float(iotype='in', desc='')

    sc_frac_flap = Float(iotype='in', desc='') # sparcap fraction of chord flap
    sc_frac_edge = Float(iotype='in', desc='') # sparcap fraction of thickness edge

    SF_blade = Float(iotype='in', desc='') #[factor]
    Slim_ext_blade = Float(iotype='in', units='MPa', desc='')
    Slim_fat_blade = Float(iotype='in', units='MPa', desc='')

    AddWeightFactorBlade = Float(iotype='in', desc='') # Additional weight factor for blade shell
    BladeDens = Float(iotype='in', units='kg/m**3', desc='density of blades') # [kg / m^3]
    BladeCostPerMass = Float(iotype='in', desc='') #[e/kg]
    HubCostPerMass = Float(iotype='in', desc='') #[e/kg]
    SpinnerCostPerMass = Float(iotype='in', desc='') #[e/kg]

    BladeWeight = Float(iotype = 'out', units = 'kg', desc = 'BladeMass' )
49
    RootChord = Float(iotype = 'out', units = 'm', desc = 'blade root chord') # 07/09/2015 added for HubSE model
Frederik Zahle's avatar
Frederik Zahle committed
50
51

    def execute(self):
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
52
53
54

        volumen = 0.
        r = np.linspace(0, self.rotor_diameter/2., self.Nsections)
Frederik Zahle's avatar
Frederik Zahle committed
55
56
        C = np.zeros(r.shape)
        thick = np.zeros(r.shape)
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
57
        RootChord = (self.rotor_diameter/2.)/25.
Frederik Zahle's avatar
Frederik Zahle committed
58
59
60
        if RootChord < 1:
            RootChord = 1

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
61
        MaxChord = 0.001*(self.rotor_diameter/2.)**2-0.0354*(self.rotor_diameter/2.)+1.6635  #Empirical
Frederik Zahle's avatar
Frederik Zahle committed
62
63
64
        if MaxChord < 1.5:
            RootChord = 1.5

65
66
        self.RootChord = RootChord

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
67
        for i in range(self.Nsections):
Frederik Zahle's avatar
Frederik Zahle committed
68
69

            # Calculating the chord and thickness
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
70
71
72
73
74
75
76
77
            if r[i] > (self.MaxChordrR*self.rotor_diameter/2.):
                C[i] = MaxChord - (MaxChord)/((self.rotor_diameter/2.) - \
                                   self.MaxChordrR*self.rotor_diameter/2.)*(r[i] - \
                                  (self.MaxChordrR*self.rotor_diameter/2.))
                thick[i] = 0.4*MaxChord - \
                          (0.4*MaxChord-0.05*MaxChord)/(self.rotor_diameter/2. - \
                          self.MaxChordrR*self.rotor_diameter/2.)*(r[i] - \
                          self.MaxChordrR*self.rotor_diameter/2.)
Frederik Zahle's avatar
Frederik Zahle committed
78
79
80
81
82
83

                if thick[i]<0.001:
                    thick[i] = 0.001
                if C[i]<0.001:
                    C[i] = 0.001
            else:
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
84
85
86
                C[i] = MaxChord-(MaxChord)/((self.rotor_diameter/2.) - \
                       self.MaxChordrR*(self.rotor_diameter/2.)) * \
                       (self.MaxChordrR*(self.rotor_diameter/2.)-r[i])
Frederik Zahle's avatar
Frederik Zahle committed
87
88
89
90
91
92
93
94
95
96
97
                thick[i] = RootChord

                if thick[i]<0.001:
                    thick[i] = 0.001
                if C[i]<0.001:
                    C[i] = 0.001

        self.thick = thick
        self.chord = C

        # Calculating load flapwise extreme
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
98
99
100
        Mext_flap = (self.overallMaxFlap-1.75*self.overallMaxFlap*r/(self.rotor_diameter/2.) +\
                     0.75*self.overallMaxFlap*r/(self.rotor_diameter/2.) * r / \
                     (self.rotor_diameter/2.))*self.TIF_FLext
Frederik Zahle's avatar
Frederik Zahle committed
101
102

        # Calculating load edgewise extreme
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
103
104
105
        Mext_edge = (self.overallMaxEdge-1.75*self.overallMaxEdge*r/(self.rotor_diameter/2.) + \
                    0.75*self.overallMaxEdge*r/(self.rotor_diameter/2.) * r / \
                    (self.rotor_diameter/2.))*self.TIF_EDext
Frederik Zahle's avatar
Frederik Zahle committed
106
107

        # Calculating load flapwise fatigue
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
108
109
110
        Mfat_flap = (self.FlapLEQ-1.75*self.FlapLEQ*r/(self.rotor_diameter/2.) + \
                    0.75*self.FlapLEQ*r/(self.rotor_diameter/2.)*r / \
                    (self.rotor_diameter/2.))*self.TIF_FLfat
Frederik Zahle's avatar
Frederik Zahle committed
111
112

        # Calculating load edgewise fatigue
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
113
114
115
116
117
118
119
        Mfat_edge = (self.EdgeLEQ-1.75*self.EdgeLEQ*r/(self.rotor_diameter/2.) + \
        0.75*self.EdgeLEQ*r/(self.rotor_diameter/2.)*r/(self.rotor_diameter/2.))*self.TIF_EDext

        self.Mext_flap = Mext_flap
        self.Mext_edge = Mext_edge
        self.Mfat_flap = Mfat_flap
        self.Mfat_edge = Mfat_edge
Frederik Zahle's avatar
Frederik Zahle committed
120
121
122

        # Calculating thickness flapwise extreme
        text_flap = np.zeros(self.Nsections)
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
123
124
125
        for i in range(self.Nsections):
            norm = self.solve_text_flap(0.01, C[i], thick[i], Mext_flap[i], 1.)
            res = minimize(self.solve_text_flap, 0.01, args = (C[i], thick[i], Mext_flap[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
Frederik Zahle's avatar
Frederik Zahle committed
126
            text_flap[i] = res['x']
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
127
128

            if not res['success']: print 'WARNING enter text_flap', i,res
Frederik Zahle's avatar
Frederik Zahle committed
129
130
131

        # Calculating thickness edgewise extreme
        text_edge = np.zeros(self.Nsections)
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
132
133
134
        for i in range(self.Nsections):
            norm = self.solve_text_edge(0.01, C[i], thick[i], Mext_flap[i], 1.)
            res = minimize(self.solve_text_edge, 0.01, args = (C[i], thick[i], Mext_edge[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
Frederik Zahle's avatar
Frederik Zahle committed
135
            text_edge[i] = res['x']
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
136
            if not res['success']: print 'WARNING solve_text_edge', i,res
Frederik Zahle's avatar
Frederik Zahle committed
137
138
139

        # Calculating thickness flapwise fatigue
        tfat_flap = np.zeros(self.Nsections)
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
140
141
142
        for i in range(self.Nsections):
            norm = self.solve_tfat_flap(0.01, C[i], thick[i], Mext_flap[i], 1.)
            res = minimize(self.solve_tfat_flap, 0.01, args = (C[i], thick[i], Mfat_flap[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
Frederik Zahle's avatar
Frederik Zahle committed
143
            tfat_flap[i] = res['x']
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
144
            if not res['success']: print 'WARNING solve_tfat_flap', i, res
Frederik Zahle's avatar
Frederik Zahle committed
145
146
147

        # Calculating thickness edgewise fatigue
        tfat_edge = np.zeros(self.Nsections)
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
148
149
150
        for i in range(self.Nsections):
            norm = self.solve_tfat_edge(0.01, C[i], thick[i], Mext_flap[i], 1.)
            res = minimize(self.solve_tfat_edge, 0.01, args = (C[i], thick[i], Mfat_edge[i], norm), bounds=[(1.e-6, 0.5)], method = 'SLSQP', tol = 1.e-8)
Frederik Zahle's avatar
Frederik Zahle committed
151
            tfat_edge[i] = res['x']
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
152
153
            if not res['success']: print 'WARNING solve_tfat_edge', i,res

Frederik Zahle's avatar
Frederik Zahle committed
154
155
156
157
158
159
160
161
162
        self.text_flap = text_flap
        self.text_edge = text_edge
        self.tfat_flap = tfat_flap
        self.tfat_edge = tfat_edge
        self.r = r

        # Calculating the final thickness for flap and edge direction
        tfinal_flap = np.maximum(text_flap, tfat_flap)
        tfinal_edge = np.maximum(text_edge, tfat_edge)
Frederik Zahle's avatar
Frederik Zahle committed
163
164
        self.tfinal_flap = tfinal_flap
        self.tfinal_edge = tfinal_edge
Frederik Zahle's avatar
Frederik Zahle committed
165
166
167
168
169
170
171
        # Calculating different costs and mass
        for i in range(0,int(self.Nsections)):
            if i>0:
                volumen = volumen+self.AddWeightFactorBlade*(r[i]-r[i-1])*((2*self.sc_frac_flap*C[i]*tfinal_flap[i]
                            +2*self.sc_frac_edge*thick[i]*tfinal_edge[i])+(2*self.sc_frac_flap*C[i-1]*tfinal_flap[i-1]
                            +2*self.sc_frac_edge*thick[i-1]*tfinal_edge[i-1]))/2.

Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
172
        self.BladeWeight = self.BladeDens*volumen
Frederik Zahle's avatar
Frederik Zahle committed
173
174
175

        print 'volumen', volumen
        print 'BladeWeigth', self.BladeWeight
Frederik Zahle's avatar
Frederik Zahle committed
176
        #
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
        # rotor = self.BladeWeight*self.BladeCostPerMass/1e6  # Meuro
        # print 'rotor', rotor

        # # Scaling laws for hub, spinner, pitch
        # HubMass = 0.954*self.BladeWeight+5680.3
        # print 'HubMass', HubMass
        #
        # HubCost = HubMass*self.HubCostPerMass/1e6
        # print 'HubCost', HubCost
        #
        # #PitchCost
        #
        # SpinnerMass = 18.5*self.rotor_diameter-520.5
        # print 'SpinnerMass', SpinnerMass
        #
        # SpinnerCost = SpinnerMass*self.SpinnerCostPerMass/1e6
        # print 'SpinnerCost', SpinnerCost
Frederik Zahle's avatar
Frederik Zahle committed
194
195
196
197
198
199
200
201
        #Rotor = 3*(self.BladeWeight*self.BladeCostPerMass/1e6)+HubCost+PitchCost+SpinnerCost




    # Defining functions in order to solve for the thickness

    # Solving for t in flap direction, extremes
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
202
    def solve_text_flap(self, t, C, thick, Mext_flap, norm):
Frederik Zahle's avatar
Frederik Zahle committed
203
204
205
        Ine = (2./3.)*self.sc_frac_flap*C*t**3-self.sc_frac_flap*C*thick*t**2+self.sc_frac_flap*C*thick**2/2.*t
        W = Ine/(thick/2.)
        sext = self.SF_blade*1e3*Mext_flap/W/1e6
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
206
207
        # print 'text_flap', abs(sext - self.Slim_ext_blade)
        return abs(sext - self.Slim_ext_blade) / norm
Frederik Zahle's avatar
Frederik Zahle committed
208
209

    #Solving for t in edge direction, extremes
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
210
    def solve_text_edge(self, t, C, thick, Mext_edge, norm):
Frederik Zahle's avatar
Frederik Zahle committed
211
212
213
        Ine = (2/3.)*self.sc_frac_edge*thick*t**3-self.sc_frac_edge*thick*C*t**2+self.sc_frac_edge*thick*C**2/2.*t
        W = Ine/(C/2.);
        sext = self.SF_blade*1.e3*Mext_edge/W/1.e6
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
214
        return abs(sext - self.Slim_ext_blade) / norm
Frederik Zahle's avatar
Frederik Zahle committed
215
216

    # Solving for t in flap direction, fatigue
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
217
    def solve_tfat_flap(self, t, C, thick, Mfat_flap, norm):
Frederik Zahle's avatar
Frederik Zahle committed
218
219
220
221
        Ine = (2/3.)*self.sc_frac_flap*C*t**3-self.sc_frac_flap*C*thick*t**2+self.sc_frac_flap*C*thick**2/2.*t
        W = Ine/(thick/2.)
        sfat = np.maximum(1.e-6, self.SF_blade*1.e3*Mfat_flap/ W / 1.e6)
        PM = self.Neq/(pow(10, (self.Slim_fat_blade - self.WohlerExpFlap*np.log10(sfat))))
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
222
        return abs(PM - self.PMtarget) / norm
Frederik Zahle's avatar
Frederik Zahle committed
223
224

    # Solving for t in edge direction, fatigue
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
225
    def solve_tfat_edge(self, t, C, thick, Mfat_edge, norm):
Frederik Zahle's avatar
Frederik Zahle committed
226
227
228
229
        Ine = (2/3.)*self.sc_frac_edge*thick*t**3-self.sc_frac_edge*thick*C*t**2+self.sc_frac_edge*thick*C**2/2.*t
        W = Ine/(C/2.)
        sfat = np.maximum(1.e-6, self.SF_blade*1.e3*Mfat_edge/W/1.e6)
        PM = self.Neq/(pow(10, (self.Slim_fat_blade - self.WohlerExpFlap*np.log10(sfat))))
Frederik Zahle's avatar
cleanup    
Frederik Zahle committed
230
        return abs(PM - self.PMtarget) / norm
Frederik Zahle's avatar
Frederik Zahle committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273

    def plot(self, fig):
        """
        function to generate Bokeh plot for web GUI.

        Also callable from an ipython notebook

        parameters
        ----------
        fig: object
            Bokeh bokeh.plotting.figure object

        returns
        -------
        fig: object
            Bokeh bokeh.plotting.figure object
        """
        try:
            # formatting
            fig.title = 'Power curve'
            fig.xaxis[0].axis_label = 'Wind speed [m/s]'
            fig.yaxis[0].axis_label = 'Power production [kW]'

            # fatigue, ultimate and final thickness line plots
            # fig.line(self.r, self.self.text_flap, line_color='orange',
            #                             line_width=3,
            #                             legend='Extreme flap')
            # fig.line(self.r, self.text_edge, line_color='green',
            #                             line_width=3,
            #                             legend='Extreme edge')
            # fig.line(self.r, self.tfat_flap, line_color='blue',
            #                             line_width=3,
            #                             legend='Fatigue flap')
            fig.line(self.r, self.tfinal_flap, line_color='red',
                                        line_width=3,
                                        legend='tfinal_flap edge')
            fig.line(self.r, self.tfinal_edge, line_color='blue',
                                        line_width=3,
                                        legend='tfinal_edge')
        except:
            pass

        return fig