中国科研软件论坛 » 【FlexPDE讨论专区】 » question about volume flow

查看完整版本: question about volume flow

michael 2007-9-4 11:03

question about volume flow

[size=2]I've built a 3D diffusion model in which substance enters the model from a finite source (e.g. a droplet). This means, the particle source has got a starting volume which decreases over time, and if the volume reaches 0, the flow of particles will stop . The amount of substance entering the model in one timestep depends on the volume which enters the model out of the source. I've realised the decreasing volume of the particle source by implementing the global variable vol into the system. The decrease of volume is being calculated by the law of Hagen-Poiseuille (Volume flow V out of the droplet through a round pore) and realised in the EQUATIONS section by dt(vol)= V. I've monitored the development of vol in a history plot, and the curve is not as smooth as it should be. Also, in the system there are exp-functions implemented in whose exponents the variable vol is being used, e.g. exp((-vol)/6.13249e-16). At a point where vol comes close to 0, FlexPDE stops the simulation with an 'exp argument too big'- error. I've worked around this difficulty by using an IF THEN ELSE construct, but I fear that something in my system is not working properly. Also, I've tried different parameters (s. .pde-file), and sometimes the exp-error occurs, and sometimes it doesn't.
Another question: Is it possible to take the inflow of volume from the source into the system into account when calculating the particle concentration? Considering an incompressible fluid, the volume entering the system pushes the particles in the system further away from the inlet. Is there a formula to calculate this?
I hope you can help me out! [/size]
[size=2][/size]
[size=2][/size]
[code] title
  'Diffusion and convection of ACh'!K9Y#M v"PlR ^Pu
H3Q8i0g%Ji-n {9|
coordinates cartesian3
7EWBp[4q#Z[&y
variables
  u(threshold=1)  !1e-3
2B"hTD/y9{iEP.@
global variables
  vol(threshold=1e-18)  !droplet volume
-T-iU~ N

definitionsA!A1@b8rJy)ZB
D = 12.2e-10             { Diffusionskoeffizient [m^2/s]  ACh}

{ System variables }(k&Dv| ICq
  r_ben = 3e-5  !Inlet area radius [m];~9ULa'q5E`
  height = 1.8e-3 !cylinder height [m]
  r = 9e-4  !radius of cylinder [m]!L;SGPz-k9a
a!V%Hm+f3H
{ Membrane parameters }  !here, different parameters have been used
uC0p&q#t~@*t
!diameter of pores!k$gx-ALG[yx~;v
d_pore = 2e-7  !1e-7, 5e-8, 1e-6iRf1H"WFs/wJ
!thickness of membrane [m]TDF8D8e
d_mem = 10e-6}j8_!W@ @$l6mg)v
!freie Fl鋍he der Membran
ff = 9.43/100 ! 5.54, 1.54, 15.7
N8B`L+_
!Stoffmenge, die insgesamt im Tropfen zur Verf黦ung steht [mol] = Ausgangskonz. [mol/m^3] * Start-Tropfenvolumen [m^3]
v0 = 26e-15  !starting volume of droplet
a_drop = 1 *  v0 !concentration of substance in droplet (1mMol)/\3] ]p'R"T

!Fluss aus der Quelle:7\_B.J.C5B
!J =  (D*(((a_drop-integral(u))/vol)-val(u,0,0,1e-6))/d_mem)*ff  { Einheit: [mol/m^2*s] }
J =  IF  vol <= 20e-20 THEN 0 ELSE (D*(((a_drop-integral(u))/vol)-val(u,0,0,1e-6))/d_mem)*ff  !SWAGE(vol-(2e-14), 0, (D*(((a_drop-integral(u))/vol)-val(u,0,0,1e-6))/d_mem)*ff, 2e-3)  { Einheit: [mol/m^2*s] }
(Mgzr:DR3p
!************************Oberfl鋍henspannung, Konvektiont)m[Iai$n|2c;i?
!z7m1QxKA8`
! surface tension(H2O) = 0,073 N/m\QO/h2B:sw
s_drop = 0.0734cqNXud(F/Yu
!dynamic viscosity = 0,001 Pa*s = N*s/m^2+T;Z7c0y%yW
dV = 0.001] X6Aiq#_E(b
  !number of pores in inlet area
N_pore = 8.49e3A*X%]:IY7p3F
%?jd OR$tV]
!Calculation of contact angle teta and drop radius r_drop via v_drop and r_ben
!Formel aus Origin: x = Tropfenvolumen v_drop, y = Tropfenradius r_drop
_{7@3zZ0RLO
a1 = IF vol <= 20e-20 THEN 0 ELSE (-vol)/3.3553e-15
a2 = IF vol <= 20e-20 THEN 0 ELSE (-vol)/6.13249e-16 \5Adg|7}

a3 = (3.73568e-4)*exp(a1) + 0.00162*exp(a2)
r_drop =  a3 + 4.74773e-5
2e&jqS[1y{
!r_drop = 4.243e-5
nJx9PT{x4x G%B
!contact angle: Kontaktwinkel = Zentriwinkel/2 + Umrechnung ins Bogenma?(->rad) f黵 cos(x) ABER: arccos(x) braucht x in Grad(?!: eAD6oy&?l `UJI
teta = (180/pi)*arcsin((2*r_ben)/(2*r_drop))n2` h&V&`%^q M$q _-M
!teta = (180/pi)*arcsin((2*r_ben)/(2*4.243e-5))

!droplet pressure: Berechnung des Drucks innerhalb des Tropfens
pcap = (2*s_drop/r_drop)*cos(teta*(pi/180))
!pcap = (2*s_drop/4.243e-5)*cos(teta*(pi/180))o:Y0Mz J C(ND k3P#B

`(M!RPI

!Gesetz von Hagen-Poiseuille (ohne Kinetik des Aufpralls des Tropfens)N8Ey{C-Z
V_hp =  -((pi*(d_pore/2)^4)/(8*dV))*((pcap*N_pore)/d_mem)#H5wc(Ehe
!Fluss [mol/s] = Konz. im Tropfen [mol/m^3]*V[m^3/s]  ->Fluss durch ALLE Poren (ges. Ben.fl鋍he) ->/Ben.fl鋍he: [mol/m^2*s]
!J_hp = (((a_drop-integral(u))/vol)*(-V_hp))/(pi*r_ben^2)#yk([_2d#Qm r
J_hp = IF  vol <= 20e-20 THEN 0 ELSE (((a_drop-integral(u))/vol)*(-V_hp))/(pi*r_ben^2)  !SWAGE(vol-(2e-14), 0, (((a_drop-integral(u))/vol)*(-V_hp))/(pi*r_ben^2), 2e-3)M{*g;P0w!p.G

J_ges =  J + J_hp  !SWAGE(-vol, J+J_hp, 0, 2e-3)
!**************************************9OiY%o.{"Y

!Zeitschritt der Simulation [s]
timesteps1 = 0.0005
timesteps2 = 0.2
timesteps3 = 5e-7GO4s(j9oK@;A
timesteps4 = 0.057W'H:z/D(f
H)ph.{d
initial values  { Anfangswerte der Variablen }Bc p`|&aN
  u = 0
  vol = 26e-15vyz.v]8sHV@
'c2q+{n6m9I
equations
u: div(D*grad(u)) = dt(u) { Diffusionsgleichung }uz,[qJ9]jy
vol: dt(vol) = v_hp  !膎derung des Tropfenvolumens 黚er die Zeit: Vdrop(n+1) = Vdrop(n) - Volumenstrom(n+1)6N~/DGR
F e|G \"[(\#P2Sp
extrusion  { 3D-Struktur }X"b'{}:mc As

surface 'fluid_bottom' z = 0^a(LDb u
layer 'fluid' RYV~4M/F_
surface 'cylinder_top' z = height

{ Erstellung des Gef溥es und der Quelle }
boundariesa0B'su;o4|6F#d%I
region "cylinder"  { Rundes Gef溥 }.i|[]T:C{2[ s
start(0,r)
arc(center=0,0) angle 360 nYG6InG"g*A
E'o6Z#E D2R
limited region "bottom_pore"  { Beschr鋘kung von 'bottom_pore' auf 'fluid_bottom'. }
  surface 'fluid_bottom'
  natural(u) = J_ges !J+J_hp  { Randbed.: Flux an 'fluid_bottom' in der Region 'bottom_pore' = D*grad(u) }

start(0,r_ben) { runde Quelle, Radius = Rad. der Benetzungsfl鋍he }4L#z;j2M$s.Ni!A g5^9Bq
arc(center=0,0) angle 360^U[x3V| ^yp

time from 0 to 3e-6 by timesteps3 to 0.025 by timesteps1 to 1.5 by timesteps4 to 30 by timesteps2*]/B_;_%?zbi0s(a'xi%e
yCZt7ysW
plots  { Ausgabe der Schaubilder }5f/@S"SH
  for t = 0 by 0.06 to 0.3 by 0.1 to 1.5 by 0.2 to endtime
contour(u) as 'Konzentration ACh' on y = 0  painted nominmax notags noheader zoom (0,0,4e-4,7.5e-4) !points = (14,70) export format '#t#x#z#1' file = 'Contour_1mmol_1祄_IFTHEN_20e-18_diff+konv_dz=5祄.txt' merge noheader  !range (0.01e-2,3.5e-2) PNG
contour(u) as 'Konzentration ACh 100nm-Membran' on y = 0 painted nominmax notags noheader
elevation(u) from (0,0,5e-6) to (7.5e-4,0,5e-6) points = (150,0,0) export format '#t#x#1' file = 'ACh_5祄_1mmol_200nm_flow+kin_dx=5祄_02.txt' merge noheaderXik N&S3Q9R4P
elevation(u) from (0,0,0) to (7.5e-4,0,0) points = (150,0,0) export format '#t#x#1' file = 'ACh_0祄_1mmol_200nm_flow+kin_dx=5祄_02.txt' merge noheader
elevation(u) from (0,0,7e-6) to (7.5e-4,0,7e-6) points = (150,0,0) export format '#t#x#1' file = 'ACh_7祄_1mmol_200nm_flow+kin_dx=5祄_02.txt' merge noheader4l X~1H?V
elevation(u) from (0,0,1e-5) to (7.5e-4,0,1e-5) points = (150,0,0) export format '#t#x#1' file = 'ACh_10祄_1mmol_200nm_flow+kin_dx=5祄_02.txt' merge noheader
elevation(u) from (0,0,15e-6) to (7.5e-4,0,15e-6) points = (150,0,0) export format '#t#x#1' file = 'ACh_15祄_1mmol_200nm_flow+kin_dx=5祄_02.txt' merge noheader
6Gi:g"zc6R4W:J%O
Summary export file='Summary_ACh_1mmol_200nm_diff+konv_dx=5祄_02.txt' merge noheader
  report(t) as 'Zeitpunkt'
  report(J) as 'Fluss aufgrund Diffusion'
  report(J_hp) as 'Fluss aufgrund Konvektion'
  report(J + J_hp) as 'total input flux (j + j_hp)qr3@U\7uX'a
  report(J_ges) as 'J_ges',u:m0Fn/DdW4Q
  report(V_hp) as 'Volumenstrom aus Tropfen (Konvektion)'
  report(integral(u)) as 'Teilchenmenge im Zylinder'3Y-R!^+}C0K D,qK
  report(a_drop-integral(u)) as 'Teilchenmenge im Tropfen'
sm*f"d!`7y g N
report(vol) as 'Volumen des Tropfens (vol)'xj0x,[ Z'g:O)o

  report(r_drop) as 'Tropfenradius'
  report(teta) as 'Kontaktwinkel'
  report(pcap) as 'Kapillardruck'!x1n3FQ)?

  report ((a_drop-integral(u))/vol) as 'Konzentration im Tropfen (/vol)'
2? G-? w%`%I'mf
  report(Tintegral(surf_integral(J,'Fluid_Bottom','Bottom_Pore'))) as 'total J (surface integral,tintegral)'/F5Y2iYG[.}E`m
  report(Tintegral(J*pi*r_ben*r_ben)) as 'total J*pi*r_ben^2(tintegral)'
report(surf_integral(J,'Fluid_Bottom','Bottom_Pore')) as 'total J (surface integral)'
report(J*pi*r_ben*r_ben) as 'total J*pi*r_ben^2'L`Xj3T E\

  report(Tintegral(surf_integral(J_hp,'Fluid_Bottom','Bottom_Pore'))) as 'total J_hp (surface integral,tintegral)'
  report(Tintegral(J_hp*pi*r_ben*r_ben)) as 'total J_hp*pi*r_ben^2(tintegral)'
  report(surf_integral(J_hp,'Fluid_Bottom','Bottom_Pore')) as 'total J_hp (surface integral)'
  report(J_hp*pi*r_ben*r_ben) as 'total J_hp*pi*r_ben^2'0t\VpX{
2q7r~e-P/l
History(a_drop-integral(u))mO/o$s_Jue
history(pcap)
history(teta)
History(V_hp)
History(vol)
history(j)
history(j_hp)OOW"[0Zg
history(J_ges)2x9k6Tn7kp1m:Z
!history(SWAGE(vol, 0, J+J_hp, 2e-3))
m&D G6~l3L_!K
! table(u) zoom(0,d_mem,50,50)  points=50 merge noheader format "#t#x#b#y#b#1" file="table_dmem_.tbl"Bu%k `~5B
;Tg|%QR1JG
|u2`^{]
  {report(integral(u,"drop_cyl","fluid")) as 'Volumenintegral im zylinder 黚er der Benetzungsfl鋍he'
  report(Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom","Bottom_Pore") )) as 'total input flux(dz)'u3PC.z(Ux
report(Tintegral(Surf_Integral(-D*dz(U),"Fluid_Bottom"))) as 'tif(dz, Fluid_bottom)
g/Ss3@Y.c
  report(surf_integral(-d*dz(u),'Fluid_Bottom','Bottom_Pore')) as 'total dz(u)(surface integral)
  report(Tintegral(surf_integral(-d*dz(u),'Fluid_Bottom','Bottom_Pore'))) as 'total dz(u)(surface integral,tintegral),~2z^hL#K
  report(val(-D*dz(u),0,0,0)) as 'Normal flux at (0,0,0)'l9k#wx:pa
  report(Tintegral(val(-D*dz(u),0,0,0))) as 'time integral of normal flux'3QE#X&R iJh8H
  report(val(-D*dz(u),2e-5,0,0)) as 'Normal flux at (2e-5,0,0)'$P)})G n(\,n

!TRANSFER: schreiben des Meshes in Datei
  transfer(u) file = "eFQ_mesh.dat"]yBdI W?
}8Rx5f,Y1g$D
Fd(T7EQ

end?VYJYTp-_2zS
z.HSD.k1eLWK
8v#fq7NSRgP

[/code]

michael 2007-9-4 11:04

[size=2]1.
You have set the Threshold for U at 1, while the actual values of U are significantly smaller than this. THRESHOLD is "the smallest variation for which you wish FlexPDE to maintain accuracy". In declaring a threshold of 1, you have essentially allowed 100 percent oscillations in some areas. I believe this is what is causing the ragged behavior of Vol. Restore the Threshold of 0.001 that you have commented out, or make it even smaller, so that accuracy will be maintained in the small wings of the U distribution. R"[-}+uv`

2. t_3I+SCRN4H:m
Your input flux is clamped at zero when Vol falls below 20e-20, but the time rate of change of Vol is not. So vol continues to go negative. This may not matter, since the flux is zero, but it looks funny on the plot.

3.
In principle, you could model the pressure increase and fluid motion due to incoming material, but my guess is the diffusion dominates the process. [/size]

页: [1]

查看完整版本: question about volume flow