Merge pull request #75 from sandialabs/dev

Update 34m Example, Optionally Calculate Fatigue
sandialabs · Sep 4, 2024 · 926ad39 · 926ad39
2 parents ff1fc4b + 5231c6b
commit 926ad39
Show file tree

Hide file tree

Showing 3 changed files with 67 additions and 44 deletions.
diff --git a/examples/SNL34m/SNL34mVAWTNormalOperation.jl b/examples/SNL34m/SNL34mVAWTNormalOperation.jl
@@ -160,6 +160,27 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng,
     angularOffset = pi/2,
     meshtype = turbineType)
 
+PyPlot.figure()
+for icon = 1:length(mymesh.conn[:,1])
+    idx1 = mymesh.conn[icon,1]
+    idx2 = mymesh.conn[icon,2]
+    PyPlot.plot3D([mymesh.x[idx1],mymesh.x[idx2]],[mymesh.y[idx1],mymesh.y[idx2]],[mymesh.z[idx1],mymesh.z[idx2]],"k.-")
+    PyPlot.plot3D([1,1],[1,1],[1,1],"k.-")
+    PyPlot.text3D(mymesh.x[idx1].+rand()/30,mymesh.y[idx1].+rand()/30,mymesh.z[idx1].+rand()/30,"$idx1",ha="center",va="center")
+    # sleep(0.1)
+end
+
+for ijoint = 1:length(myjoint[:,1])
+    idx2 = Int(myjoint[ijoint,2])
+    idx1 = Int(myjoint[ijoint,3])
+    PyPlot.plot3D([mymesh.x[idx1],mymesh.x[idx2]],[mymesh.y[idx1],mymesh.y[idx2]],[mymesh.z[idx1],mymesh.z[idx2]],"r.-")
+    PyPlot.text3D(mymesh.x[idx1].+rand()/30,mymesh.y[idx1].+rand()/30,mymesh.z[idx1].+rand()/30,"$idx1",color="r",ha="center",va="center")
+    PyPlot.text3D(mymesh.x[idx2].+rand()/30,mymesh.y[idx2].+rand()/30,mymesh.z[idx2].+rand()/30,"$idx2",color="r",ha="center",va="center")
+    # sleep(0.1)
+end
+PyPlot.xlabel("x")
+PyPlot.ylabel("y")
+PyPlot.zlabel("z")
 
 RPM = 34.0
 omega = RPM*2*pi/60
@@ -209,7 +230,7 @@ PyPlot.ylabel("Torque (kN-m)")
 PyPlot.xlim([0.0,25.0])
 PyPlot.ylim([-25.0,150.0])
 PyPlot.legend()
-# PyPlot.savefig("$(path)/../figs/34m_fig3_6.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig3_6.pdf",transparent = true)
 
 ##########################################
 #### Power Plot
@@ -227,7 +248,7 @@ PyPlot.xlim([0.0,25.0])
 PyPlot.ylim([-100.0,600.0])
 PyPlot.grid(linestyle="--")
 PyPlot.legend()
-# PyPlot.savefig("$(path)/../figs/34m_fig3_7.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig3_7.pdf",transparent = true)
 
 ##########################################
 #### CP Plot
@@ -254,7 +275,7 @@ PyPlot.ylabel("Power Coefficient")
 PyPlot.xlim([0,13])
 PyPlot.ylim([0,0.8])
 PyPlot.legend()
-# PyPlot.savefig("$(path)/../figs/34m_fig3_8.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig3_8.pdf",transparent = true)
 
 ##########################################
 #### CT Plot
@@ -268,7 +289,7 @@ PyPlot.ylabel("Thrust Coefficient")
 PyPlot.xlim([0,15])
 PyPlot.ylim([0,0.9])
 PyPlot.legend()
-# PyPlot.savefig("$(path)/../figs/34m_figCT_34RPM.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_figCT_34RPM.pdf",transparent = true)
 
 ##########################################
 ############ AeroElastic #################
@@ -392,43 +413,43 @@ SNL34m_5_4_FlatwiseStress = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Da
 
 # Plots
 PyPlot.figure()
-PyPlot.plot(t.-offsetTime,flatwise_stress1[:,end-5]./1e6,"-",color=plot_cycle[1],label = "OWENS Blade 1")
+PyPlot.plot(t.-offsetTime,flatwise_stress1[:,end-6]./1e6,"-",color=plot_cycle[1],label = "OWENS Blade 1")
 PyPlot.plot(SNL34m_5_4_FlatwiseStress[:,1].-0.8,SNL34m_5_4_FlatwiseStress[:,2],"k-",label = "Experimental")
 PyPlot.xlabel("Time (s)")
 PyPlot.ylabel("Flapwise Stress (MPa)")
 PyPlot.xlim([0,SNL34m_5_4_FlatwiseStress[end,1]])
 # PyPlot.ylim([-60,0])
 PyPlot.legend(loc = (0.06,1.0),ncol=2)
-# PyPlot.savefig("$(path)/../figs/34m_fig5_4_NormalOperation_flapwise_Blade2Way.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig5_4_NormalOperation_flapwise_Blade2Way.pdf",transparent = true)
 
 println("here")
 exp_std_flap = Statistics.std(SNL34m_5_4_FlatwiseStress[:,2])
 println("exp_std_flap $exp_std_flap")
 exp_mean_flap = Statistics.mean(SNL34m_5_4_FlatwiseStress[:,2])
 println("exp_mean_flap $exp_mean_flap")
-sim_std_flap = Statistics.std(flatwise_stress1[:,end-5]./1e6)
+sim_std_flap = Statistics.std(flatwise_stress1[:,end-6]./1e6)
 println("sim_std_flap $sim_std_flap")
-sim_mean_flap = Statistics.mean(flatwise_stress1[:,end-5]./1e6)
+sim_mean_flap = Statistics.mean(flatwise_stress1[:,end-6]./1e6)
 println("sim_mean_flap $sim_mean_flap")
 
 SNL34m_5_4_LeadLagStress = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.4AML.csv",',',skipstart = 1)
 
 PyPlot.figure()
-PyPlot.plot(t.-offsetTime,lag_stress1[:,end-5]./1e6,"-",color=plot_cycle[1],label = "OWENS Blade 1")
+PyPlot.plot(t.-offsetTime,lag_stress1[:,end-6]./1e6,"-",color=plot_cycle[1],label = "OWENS Blade 1")
 PyPlot.plot(SNL34m_5_4_LeadLagStress[:,1],SNL34m_5_4_LeadLagStress[:,2],"k-",label = "Experimental")
 PyPlot.xlabel("Time (s)")
 PyPlot.ylabel("Lead-Lag Stress (MPa)")
 PyPlot.xlim([0,SNL34m_5_4_LeadLagStress[end,1]])
 PyPlot.legend(loc = (0.06,1.0),ncol=2)
-# PyPlot.savefig("$(path)/../figs/34m_fig5_4_NormalOperation_LeadLag_Blade2Way.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig5_4_NormalOperation_LeadLag_Blade2Way.pdf",transparent = true)
 
 exp_std_lag = Statistics.std(SNL34m_5_4_LeadLagStress[:,2])
 println("exp_std_lag $exp_std_lag")
 exp_mean_lag = Statistics.mean(SNL34m_5_4_LeadLagStress[:,2])
 println("exp_mean_lag $exp_mean_lag")
-sim_std_lag = Statistics.std(lag_stress1[:,end-5]./1e6)
+sim_std_lag = Statistics.std(lag_stress1[:,end-6]./1e6)
 println("sim_std_lag $sim_std_lag")
-sim_mean_lag = Statistics.mean(lag_stress1[:,end-5]./1e6)
+sim_mean_lag = Statistics.mean(lag_stress1[:,end-6]./1e6)
 println("sim_mean_lag $sim_mean_lag")
 
 ##########################################
@@ -447,7 +468,7 @@ PyPlot.xlabel("Time (s)")
 PyPlot.xlim([0,100])
 PyPlot.ylabel("Torque (kN-m)")
 PyPlot.legend()#loc = (0.06,1.0),ncol=2)
-# PyPlot.savefig("$(path)/../figs/34m_fig5_32Way.pdf",transparent = true)
+PyPlot.savefig("$(path)/figs/34m_fig5_32Way.pdf",transparent = true)
 
 # #Extract and Plot Frequency - Amplitude
 
@@ -583,4 +604,4 @@ end
 
 azi=aziHist#./aziHist*1e-6
 saveName = "$path/vtk/SNL34mVAWTNormalOperation"
-OWENS.OWENSFEA_VTK(saveName,t,uHist,system,assembly,sections;scaling=1,azi,userPointNames,userPointData)
+# OWENS.OWENSFEA_VTK(saveName,t,uHist,system,assembly,sections;scaling=1,azi,userPointNames,userPointData)
diff --git a/examples/SNL34m/SNL34m_Inputs.yml b/examples/SNL34m/SNL34m_Inputs.yml
@@ -23,7 +23,7 @@ turbulentInflow:
 controlParameters:
   controlStrategy: constantRPM # TODO: incorporate the others
   RPM: 34.0 #RPM
-  numTS: 5000 #
+  numTS: 1000 #
   delta_t: 0.05 # s
 
 
@@ -35,7 +35,7 @@ AeroParameters:
   adi_rootname: "/SNL34m"
 
 structuralParameters:
-  structuralModel: GX #GX, TNB, ROM
+  structuralModel: ROM #GX, TNB, ROM
   nonlinear: false #TODO: propogate
   ntelem: 20 #tower elements in each 
   nbelem: 60 #blade elements in each 

diff --git a/src/PostProcessing.jl b/src/PostProcessing.jl
@@ -314,7 +314,7 @@ end
 
 function calcSF(stress,SF_ult,SF_buck,lencomposites_span,plyprops,
     precompinput,precompoutput,lam_in,eps_x,eps_z,eps_y,kappa_x,
-    kappa_y,kappa_z,numadIn;failmethod = "maxstress",CLT=false,upper=true,layer=-1)
+    kappa_y,kappa_z,numadIn;failmethod = "maxstress",CLT=false,upper=true,layer=-1,calculate_fatigue=true)
     topstrainout = zeros(length(eps_x[1,:,1]),lencomposites_span,length(lam_in[1,:]),9) # time, span, lam, x,y, Assumes you use zero plies for sections that aren't used
 
     damage = zeros(lencomposites_span,length(lam_in[1,:])) #span length, with number of laminates, with number of plies NOTE: this assumes the number of plies is constant across all span and laminate locations
@@ -441,28 +441,30 @@ function calcSF(stress,SF_ult,SF_buck,lencomposites_span,plyprops,
             end
 
             # Miner's Damage
-            damage_layers = zeros(length(lam.nply))
-            for ilayer = 1:length(lam.nply)
-                #TODO: multiple mean ranges and goodman correction, also non-principal stress
-                stressForFatigue = stress_eachlayer[:,ilayer,1]
-                Ncycles,meanIntervals,rangeIntervals,_ = rainflow(stressForFatigue;nbins_range=20,nbins_mean=1,m=3,Teq=1)
-                imean = 1
-                cyclesatStressRanges = Ncycles[imean,:]
-                stress_levels = (rangeIntervals[1:end-1] .+ rangeIntervals[2:end])./2 #from intervals to mean between the interval
-
-                SN_stressMpa1 = SN_stressMpa[ilayer]
-                Log_SN_cycles2Fail1 = Log_SN_cycles2Fail[ilayer]
-
-                if Log_SN_cycles2Fail1[2]>Log_SN_cycles2Fail1[1]
-                    reverse!(SN_stressMpa1)
-                    reverse!(Log_SN_cycles2Fail1)
-                end
+            if calculate_fatigue
+                damage_layers = zeros(length(lam.nply))
+                for ilayer = 1:length(lam.nply)
+                    #TODO: multiple mean ranges and goodman correction, also non-principal stress
+                    stressForFatigue = stress_eachlayer[:,ilayer,1]
+                    Ncycles,meanIntervals,rangeIntervals,_ = rainflow(stressForFatigue;nbins_range=20,nbins_mean=1,m=3,Teq=1)
+                    imean = 1
+                    cyclesatStressRanges = Ncycles[imean,:]
+                    stress_levels = (rangeIntervals[1:end-1] .+ rangeIntervals[2:end])./2 #from intervals to mean between the interval
+
+                    SN_stressMpa1 = SN_stressMpa[ilayer]
+                    Log_SN_cycles2Fail1 = Log_SN_cycles2Fail[ilayer]
+
+                    if Log_SN_cycles2Fail1[2]>Log_SN_cycles2Fail1[1]
+                        reverse!(SN_stressMpa1)
+                        reverse!(Log_SN_cycles2Fail1)
+                    end
 
-                logcycles2fail = safeakima(SN_stressMpa1.*1e6,Log_SN_cycles2Fail1,stress_levels)
-                cycles2fail = 10.0 .^ logcycles2fail
-                damage_layers[ilayer] = sum(cyclesatStressRanges./cycles2fail)
+                    logcycles2fail = safeakima(SN_stressMpa1.*1e6,Log_SN_cycles2Fail1,stress_levels)
+                    cycles2fail = 10.0 .^ logcycles2fail
+                    damage_layers[ilayer] = sum(cyclesatStressRanges./cycles2fail)
+                end
+                damage[i_station,j_lam] = maximum(damage_layers)
             end
-            damage[i_station,j_lam] = maximum(damage_layers)
         end
     end
     return topstrainout,damage
@@ -625,7 +627,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
     LE_U_idx=1,TE_U_idx=6,SparCapU_idx=3,ForePanelU_idx=2,AftPanelU_idx=5,
     LE_L_idx=1,TE_L_idx=6,SparCapL_idx=3,ForePanelL_idx=2,AftPanelL_idx=5,
     Twr_LE_U_idx=1,Twr_LE_L_idx=1,throwawayTimeSteps=1,delta_t=0.001,AD15bldNdIdxRng=nothing,AD15bldElIdxRng=nothing,
-    strut_precompoutput=nothing,strut_precompinput=nothing,plyprops_strut=nothing,numadIn_strut=nothing,lam_U_strut=nothing,lam_L_strut=nothing)
+    strut_precompoutput=nothing,strut_precompinput=nothing,plyprops_strut=nothing,numadIn_strut=nothing,lam_U_strut=nothing,lam_L_strut=nothing,calculate_fatigue=true)
 
     # Linearly Superimpose the Strains
     epsilon_x_hist = epsilon_x_hist_ps#copy(epsilon_x_hist_ps).*0.0
@@ -797,7 +799,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
     topstrainout_blade_U,topDamage_blade_U = calcSF(stress_U,SF_ult_U,SF_buck_U,length(bld_precompinput),plyprops_bld,
     bld_precompinput,bld_precompoutput,lam_U_bld,eps_x_bld,eps_z_bld,eps_y_bld,kappa_x_bld,
-    kappa_y_bld,kappa_z_bld,numadIn_bld;failmethod = "maxstress",upper=true)
+    kappa_y_bld,kappa_z_bld,numadIn_bld;failmethod = "maxstress",upper=true,calculate_fatigue)
 
     if verbosity>0
         println("Composite Ultimate and Buckling Safety Factors")
@@ -816,7 +818,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
     topstrainout_blade_L,topDamage_blade_L = calcSF(stress_L,SF_ult_L,SF_buck_L,length(bld_precompinput),plyprops_bld,
     bld_precompinput,bld_precompoutput,lam_L_bld,eps_x_bld,eps_z_bld,eps_y_bld,kappa_x_bld,
-    kappa_y_bld,kappa_z_bld,numadIn_bld;failmethod = "maxstress",upper=false)
+    kappa_y_bld,kappa_z_bld,numadIn_bld;failmethod = "maxstress",upper=false,calculate_fatigue)
 
     if verbosity>0
         println("\n\nLOWER BLADE SURFACE")
@@ -838,7 +840,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
         topstrainout_strut_U,topDamage_strut_U = calcSF(stress_U_strut,SF_ult_U_strut,SF_buck_U_strut,length(strut_precompinput),plyprops_strut,
         strut_precompinput,strut_precompoutput,lam_U_strut,eps_x_strut,eps_z_strut,eps_y_strut,kappa_x_strut,
-        kappa_y_strut,kappa_z_strut,numadIn_strut;failmethod = "maxstress",upper=true)
+        kappa_y_strut,kappa_z_strut,numadIn_strut;failmethod = "maxstress",upper=true,calculate_fatigue)
 
         if verbosity>0
             println("Composite Ultimate and Buckling Safety Factors")
@@ -857,7 +859,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
         topstrainout_strut_L,topDamage_strut_L = calcSF(stress_L_strut,SF_ult_L_strut,SF_buck_L_strut,length(strut_precompinput),plyprops_strut,
         strut_precompinput,strut_precompoutput,lam_L_strut,eps_x_strut,eps_z_strut,eps_y_strut,kappa_x_strut,
-        kappa_y_strut,kappa_z_strut,numadIn_strut;failmethod = "maxstress",upper=false)
+        kappa_y_strut,kappa_z_strut,numadIn_strut;failmethod = "maxstress",upper=false,calculate_fatigue)
 
         if verbosity>0
             println("\n\nLOWER STRUT SURFACE")
@@ -888,7 +890,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
     topstrainout_tower_U,topDamage_tower_U = calcSF(stress_TU,SF_ult_TU,SF_buck_TU,length(twr_precompinput),plyprops_twr,
     twr_precompinput,twr_precompoutput,lam_U_twr,eps_x_twr,eps_z_twr,eps_y_twr,kappa_x_twr,
-    kappa_y_twr,kappa_z_twr,numadIn_twr;failmethod = "maxstress",upper=true)
+    kappa_y_twr,kappa_z_twr,numadIn_twr;failmethod = "maxstress",upper=true,calculate_fatigue)
 
     println("\n\nUPPER TOWER")
     printsf_twr(verbosity,lam_U_twr,SF_ult_TU,SF_buck_TU,length(twr_precompinput),Twr_LE_U_idx,topDamage_tower_U,delta_t,total_t)
@@ -899,7 +901,7 @@ function extractSF(bld_precompinput,bld_precompoutput,plyprops_bld,numadIn_bld,l
 
     topstrainout_tower_L,topDamage_tower_L = calcSF(stress_TL,SF_ult_TL,SF_buck_TL,length(twr_precompinput),plyprops_twr,
     twr_precompinput,twr_precompoutput,lam_U_twr,eps_x_twr,eps_z_twr,eps_y_twr,kappa_x_twr,
-    kappa_y_twr,kappa_z_twr,numadIn_twr;failmethod = "maxstress",upper=false)
+    kappa_y_twr,kappa_z_twr,numadIn_twr;failmethod = "maxstress",upper=false,calculate_fatigue)
 
     println("\n\nLower TOWER")
     printsf_twr(verbosity,lam_L_twr,SF_ult_TL,SF_buck_TL,length(twr_precompinput),Twr_LE_L_idx,topDamage_tower_L,delta_t,total_t)