77r_coord = ones(lR ,lr_N ).*logspace(-0.5 ,lr_length ,lr_N );
88% calculating the shear
99varsigmadot_r = f_gammadot_r(r_coord ,R ,U ,lR ,lr_N );
10- f_r = sf_carreau(v_nc ,v_lambda ,varsigmadot_r );
1110[xcon ,ycon ] = meshgrid(t ,r_coord(1 ,: ));
1211ycon = log10(ycon );
1312clevels = 50 ;
1615ntau_r = tau_r / max_tau_r ;
1716
1817% f contour figure
19- figure(5 )
18+ figure(1 )
2019hold on ;
2120xlabel(' $t / t_c$' , ' Interpreter' , ' Latex' , ' FontSize' , 20 );
2221ylabel(' log$_{10}(\it{r}/R_o)$' ,' Interpreter' ,' Latex' ,' FontSize' ,24 );
3938xa.LineWidth = 1.5 ;
4039box on ;
4140plot(t , log10(R ),lm ,' LineWidth' ,3 );
42- contourf(xcon ' ,ycon ' ,f_r ,clevels ,' edgecolor' ,' none' )
41+ contourf(xcon ' ,ycon ' ,ntau_r ,clevels ,' edgecolor' ,' none' )
4342saveas(gcf ,' ./figs/baseline/fcon_T' ,' png' )
4443
45- % shear stress contour figure
46- figure(6 )
47- hold on ;
48- xlabel(' $t / t_c$' , ' Interpreter' , ' Latex' , ' FontSize' , 20 );
49- ylabel(' log$_{10}(\it{r}/R_o)$' ,' Interpreter' ,' Latex' ,' FontSize' ,24 );
50- colormap jet ;
51- cbar = colorbar ;
52- cbar.Label.String = ' $\t au_{rr}/\mathrm{max}(\t au_{rr})$' ;
53- set(cbar ,' TickLabelInterpreter' ,' latex' );
54- pos = get(cbar ,' Position' );
55- cbar.Label.Position = [1.5 * pos(1 ) - 1.2 ];
56- cbar.Label.Rotation = 0 ;
57- cbar.Label.Interpreter = ' latex' ;
58- clim([-1 1 ])
59- xlim([0 fig_tend ]);
60- xticks(tickrange )
61- set(gcf ,' color' ,' w' );
62- set(gca ,' FontName' ,' Times' ,' FontSize' ,20 );
63- set(gca ,' TickLabelInterpreter' ,' latex' )
64- xa = gca ;
65- xa.TickLength = [.015 .015 ];
66- xa.LineWidth = 1.5 ;
67- box on ;
68- plot(t , log10(R ),lm ,' LineWidth' ,3 );
69- contourf(xcon ' ,ycon ' ,ntau_r ,clevels ,' edgecolor' ,' none' );
70- saveas(gcf ,' ./figs/baseline/taucon_T' ,' png' )
44+ % shear as a function of r (radial coordinate) calculation
45+ function sofr = f_gammadot_r(r ,R ,Rdot ,N ,M )
46+
47+ sofr = zeros(N ,M );
48+ for i = 1 : N
49+ for j = 1 : M
50+ if r(i ,j ) >= R(i ,1 )
51+ sofr(i ,j ) = - 2 * Rdot(i ,1 )*(R(i ,1 )^2 )/((r(i ,j ))^3 );
52+ else
53+ sofr(i ,j ) = NaN ;
54+ end
55+ end
56+ end
57+ end
58+
59+ function [feps_r ] = f_f_filter(f_r ,N ,M )
60+ eps = 0.01 ;
61+ feps_r = zeros(size(f_r ));
62+ for i = 1 : N
63+ for j = 1 : M
64+ if f_r(i ,j ) > 1 - eps
65+ feps_r(i ,j ) = NaN ;
66+ elseif f_r(i ,j ) < eps
67+ feps_r(i ,j ) = NaN ;
68+ else
69+ feps_r(i ,j ) = f_r(i ,j );
70+ end
71+ end
72+ end
73+
74+ end
0 commit comments