- Get link
- X
- Other Apps
Featured Post
- Get link
- X
- Other Apps
Hello, everyone. Today we'll continue our discussion about Brownian Carnot engine, by directly build a virtual one. To understand the code, you need to know how to do a computer simulation of a Langevin equation. Please refer to our previous post to check it out.
/*------------Divider--------------*/
So let's see our code written in MATLAB...
/*------------Divider--------------*/
So let's see our code written in MATLAB...
function Result = Brownian_Carnot_Module %% permanent setting tic; global dt Dt Dt_obs kB R eta zeta T0 A perN_hit perN_obs N_multiple m N_record_init dt = 1E-7; Dt = 1E-5; Dt_obs = 1E-4; kB = 1.38E-23; R = 500E-9; eta = 8.90E-4; zeta = 6*pi*eta*R; rho = 1040; m = 4/3*pi*R^3*rho; T0 = 300; A = sqrt(2*kB*T0*zeta*dt); perN_hit = round(Dt/dt); perN_obs = round(Dt_obs/dt); N_multiple = 300; %initiation delT_init = 0.02; T_set_init = 300; kappa_init = 1E-6; IsoComp = 4; AdiaComp = 4; delT_isoComp = 0.03; delT_adiaComp = 0.03; delT_isoExp = 0.03; delT_adiaExp = 0.03; State.x_record = []; State.T_record= []; State.T_set = []; State.kappa_record = []; State.sign = []; State.time = 0; State.ProcessRecord = {}; State.QRecord = []; %% Carnot Engine State = Carnot_Initiation(State, delT_init, T_set_init , kappa_init); toc; State = Carnot_Process(State, 'I', delT_isoComp , IsoComp); toc; State = Carnot_Process(State, 'A', delT_adiaComp , AdiaComp); toc; State = Carnot_Process(State, 'I', delT_isoExp , 1/IsoComp); toc; State = Carnot_Process(State, 'A', delT_adiaExp , 1/AdiaComp); toc; %% time course figure summary x_record = State.x_record; kappa_summary = State.kappa_record; v_record = abs((x_record(:, 2:end) - x_record(:, 1:(end-1)))/Dt_obs).*repmat(sqrt(Lt(kappa_summary(1:(end-1)))), N_multiple ,1); Sx = EntropyCal(x_record); Sv = EntropyCal(v_record); S = Sx(1:(end-1)) + Sv; %calculate entropy Time = Dt_obs:Dt_obs:(State.time - Dt_obs); Temp = State.T_record; T_set = State.T_set; Result.x_record = x_record; Result.v_record = v_record; Result.kappa_summary = kappa_summary; Sign = State.sign; Result.Sx = Sx; Result.Sv = Sv; Result.S = S; Result.Time = Time; Result.Temp = Temp; Result.ProcessRecord = State.ProcessRecord; figure; [ax,h1,h2] = plotyy(Time, S, Time, Temp(1:(end-1))); hold on; set(h1,'linewidth',3); set(h2,'linewidth',2); set(ax, 'FontSize', 14); xlabel('Time','FontSize',16); ylabel('Entropy','FontSize',16); yLim = get(gca, 'YLim'); [nProcess,~]=size(State.ProcessRecord); for i = 1:1:nProcess line([State.ProcessRecord{i,2}, State.ProcessRecord{i,2}], yLim,'LineWidth',2); text(State.ProcessRecord{i,2},yLim(2),State.ProcessRecord{i,1},'FontSize',20); end toc; %% Other Plots %T-kappa figure; scatter(kappa_summary, Temp,[],Sign); hold on; scatter(kappa_summary, T_set ,[] ,Sign,'filled'); xlabel('\kappa (N/m)','FontSize',16); ylabel('Temperature(K)','FontSize',16); set(gca, 'FontSize', 14); %F-kappa figure; F_kappa_exp = mean(x_record.^2, 1)/2; F_kappa_theo = 1/2*kB*T_set./kappa_summary; scatter(kappa_summary, F_kappa_exp,[],Sign); hold on; scatter(kappa_summary, F_kappa_theo,[],Sign,'filled'); ylabel('F_{\kappa} (m^2)','FontSize',16); xlabel('\kappa (N/m)','FontSize',16); set(gca, 'FontSize', 14); %T-S figure; S_baseline = mean(S(1:N_record_init)); S_theo = S_baseline + log(T_set/T_set_init)-1/2*log(kappa_summary/kappa_init); scatter(S, Temp(1:end-1),[],Sign(1:end-1,:)); hold on; scatter(S_theo, T_set ,[] ,Sign,'filled'); xlabel('S/k_B (1/K)','FontSize',16); ylabel('Temperature(K)','FontSize',16); set(gca, 'FontSize', 14); S_cal = 0.5*log(2*pi*kB*T0/kappa_init)+1/(2*pi)+0.5*log(2*pi*kB*T0/m)+1/(2*pi); display([S_baseline, S_cal]); %efficiency dW = F_kappa_exp(1:end-1).*(kappa_summary(2:end)-kappa_summary(1:end-1)); dQ = kB*(Temp(2:end)-Temp(1:end-1))- dW; recordOrNot = logical(State.QRecord(1:end-1)); efficiency = - sum(dW)/sum(dQ(recordOrNot)); efficiency_theo = 1- 1/sqrt(AdiaComp); display([efficiency, efficiency_theo]); Result.efficiency = efficiency; Result.efficiency_theo = efficiency_theo; Result.dQ =dQ; Result.dW =dW; Result.recordOrNot = recordOrNot; %% end beep; end function correction = Lt(kappa) global Dt_obs zeta m f = 1/Dt_obs; fp = zeta/m; fk = kappa/(2*pi*zeta); f0 = sqrt(fp*fk); f1 = sqrt(fp^2/4 - f0.^2); correction = 1./(2.*f.^2).*(1./f0.^2 + 1./f1.*(exp(-f1/f - fp/(2*f))./(fp+2*f1) - exp(f1/f - fp/(2*f))./(fp-2*f1))).^(-1); end function S = EntropyCal(record) [~, nLength] = size(record); S = zeros(1, nLength); for i = 1:1:nLength h = histogram(record(:,i), 25,'Normalization','pdf'); pdf = h.Values; d = h.BinWidth; close all; S(i) = - nansum(pdf.*log(pdf))*d; end end function State = Carnot_Initiation(State, delT, TempSet, kappa_init) %% global dt perN_obs T0 kB N_multiple A zeta perN_hit Dt N_record_init display('initiation...'); State.ProcessRecord{1,1}='initiation';State.ProcessRecord{1,2}=State.time; N_init = delT/dt; N_record_init = round(N_init/perN_obs); if TempSet < T0 display('Error!');display('The setted temperature should always be higher than 300K.') return end T_set_init = TempSet - T0; A1_init = sqrt(2*kB*T_set_init*zeta*Dt); %% x_record_init = zeros(N_multiple, N_record_init); for j = 1:1:N_multiple x_temp = zeros(1, N_init); x_temp(1) = sqrt(kB*TempSet/kappa_init)*randn(1,1); for i = 2:1:N_init x_old = x_temp(i-1); x_temp(i) = x_old + 1/zeta*(-dt*(kappa_init*x_old)+A*randn(1,1)); if mod(i, perN_hit) == 0 x_temp(i) = x_temp(i) + 1/zeta*(A1_init*randn(1,1)); end if mod(i, perN_obs) == 0 x_record_init(j, i/perN_obs) = x_temp(i); end end end %% T_eff_init = kappa_init*mean(x_record_init.^2,1)/kB; State.x_record = [State.x_record x_record_init]; State.T_record = [State.T_record T_eff_init]; State.T_set = [State.T_set TempSet*ones(1, N_record_init)]; State.kappa_record = [State.kappa_record kappa_init*ones(1, N_record_init)]; State.sign = vertcat(State.sign, repmat([0 0 0],N_record_init,1)); State.time = State.time + delT; State.QRecord = [State.QRecord zeros(1, N_record_init)]; End function State = Carnot_Process(State, IsoAdia, delT, multiply) %% global dt perN_obs T0 kB N_multiple A zeta perN_hit Dt N = delT/dt; N_record_Process = round(N/perN_obs); kappa_init = State.kappa_record(end); kappa_endProcess = multiply*kappa_init; kappa_Process = kappa_init:((kappa_endProcess-kappa_init)/(N-1)):kappa_endProcess; T_init = State.T_set(end); %% if strcmp(IsoAdia, 'I') if multiply >1; display('isothermal compression ...'); sign = [0 0 1]; State.ProcessRecord{end+1,1}='isothermal compression'; State.ProcessRecord{end,2}=State.time; State.QRecord = [State.QRecord zeros(1, N_record_Process)]; else display('isothermal expansion ...'); sign = [1 0 0]; State.ProcessRecord{end+1,1}='isothermal expansion'; State.ProcessRecord{end,2}=State.time; State.QRecord = [State.QRecord ones(1, N_record_Process)]; end T_set_Process = (T_init - T0)*ones(1, N); A1_Process = sqrt(2*kB*T_set_Process*zeta*Dt); elseif strcmp(IsoAdia, 'A') if multiply >1; display('adiabatic compression ...'); sign = [0 1 0]; State.ProcessRecord{end+1,1}='adiabatic compression'; State.ProcessRecord{end,2}=State.time; State.QRecord = [State.QRecord ones(1, N_record_Process)]; else display('adiabatic expansion ...'); sign = [1 0 1]; State.ProcessRecord{end+1,1}='adiabatic expansion'; State.ProcessRecord{end,2}=State.time; State.QRecord = [State.QRecord ones(1, N_record_Process)]; end Const = T_init^2/kappa_init; T_set_Process = round(sqrt(Const*kappa_Process) - T0, 5); A1_Process = sqrt(2*kB*T_set_Process*zeta*Dt); else display('incomprehensible input'); end x_record_Process = zeros(N_multiple, N_record_Process); for j = 1:1:N_multiple x_temp = zeros(1, N); x_old = State.x_record(j, end); x_temp(1) = x_old + 1/zeta*(-dt*(kappa_Process(1)*x_old)+A*randn(1,1)); for i = 2:1:N x_old = x_temp(i-1); x_temp(i) = x_old + 1/zeta*(-dt*(kappa_Process(i)*x_old)+A*randn(1,1)); if mod(i, perN_hit) == 0 x_temp(i) = x_temp(i) + 1/zeta*(A1_Process(i)*randn(1,1)); end if mod(i, perN_obs) == 0 x_record_Process(j, i/perN_obs) = x_temp(i); end end end kappa_Process_abre = kappa_Process(perN_obs: perN_obs: N); T_eff_Process = kappa_Process_abre.*mean(x_record_Process.^2,1)/kB; T_set_Process_abre = T_set_Process(perN_obs: perN_obs: N); %% State.x_record = [State.x_record x_record_Process]; State.T_record = [State.T_record T_eff_Process]; State.kappa_record = [State.kappa_record kappa_Process_abre]; State.T_set = [State.T_set T_set_Process_abre+T0]; State.sign = vertcat(State.sign, repmat(sign,N_record_Process,1)); State.time = State.time + delT; end
It looks a little complicated. Our simulation is composed of an initial waiting period and a full Carnot cycle. This is run for 300 times to take average and calculate entropy. The entropy is calculated from the probability density function of the location and speed of the microsphere using Shannon formula (that is,
). After these simulation, we made some plot to display our results.
In fact, it would be great if our readers try to modify these codes. For example, by adjusting the rate of compression and expansion, or change the order of each process.
Before we delve to our results, I have to briefly raise 3 issues. First, in the above code, we estimate the baseline entropy by taking the average of the entropy measured in initial waiting period. This is in fact not necessary. The entropy of the system could be calculated explicitly as
In fact, it would be great if our readers try to modify these codes. For example, by adjusting the rate of compression and expansion, or change the order of each process.
Before we delve to our results, I have to briefly raise 3 issues. First, in the above code, we estimate the baseline entropy by taking the average of the entropy measured in initial waiting period. This is in fact not necessary. The entropy of the system could be calculated explicitly as
This analytical solution reads -20.3568, which is basically the same as our average in initial waiting period, -20.2964. This is left for the readers to prove that.
Secondly, the evaluation of efficiency is a little bit different from classical Carnot engine. In classical Carnot engine, the efficiency is calculated as:
However, in its stochastic counterpart, the heat flow during the namely "adiabatic" process must be taken into consideration. Thus the efficiency of a Brownian Carnot engine should be calculated as:
Thirdly, since we observe the microsphere and estimate its speed by a rather long time interval, it is necessary to correct our average speed to approximate the instantaneous speed. The authors have proved the validity of their correction formula.
And finally, let's see our results:
This is the summary of our experiment time course. As expected, the entropy decreases during isothermal compression, and nearly holds the same during 2 adiabatic processes.
The filled color lines indicate what a classical Carnot engine would behave in our experimental setting. There are some fluctuations but our stochastic engine acts like his classical counterpart. However, this would not be the case if we hasten the rate of compression and expansion. What do you expect to change? What is the source of irreversibility in Brownian Carnot engine?
This is the F-κ plot, which you could interpret as a P-V plot. The area enclosed in the cycle is the work our engine done.
And this is a T-S plot. In our simulation here, the efficiency is 0.4427, which is close to the efficiency of its classical counterpart, 0.50, because we compress and expand the system rather slowly and take average from 300 cycles. Efficiency above classical limit is possible, however, due to the stochastic nature of our system, and it is in fact a big issue to investigate the distribution of efficiency of our system. But we will stop here for it takes time to explore.
- Get link
- X
- Other Apps
Comments
Such a nice post, thanks. iam.
ReplyDeletehow to calculate the efficiency?
ReplyDelete