MATLAB Answers

Help with a complex 2nd order, nonlinear ODE BVP with complex boundary conditions

3 ビュー (過去 30 日間)
Whitewater
Whitewater 2015 年 5 月 13 日
編集済み: Whitewater 2015 年 5 月 13 日
I am attempting to solve the Poisson-Boltzmannn equation, which describes the distribution of potential and concentration in relation to a charged surface in a liquid, for a cylinder. The equation this presents is a 2nd order nonlinear ODE that cannot be represented as an IVP, but instead must be solved as a BVP. The initial ODE from this relationship is:
y''=2*z*F*C_bulk/ereo/Vt*exp(mu_bulk-mu_x)*sinh(z*yd)-y'/x
I know in order to use bvp4c in Matlab I have to break this into a system of first order ODEs, so I have set y1=y' to generate the following system:
(1) y1=y'
(2) y1'=y2=2*z*F*c_bulk/ereo/Vt*exp(mu_bulk-mu_x)*sinh(z*y_d)-y1/x
with the boundary cnoditions:
(1) y1 @ x(1,end) = 0
(2) y @ x(1,1) = y_d
There are additional equations that define the system, which include:
mu_x = eta_x.*(8-9*eta_x+3*eta_x.^2)./((1-eta_x).^3)
eta_x = vion*(c_cation+c_anion)
c_cation = c_bulk*exp(-z*y+mu_bulk-mu_x)
c_anion = c_bulk*exp(z*y+mu_bulk-mu_x)
c_x = c_cation+c_anion
q_x = z*c_cation-z*c_anion
c_tot = sum(c_x*pi*(x(1,n+1).^2-x(1,n).^2))
q_tot = sum(q_x*pi*(x(1,n+1).^2-x(1,n).^2))
y_d = Vcell/2/Vt-q_tot*dion*F/2/ereo/Vt
In these equations, variables directly related to distance from the charged surface, x, are row vectors with identical length to the vector, x. These variables include: y1, y2, c_cation, c_anion, c_x, q_x, mu_x, and eta_x.
Some values are constant and they are included in the code I have already generated which will be posted below.
I have generated a code that will predict reasonable values for y if given reasonable mu_x and y_d values (see the original system of ODEs), but the issues I am having are:
(1) How do I include the additional equations in the existing model,and
(2) These equations must be solved iteratively and both depend on and affect the solution of the system of ODEs. How do I perform this iterative calculation?
Here is the code I have generated:
-----------------------------------------------------------------------------------------
function dydx = mybvp(x,y)
global F z c_bulk ereo Vt mu_bulk mu_x
dydx = [ y(2,:).*10^-12; 2*z*F*c_bulk/ereo/Vt*exp(mu_bulk-mu_x(1,:))*sinh(y_d)];
end
function in = mybvpin (x)
in=[(9*10^-6.*((x*10^12)^2))-(0.0137.*(x*10^12))+(y_d)
(1.8*10^-5.*(x*10^12))-0.0137];
end
clear; clc;
global MW %make the variable, MW, usable in all functions
MW=60;%input('Enter the molar mass of the salt ion in g/mol '); %retrieve the molar mass from user
global z %make the variable, z, usable in all functions
z=1;%input('Enter the valence number of the salt ions '); %retrieves the valence number of the salt
global c_bulk %make the variable, c_bulk, usable in all functions
c_bulk=500/MW; %input('Enter the bulk salt concentration in mg/L ')/MW; %retrieve the bulk concentration from user and convert to mol/m3
global dion %make the variable, dion, usable in all functions
dion=700*10^-12; %input('Enter the hydrated ionic diameter in pm ')/10^12; %retrieve the hydrated ionic diameter and convert it to m
global rion %make the variable, rion, usable in all functions
rion=dion/2; %calculate the hydrated ionic radius
global del %make the variable, del, usable in all functions
del=rion; %set the Stern layer thickness
global d %make the variable, d, usable in all functions
d=1*10^-12; %input('Enter the pore diameter in nm ')/10^9; %retrieve the pore diameter and store it as d in units of m
global Vcell %make the variable, Vcell, usable in all functions
Vcell=1.2; %input('Enter the applied voltage in volts '); %retrieve the applied voltage
global vion %make the variable, volume, usable in all functions
vion=pi/6*dion^3; %calculate the volume of one ion from its diameter
global Av %make the variable, Av, usable in all functions
Av=6.022*10^23; %set Avogardo's number
global eta_bulk %make the variable, eta_b, usable in all functions
eta_bulk=vion*c_bulk*2*Av; %calculate the volume fraction of ions
global mu_bulk %make the variable, mu_bulk, usable in all functions
mu_bulk=eta_bulk*(8-9*eta_bulk+3*eta_bulk^2)/((1-eta_bulk)^3); %calculate the bulk electrochemical potential correction in accordance with the CS correction
global e %make the variable, e, usable in all functions
e=-1.60218*10^-19; %set the charge of an electron
global ereo %make the variable, ereo, usable in all functions
ereo=78.64*8.854*10^-12; %set the dielectric permittivity of water
global Vt %make the variable, Vt, usable in all functions
Vt=0.025692577; %set the thermal voltage
global F %make the variable, F, usable in all functions
F=96485; %set the Faraday costant
global y_0 %make the variable, y_0 usable in all functions
y_0=Vcell/2/Vt; %Input the equation for the potential at the charged surface
solinit = bvpinit(linspace(0,155*10^-12,100),@mybvpin);
sol = bvp4c(@mybvp,@mybvpbc,solinit);
x=linspace(0,155*10^-12,20);
y=deval(sol,x);
plot(x,y(1,:));
Thanks in advance for your assistance!

  2 件のコメント

Torsten
Torsten 2015 年 5 月 13 日
Could you clarify where the additional equations depend on the solution variables y1 and y2 ?
Best wishes
Torsten.
Whitewater
Whitewater 2015 年 5 月 13 日
I'm sorry. There was a typo in my original equations. y_d in the equations for c_cation and c_anion should instead be y, which is also a vector with the same length as x.
c_cation = c_bulk*exp(-z*y+mu_bulk-mu_x)
c_anion = c_bulk*exp(z*y+mu_bulk-mu_x)
Also, perhaps I am a little confused about what is being returned by the function. Above the equations are dependent on the local potential, y, which I believe is being returned by the bvp4c routine in the first row with how I have the problem set up. This is y where y1=y' and y2=y''.

サインインしてコメントする。

回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by