From 777d75baa557ee554f59d1284eb2c91f9003307c Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Wed, 20 Apr 2022 09:34:32 +0200 Subject: [PATCH] Checkout files from age --- doc/ADI_scheme.org | 94 +++++++++++++++++++++++++++++++++++++++++++++ doc/grid_pqc.pdf | Bin 0 -> 5809 bytes 2 files changed, 94 insertions(+) create mode 100644 doc/grid_pqc.pdf diff --git a/doc/ADI_scheme.org b/doc/ADI_scheme.org index 6a86caf..8a3fd92 100644 --- a/doc/ADI_scheme.org +++ b/doc/ADI_scheme.org @@ -1,6 +1,7 @@ #+TITLE: Adi 2D Scheme #+LaTeX_CLASS_OPTIONS: [a4paper,10pt] #+LATEX_HEADER: \usepackage{fullpage} +#+LATEX_HEADER: \usepackage{amsmath} #+OPTIONS: toc:nil * Input @@ -104,3 +105,96 @@ bc(i,j).\text{value} & \text{if } bc(i,N-1) = \text{constant} \\ \end{cases}$ [fn:1] $p$ is called =t0_c= inside code + +** Finite differences with nodes as cells' centres + +*** The explicit FTCS scheme as in PHREEQC + +The 1D diffusion equation is: + +#+NAME: eqn:1 +\begin{align} +\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha \frac{\partial C }{\partial x} \right) \nonumber \\ + & = \alpha \frac{\partial^2 C}{\partial x^2} +\end{align} + +We discretize it following a Forward Time, Centered Space finite +difference scheme where the nodes correspond to the centers of a grid +such as: + +[[./grid_pqc.pdf]] + +The left boundary is defined on $x=0$ while the center of the first +cell is in $x=dx/2$, with $dx=L/n$. + +We discretize [[eqn:1]] as following, for each index i in 1, \dots, n-1 +and assuming constant $\alpha$: + +#+NAME: eqn:2 +\begin{equation}\displaystyle + \frac{C_i^{j+1} -C_i^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{i+1}-C^{j+1}_{i}}{\Delta x}-\frac{C^{j+1}_{i}-C^{j+1}_{i-1}}{\Delta x}}{\Delta x} +\end{equation} + +In practice, we evaluate the first derivatives of $C$ w.r.t. $x$ on +the boundaries of each cell (i.e., $(C_{i+1}-C_i)/\Delta x$ on the +right boundary of the i-th cell and $(C_{i}-C_{i-1})/\Delta x$ on its +left cell boundary) and then repeat the differentiation to get the +second derivative of $C$ on the the cell centre $i$. + +This discretization works for all internal cells, but not for the +boundaries. To properly treat them, we need to account for the +discrepancy in the discretization. + +For the first (left) cell, whose center is at $x=dx/2$, we can +evaluate the left gradient with the left boundary using such distance, +calling $l$ the numerical value of a constant boundary condition: + +#+NAME: eqn:3 +\begin{equation}\displaystyle +\frac{C_0^{j+1} -C_0^{j}}{\Delta t} = \alpha\frac{\frac{C^{j+1}_{1}-C^{j+1}_{0}}{\Delta x}- +\frac{C^{j+1}_{0}-l}{\frac{\Delta x}{2}}}{\Delta x} +\end{equation} + +This expression, once developed, yields: + +#+NAME: eqn:4 +\begin{align}\displaystyle +C_0^{j+1} & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}-C^{j+1}_{0}- 2 C^{j+1}_{0}+2l \right) \nonumber \\ + & = C_0^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( C^{j+1}_{1}- 3 C^{j+1}_{0} +2l \right) +\end{align} + +Which leads to the following term applicable to our approach: + +#+NAME: eqn:5 +\begin{equation}\displaystyle + -C^j_0} = \left[\frac{\alpha \cdot \Delta t}{\Delta x^2}\right] \cdot \left(C^{j+1}_1 + 2l \right) + \left[ -1 - 3 \cdot \frac{\alpha \cdot \Delta t}{\Delta x^2} \right] \cdot C^{j+1}_0 +\end{equation} + +In case of constant right boundary, the finite difference of point +$C_n$ - calling $r$ the right boundary value - is: + +#+NAME: eqn:6 +\begin{equation}\displaystyle +\frac{C_n^{j+1} -C_n^j}{\Delta t} = \alpha\frac{\frac{r - C^j_{n}}{\frac{\Delta x}{2}}- +\frac{C^j_{n}-C^j_{n-1}}{\Delta x}}{\Delta x} +\end{equation} + +Which, developed, gives +#+NAME: eqn:7 +\begin{align}\displaystyle +C_n^{j+1} & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 2 C^j_{n} -C^j_{n} + C^j_{n-1} \right) \nonumber \\ + & = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot \left( 2 r - 3 C^j_{n} + C^j_{n-1} \right) +\end{align} + +If on the right boundary we have closed or Neumann condition, the left derivative in eq. [[eqn:5]] +becomes zero and we are left with: + + +#+NAME: eqn:8 +\begin{equation}\displaystyle +C_n^{j+1} = C_n^{j} + \frac{\alpha \cdot \Delta t}{\Delta x^2} \cdot (C^j_{n-1} - C^j_n) +\end{equation} + + + +A similar treatment can be applied to the BTCS implicit scheme. diff --git a/doc/grid_pqc.pdf b/doc/grid_pqc.pdf new file mode 100644 index 0000000000000000000000000000000000000000..b970d186492f1e3f6c90268bc2a56447fc4cabe8 GIT binary patch literal 5809 zcma)A2{@Ep*p?d0WG_2;?IFw=mSt47D>`}JJ zk|k@1lBHyee}?M&zU%Az|9{?@_q^vh&pGEk&vQR>UH2_$bVgGgDj^FNto&TE4VDHV z0M6ke7=ZvJ^-)+Cf-3+gnoPg|07z=0-3Ta8V(W$^pw6HiaZV_(k`fqC@I)aoV4t*W zrUr4xpfu0tEH;MZV=gI-Q%Y!4p+S`P*mGgGKUW0~A4Umot$9eh{kFFI76#mySQNMj<0=Z7Si1tmoY-}345?G@#Fm~24nYRT zF(A&1F>r@Z)N5LIGQ8Cfn~O$f6gI?O`b1vwKuAzdO#7L1paWcWv6WBEeq9>4_d!Z< z#_YC9l%i09Cg{;pN_0`PZauSml1j{lyFACjq{{7JVM4U{^(j}j@}`;D-rGlPJ;l`T z&6YtrVv`#3G@o{sJ{}o0D{W1!PU|OdA!wa1mb06_6}s=oCUWG)xL_KRLf|;)Z5H~0 z@%M@^_Xs!Qa$aY$GQNH|yyk$aD7oqR@<}tNpgjMFwjA#V~o=}`6Zo?mZ(qg?9P?1-tb$J)11onb_B_e!o@<_JA~c#i@duPEom zv7AlH&Y+RIC3DZRz~!M#j$J1%1PM(o#;N&j?V8+9fzR%h`R8X|$gW~9j`{);xeexD zs=9w|kTN|F>udLhP9w+R@B%kqIsdx*VzyE5W$Oh!&n|Hh#f#uT@`==l8#8?oro9W zz`u}MOfZ5;rCsO5p&6G3* zp!)-wXiq!=kpB7nzoS>N^EtQ?avU`A{HzDdoKWTIFgi?@9F@sfGjDS=K}XlWx$e1+ zKRqgF@kSdH~Bu*ivb|L;gAqymwVnE-<7onMW&{pO4kj>xT{<->yMUt zcDEDZ$S?+xoSl)j^Ya}2x}OIg=TQ^z2>XE~g%hg&k4jEAr?xaZ z$Eh^_<4$c~U=3ZPeup5Y8)N``olEX*L&nU=pn6D^O@`5m!r9)Ef`Y<{$L!2IRfdjN z;;^K*uWT7B+gR}J9aZ>`Q{w?|cDJq+KDbvmS34&^H~K|kZo8}Y`!zejT${{(D7YuG zzQ!`HEt0I$H7+y0W;u|J+&q+lA*0|~)N&kosk0BA?-em*RN3@{nT%n<*ida}nol?_ z`URXLMYMRzNtyC#bssIxDW{M#Z_+R9vR%Z~&)pNv1-v)yT$N;_dT_f_xmo-*ZO_D6={ zz#~2;KEYErZ)!!RSsvq*3LHuZa$-nPpH1*>zR@9kv`wu%vgR25lnzUxpp(c7-|I+O zE5oaik}LdTk;>G9NGEY>Va7A+d^)nV`X%v2GyJR2kGvW2mu5^?**DX!s+J^5w$pkh zj`CHvTg=dHrj80epAqJZSMto*iP)j}QQbrq%&e^fKGinC7#{*3swXQ83xZHlG%>J& zn^X>YLAX9O9J-YmtkuBywm#@sOgI$6l7{4piI_u@>xZ)1vvh>ZLs(~nI6ERhWehAq z3O<3tr@{rxjz6lyEC-?{rbITx`2+2u`(NJr5&bnZQW%y_qZsNq%qvH+Vmk;ky9ru|J%TNha zs-fQt%3yCKYl5^{LYWrFImbq}AF8>I0g4$`w-e_07OkcGZsJhi-wX;sO~x6;?t?UL zyHyd+U-0ucyc-i|zt)wU?0pwC5P*I7Mo?ohL2zXG^+rCQupqyzUZqemf&S%-Zs*xe zp|FKnfztwy^0I?R9s0W*t^vYw6`iWaxkYB=3Yzg${GmuzAtS`-g-UNq**LgM-vc_` z(7E@fY%~Ra>8F*T=Sz(Z&U=Q@Wt(qH61a!wp5EGyzB>B#K|x+hIVNCbm)@B{T5J9~ z2giW1(cZ{Iz6H*g3pHD<3k%`lm-{)o*1xcdCilhSc0ZQ9NR5xT{#I6vwauATrR@6N z15$Jyv+Iz*(7XJ=gT>3UuH#n1>gZ9i&_I&!@&Uaef`kZ)od zpG!-?DQbx(+LgPBu)}PiHI8HL(vT}Du#He7hf6)OYXia)1o|6$G#}^R7esYDQ^A;} zo{b6-vgPw<4%@8ND70-sAF{C%->{W5oxzg6|@G-wR?DGFsYL< z^QHn@2$3fpNzrQ#@G0_c>p+=z#>&ysh`QBIcSAPKZ3!)c!Ay>sO!d9~lvtgc#MxO+YY_g5DL407*_)gys-lGGcdf6Ej|(m3E!P0d zLRv=_Xp!&e+P_jCEqcRz3wP>N>Ch|Nyu+nO4Vmzo7xGk7?AySwLQI^|UHXDm^@+Pa z3K|(`I$Q!+wF38RP4Q@Sw9finSo(07Lg-`^ z8%OVmNNBT-qxX!Nz+E+`@4{vBmo7gkPR5`Q*bQYa;}c;Q_1rO!-f zu~lZ;HOi)*x`_T!3I7pL=3pr=j}AAp2&iG{bOcAKaA_@@R`Bllyz^fN|_AfYo0@ny~kfXN^He9>NnUXyb`mr#(-Wt4lT)#uX?Rq zuiX zq{v38*`9w#i(N~z(>>Qh_>~Vamd{Y{#3p$r_v_m5=y0)1BOHL8%ICti;8ny0C4R)nq&&QgaT!B?`43v}h zT*2xdY@$7M5%(2O0}+NBY_at>H{_Ip&Q6(bNGpYOPvPajTGMP?(Oq(MFQ?b!4hc^k zm7}(6Zr@<4i!pTp&EO0SkVJ2(U3bI_PU2N`9A$tHcyF4fX($o zY>`C}ayZ*{W40Smw&+C&-Jl}^TYp5qSwLS?JkMTq!6*)WCzhh4UcaN!P1=EOHk*14 z!qaIP7*@j=UB;$YcKmf2msn45QVq52G|R)rV#^yN?g|6;;sZ!U8GC*iB&=$J&aFje z^g3?x#Ax)kyYh+p>J^|aBx2bX5u>$y(rfC94^7-O^;I6_dioWCdyRcQ6x~zuh&a|| zKCel^Ppz+gK&4Y!8_@L$VINxMmhNTz4teQxvd>}Ml=!n7c9YyY9ABc&>}5MV>DKK; zB)xnU&d(u{HP+*XkhuRc=OfQNr$Z_V!{~4Q-0XcFb(U%x&i>Y-Q*Ne`?9K9O2``g! zYjHxKanw>Tw?bS&EG90^sLRT(ut-ZOrrLs_?>KmAO|}*9p^=a<-xloJf5*c1<&%Dj zbeidG$HA8m7RoEQ2414(nUXwmpNy?sWG$U9v;WYmw_7dgHMpddkHrnY9W6mDZEqDh zS7x5Pb*epDyu{4B^6>z3LB&WF+nSkzx79*w&22xu#X~6dSs6^ldB|50=5X#Yc3e@N zs}8kCM&0@NnD4GWcLqlcK}^0G&dT(OLsBoy^Tzm2Jw3LF#}nwc z&cTCeTIAW*!kc(CsY_Sk2FhYj?Un%L;L~Ex8==j9 z`K?cEvp#u$xmP=EUP3QA37!^6&F;M3UpQQBUB!Lxk$3MZYK1o_tvX=lVRBON)GAvP z?sK8rMbkmdh{^cSrV4_2$T&pfw1Q+o$?CV6wAm zNya?u+x6;+te|Fzucm$cZqB=@Z|{6<(V;V@-xy9yw$AsvfAX5(V70kfK+4Ax6joBS z)n!YdkGYyqM95ghQ;$B0{@v7g_O!Az_>;Vk2q}#WxZ7SV=S(7ph1U{p+`mjC-9+@= z^u{F#!PXu9%I%y;BZF*xm$tSxBj>v@o0{JgUYxDP>p5PzTzjFyirdyIZ#-wX&UdX! zW^XtQd%;*NyLO)WBd@&VYoWa`rqx?VBXRN-aW!RxXIYn@*w z>|rYd9JqJJL^rnHzo6ZdtdRaj|8Ewt&kcUk2pB{P`io*grT)MpNF?w7$Np^-y0IY; zn$|DAV4kPypiuJS(D=NiMr43=6=h)aV{QI3yo$@~_^HA~FS5Ohhds zLx2qM_2w9Jk|^6&@=ja84(v|Ew5p4yQb@-Nd>lJ zk^V@q($xl2+VgETINPv~V!^}h2*(fHSzRo<=H<`n`Fc9JXx$(=m#*v4P$|$Wv96+e z5T3EdN#c5ux)V72Hrfe%OWAW*J8iM;cx0_m9Q;#O)Ct3fzF-@uMUd44n$wi)cq;#> zplftNRxcvA*cV;0lAkBN;oM86*-KjU770}y6K2w6>00+-S(xSGR)0q3v5Q19k42drr znCqMBiy5LZ4wvyJI1CaiZh~^T?1ucyQ4@(lyZHhq4}5@qCy|3?+>kE#UqqU?q?#Je z2e1ajp)zuSxRe|efI=V;;tdW1Z1*SV5Rh(YM^&tg8w!BzTU7CmC?egLk%R4bN$*MO zi^HI@`v!HSyEY2#;!47+zri-HlI_YwM+Qb5uVNa!XB z4}XQ_V7-4(-hGoA5|7%KA^y+!kfa24SEMH~2m3Sikw5!Tse>2%XY>ImDL)6Lqk$!A zWoWDmAZdZds$%izUxxjKEYMB_S3EJAQgG5){`-=Wl>^|!F#r1f$rCO^+V;OMnf>bj z?{OH6gtmh_X&=%>YNe&6f48L3U++MACz8*v_F%@r{J(4@?|oblQ-V8(MU$jMVE_Nd z({Q>?M&_|Rd$yKkz{m`&$|DZC~vTWd}WvOkjY}6fV@kZ z`TL<<(+m7dqBnlC0<;1t6ARF3f2Xw4rz<_PEeL;e5KzKHb? zm4W@P|6gj*UqH7sK{+2_^FX?lmD}(ADJJ_05ZaHJ1^ELI0@hGK3XnReUK1RSScboq zI7uC>GY)_rD8YY^rlY<6#$Lwtpv0h{>p|&3K{wkuv?0qOmW3*6^GuPI