diff --git a/examples/SNL34m/SNL34mVAWTNormalOperation.jl b/examples/SNL34m/SNL34mVAWTNormalOperation.jl index 1dca2737..3c8ad510 100644 --- 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 index 925b210a..5317e990 100644 --- 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 index 2d043892..74b9edad 100644 --- 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)