Hi Daniel and Jonathan,
I used Matlab's pdepe function and got numerical solution that is close to
my analytical solution. I've attached the code for your information.
Best regards,
Rose
On Fri, Jun 27, 2014 at 2:10 PM, yuan wang <[email protected]> wrote:
> Hi Daniel and Jonathan,
>
> I double checked the analytical solution and it is correct. Even though C
> is close to zero on the left side, it is still considerably big compared to
> Cw (C0=0.01576 and Cw=1e-3). So, it should not be neglected in computing
> the gradient. The gradient on the left boundary is close to 1 (0.9842) and
> the analytical solution matches the problem statement. If you run the
> model, there is a plot of the analytical solution.
>
> Thanks,
> Rose
>
>
>
>
> On Fri, Jun 27, 2014 at 1:33 PM, Guyer, Jonathan E. Dr. <
> [email protected]> wrote:
>
>> What Daniel is saying is to set the FiPy solution aside for the moment.
>> Your analytical solution does not appear to be consistent with your
>> boundary condition:
>>
>> A = beta*(Cl-Cw)/(D+beta*L)
>> \approx 1
>>
>> B = Cl-A*L
>> \approx 0
>>
>> Canal = A*zc+B
>> \approx zc
>>
>> (\partial Canal / \partial zc) \approx 1
>>
>> but your boundary condition states
>>
>> (\partial C / \partial z) = \beta (C - Cw) / D
>> \approx -67
>>
>> Daniel and I disagree by 3 orders of magnitude, but neither of us gets 1.
>>
>>
>> FiPy’s solution may well be wrong, but your analytical solution is not
>> consistent with your problem statement.
>>
>>
>>
>> On Jun 27, 2014, at 1:12 PM, yuan wang <[email protected]> wrote:
>>
>> > Hi Daniel,
>> >
>> > Thanks for your response. The left side boundary condition is the first
>> boundary condition shown below. So, what I did is to express dC/dz at z=0
>> as beta*(C-Cw)/D to set the constraint for the gradient on the left side. I
>> learned this from one of the examples (Convection.Robin) in the manual. I
>> noticed in the example there are square brackets around the expression,
>> which I didn't use here. If I included it, the solution doesn't sense at
>> all. Did I write the expression in a wrong way? Thanks again!
>> >
>> > <image.png>
>> >
>> > Best regards,
>> > Rose
>> >
>> >
>> > On Fri, Jun 27, 2014 at 11:45 AM, Daniel Wheeler <
>> [email protected]> wrote:
>> > Hi Rose,
>> >
>> > I was trying to debug this and it seems like there is an inconsistency
>> between the left hand boundary and the analytical solution.
>> >
>> > By my calculations A~1 and B~0, that makes the gradient of the line
>> about 1, yet the left hand boundary condition expects a gradient of about
>> 0.0666 (if C is close to 0 on the left hand side, which the analytical
>> solution wants).
>> >
>> > There may still be issues with FiPy, but I'd like to resolve that one
>> first. I may also have gotten mixed up in the calculations, but I can't see
>> where.
>> >
>> > Cheers,
>> >
>> > Daniel
>> >
>> > On Wed, Jun 25, 2014 at 4:55 PM, yuan wang <[email protected]> wrote:
>> > Hi,
>> >
>> > I tried to solve a problem with a Robin Boundary condition in fipy. The
>> physical problem is very simple, but I could not get the numerical solution
>> to match the analytical solution.
>> >
>> > <image.png>
>> >
>> > Is there something wrong with the way I set up the boundary condition?
>> I've copied my script below and bolded the boundary condition lines. Thanks
>> for your help.
>> >
>> > Best regards,
>> > Rose
>> > ----------------------------------------------
>> > from fipy import *
>> > L = 1.
>> > nx = 100
>> > dx = L/nx
>> > mesh = Grid1D(nx=nx,dx=dx)
>> > zc = mesh.cellCenters[0]
>> > Dm = 8e-6 # moleculer diffusion in water column
>> > Dmp = 5e-6 # moleculer diffusion in porewater (corrected for tortuosity)
>> > Dbw = 1e-5 # bioturbation in porewater
>> > D = Dmp+Dbw
>> >
>> > Cw=1e-3
>> > Cl=1.
>> >
>> > time = Variable()
>> > beta = 1e-3
>> > #*abs(numerix.sin(w*time))
>> >
>> > # steadystate solution
>> > A = beta*(Cl-Cw)/(D+beta*L)
>> > B = Cl-A*L
>> > Canal = CellVariable(name = "analytical solution steadystate",
>> mesh=mesh, value = A*zc+B)
>> >
>> > C = CellVariable(name="concentration", mesh = mesh,value=0.)
>> > C.faceGrad.constrain(beta*(C.faceValue-Cw)/D, mesh.facesLeft)
>> > C.constrain(Cl,mesh.facesRight)
>> > eqI =TransientTerm()==ImplicitDiffusionTerm(coeff=D)
>> >
>> > viewer = Viewer(vars= (Canal,C))
>> > dt = 864
>> > while time()<86400:
>> > time.setValue(time() + dt)
>> > eqI.solve(var = C,dt=dt)
>> > viewer.plot()
>> >
>> >
>> >
>> >
>> >
>> > _______________________________________________
>> > fipy mailing list
>> > [email protected]
>> > http://www.ctcms.nist.gov/fipy
>> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>> >
>> >
>> >
>> >
>> > --
>> > Daniel Wheeler
>> >
>> > _______________________________________________
>> > fipy mailing list
>> > [email protected]
>> > http://www.ctcms.nist.gov/fipy
>> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>> >
>> >
>> >
>> >
>> > --
>> > Yuan (Rose) Wang
>> > PhD Candidate, Tufts University
>> > Cellphone: 617-699-8006
>> > _______________________________________________
>> > fipy mailing list
>> > [email protected]
>> > http://www.ctcms.nist.gov/fipy
>> > [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>>
>> _______________________________________________
>> fipy mailing list
>> [email protected]
>> http://www.ctcms.nist.gov/fipy
>> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>>
>
>
>
> --
> Yuan (Rose) Wang
> PhD Candidate, Tufts University
> Cellphone: 617-699-8006
>
--
Yuan (Rose) Wang
PhD Candidate, Tufts University
Cellphone: 617-699-8006
function pdex1
m = 0;
x = linspace(0,1,50);
t = linspace(0,86400,100);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);
%Surface plot
figure(1)
surf(x,t,u)
title('Numerical solution computed with 50 mesh points.')
xlabel('Distance x')
ylabel('Time t')
%Solution profile after 1 day
figure(2)
plot(x,u(end,:))
title('Solution at t = 86400')
xlabel('Distance x')
ylabel('u(x,86400)')
%---------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
D = 15e-6
c = 1/D;
f = DuDx;
s = 0;
%----------------------------------------------------
function u0=pdex1ic(x)
u0 = 0;
%----------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
beta=1e-3
Cw=1e-3
Cl = 1.
D = 15e-6
pl = beta*(Cw-ul);
ql = D;
pr = Cl-ur;
qr = 0;
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
