Oseen Flow in Ink Marbling Aubrey Jaffer Abstract Ink marbling refers to techniques for creating intricate designs in colored inks floating on a liquid surface. If the marbling motions are executed slowly, then this layer of inks can be modeled as a two-dimensional incompressible Newtonian fluid. 7 In this highly constrained model many marbling techniques can be exactly represented by 1 0 closed form homeomorphisms. Homeomorphisms can be composed and compute the composite 2 mapping at any resolution. Computing homeomorphisms directly is orders of magnitude faster b than finite-element methods in solving ink marbling flows. e Mostmarblingpatternsinvolvedrawingrakesfromonesideofthetanktotheother;andthese F canbe modeledby closed-formhomeomorphisms. But pictorialdesignsforflowersandanimalsuse 1 short strokes of a single stylus; presented is a closed form velocity field for Oseen fluid flow and its application to creating short stroke marbling homeomorphisms. ] n y Keywords d - ink marbling; Oseen flow; Stokes flow; fluid dynamics u l f . s 1. Introduction c si Marbling originated in Asia as a decorative art more than 800 years ago and spread to Europe in the 1500s y where it was used for endpapers and book covers. h Themathematicalfascinationwithinkmarblingisthatwhilerakingsacrossthetankstretchanddeform p the ink boundaries, they do not break or change the topology of the surface. With mechanical guides, a [ raking can be undone by reversing the motion of the rake to its original position. Raking is thus a physical 1 manifestation of a homeomorphism, a continuous function between topological spaces (in this case between v a topological space and itself) that has a continuous inverse function. 6 0 1 2 0 . 2 0 7 1 : v i X r a 1 2. Dropping Ink First, inks are dropped onto the surface. Consider the tank as being an infinite plane covered with a film of “empty ink” initially. The first ink drop forms a circular spot with area a. If a second drop with area b is put in the center of the first drop, then the total covered area increases from a to a+b. Points near the center will move from small radius to radius a/π; and boundary points will move from radius a/π to radius (a+b)/π. p p The movements of a point on the surface do not depend on its ink color; so the motion of every point p is the same as for the concentric ink-drop case. Given a point P~ and a new ink drop of radius r centered at C~, map the point P~ to: r2 C~ + P~ −C~ 1+ 2 v P~ −C~ (cid:16) (cid:17)u u t (cid:13) (cid:13) Figure 1 shows a pattern formed by serial injection of 75(cid:13)drops o(cid:13)f random size and position. For a more (cid:13) (cid:13) complete discussion of ink-dropping see Mathematical Marbling[1]. Figure 1 Figure 2 3. Line Deformation Immediately adjacent to the line, particles of fluid are moved by distance λ. Because flow is laminar, flow must be uniform along the line. The friction between adjacent lamina results in exponential decrease with (normal) distance from the line. For a line on the y-axis, the displacement in the y-direction would be: −|x| λexp L Because the displacement depends only on distance in the perpendicular direction, the displacements from parallel lines add linearly, which speeds computation of (parallel) raking homeomorphisms. Figure 2 shows a single line deformation through the center of concentric ink circles. A line with unit direction vector M~ and point B~ on the line maps point P~ to: P~ −B~ ×M~ P~ +λM~ exp (cid:13)(cid:16) −L(cid:17) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) 2 4. Short Stroke At the 2016 Lowell Folk Festival Regina and Dan St. John (Chena River Marblers) were kind enough to let me perform an experiment on their equipment. With a marbling pattern already in the tank, I took a rod,insertedit into the tank, movedita shortdistance, andwithdrew itfromthe tank. I couldthenclearly see the effect of a short stroke on the ink contours floating in the tank. The bands perpendicular to the rod motion were compressed in the direction of motion and spread perpendicular to the motion to form a gentle curve. Behind the point where the rod was withdrawn the contours formed a sharp V leading to the 1 extraction point. I repeatedthe experiment,but stopping the motionhalfway,then resumingfor the sametotaldistance; the deformation was indistinguishable from the first stroke. A vertical rod drawn slowly through a layer of inks floating on the surface of a liquid, one of the tech- niques of ink-marbling,canbe treated asthe motionof a circulardisk in anincompressible two-dimensional liquid. Because the movement is slow, viscous forces dominate inertial forces. In the simulation of the two- dimensional Stokes model[2] shown in Figure 3, as the radius of the cylinder shrinks, the y displacements vanish. ThustheStokesflowcannotproducethespreadingIobservedaheadofthemotionoftherod. Inthe simulationoftheOseenmodel[3]showninFigure4the streamlinesareorbitsdisplacingthey coordinatesof points atanydistance fromdisk. Althoughthe asymmetryissubtle, the Stokesstreamlinesaresymmetrical around the y-axis, while the Oseen streamlines are not[4]. Stokes Flow Streamlines and V Field Oseen Flow Streamlines and V Field 2 2 1.5 1.5 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -1 -1.5 -1.5 -2 -2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Figure 3 Figure 4 In “Small Re flows, ǫ = Re ≪ 1” Lagr´ee[5] gives a derivation of the two-dimensional Oseen formula which is composed of a near-field and far-field componenets. The near field is responsible for the flow deflectingaroundthecylinder,whilethefarfieldisresonsibleforthewake. Becausethecylinderiswithdrawn at the end of the stroke, the deflection will collapse, and is not of interest for the marbling deformation. Butthe far-fieldcomponentofthe Oseenapproximationresultsinstreamlineswiththe potentialto be more marbling-like. 1 I was so engaged that I neglected to take a photograph. 3 5. Velocity Field Because the fluid is modeled as incompressible, the divergence of the velocity field ∇·F~ = 0. In polar coordinates: 1∂rF 1∂F ∇·F~(r,θ)= r + θ =0 (1) r ∂r r ∂θ At large r values, the magnitude of the velocity vanishes. lim F~(r,θ) =0 (2) r→∞ (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) The other boundary condition is the velocity at the origin F~(0,0) = [U,0] where U is the speed. Converting this constraint to polar coordinates: F~(0,θ)=F (0,θ)rˆ+F (0,θ)θˆ rˆ=[cosθ,sinθ] θˆ=[−sinθ,cosθ] (3) r θ U =F (0,θ)cosθ−F (0,θ)sinθ 0=F (0,θ)sinθ+F (0,θ)cosθ r θ r θ 2 2 sin θ cos θ U =F (0,θ)cosθ+F (0,θ) U =−F (0,θ) −F (0,θ)sinθ r r θ θ cosθ sinθ F (0,θ)=Ucosθ F (0,θ)=−Usinθ (4) r θ On the basis of the Oseen formula, it is likely that the functions F (r,θ) and F (r,θ) satisfying the r θ constraints are the product of trigonometric and exponential expressions. The system characteristic length L = ν/U, where ν is the dynamic viscosity of the liquid. A solution (5) is incompressible (1) and satisfies boundary conditions (2) and (4). −r r −r F (r,θ)=Ucosθexp F (r,θ)= −1 Usinθexp (5) r θ L L L h i Expressing(5)in Cartesiancoordinates(6)accordingto (3)andintegratingtofind the streamfunction ψ (7): rL−y2 xy r = x2+y2 F =U F =U (6) x y rLexp(r/L) rLexp(r/L) p Uy ψ(x,y)= (7) exp(r/L) 4 Figure 5 Figure 6 Velocity Field and Stream-lines L=1.5 U=0.7 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Figure 7 Figure 8 Figure 5 shows a displacementgraphresulting fromfour evenly spacedapplications ofthe velocity field 2 to a square grid. The deformation matches my description of the short stroke in the introduction; bands are compressed ahead of the stroke and a vee trails it. Figure 6 shows a displacement graph resulting from twoapplicationsofthevelocityfieldwiththesametotaldisplacement. Figure7showsadisplacementgraph resultingfromoneapplicationofthevelocityfieldtothesquaregrid. Clearly,theareaof(deformed)squares has not been preserved in Figure 7. Figure 8 shows the streamlines (contours of constant ψ values) and velocity field vectors from equations (5) in the steady-state. Because the streamlines form closed orbits, it is a rotational flow while the Stokes flow is not. 2 The light gray boxes show the horizontal extent of each stroke; the vertical grid-lines with an initial position to the right of the beginning of the first stroke are magenta. 5 6. Displacement Field So the flow velocity field is not the ultimate goal, rather the displacement after a finite time. Consider the case of a point on the horizontal line where y = 0. Along this line xf > x0. t(xf) is time as a function of distance, the inverse of the desired x (t). f −|x| xf dx F (x,0)=Uexp F (x,0)=0 t(x )= x y f L F (x,0) Zx0 x If x0 ≥0: xf |x|dx L xf x0 t(x )= exp = exp −exp f L U U L L Zx0 h i x0 tU x (t)=Lln exp + (8) f L L (cid:18) (cid:19) If x0 ≤xf ≤0: |x0| |x|dx L |x0| −xf t(x )= exp = exp −exp f L U U L L Z|xf| (cid:20) (cid:21) |x0| tU x (t)=−Lln exp − =−Lln(β) (9) f L L (cid:18) (cid:19) Otherwise x0 ≤0 and xf ≥0: |x0| |x|dx xf |x|dx L |x0| xf t(x )= exp + exp = exp −1+exp −1 f 0 L U 0 L U U L L Z Z (cid:20) (cid:21) |x0| tU x (t)=Lln 2−exp + =Lln(2−β) (10) f L L (cid:18) (cid:19) |x0| tU β =exp − L L x0 ≤ xf ≤ 0 tests the final value of xf which isn’t yet known when trying to compute it. But the transition between x (t) in (9) and x (t) in (10) is β =1; so test β >1 instead. f f When x0 ≥ 0 formula (8) is used. Formula (9) is used when β > 1. Otherwise formula (10) is used. Figure 9 shows the position versus time of five points on the x-axis with L=1.5 and U =0.7. Position x(t) Along x Axis Velocity F(x(t),0) Along x Axis 4 0.7 x =+1 3 0.6 x0=+0 0 x =-1 2 0.5 x0=-2 0 x =-3 n 1 x =+1 y 0.4 0 positio 0 xxx000===+--120 velocit 0.3 0 -1 x0=-3 0.2 -2 0.1 -3 0 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 time time Figure 9 Figure 10 The derivative of xf(t) with respect to t gives the velocity as a function of x0 and t. Figure 10 shows the velocity with L=1.5 and U =0.7. 6 tU/L+exp|x0|/L, if x0 >0; dx (t) f =U −tU/L+exp(|x0|/L), if −tU/L+exp(|x0|/L)>1; (11) dt  ,2+tU/L−exp(|x0|/L), otherwise. Equation (11) can be simplified with r =|x0| and cosθ =x0/r: dxf(t) cosθ tU/L+exp(r/L), if cosθ tU/L+exp(r/L)>1; =U dt 2+tU/L−exp(r/L), otherwise. (cid:30)(cid:26) r tU =U 1+ exp +cosθ −1 L L (cid:30)(cid:26) (cid:12) (cid:12)(cid:27) (cid:12) (cid:12) (cid:12) (cid:12) This one-dimensional case was (pi(cid:12)ece-wise) integrable b(cid:12)ecuase t(xf) is a monotonic function of a single variable. Figure 8 shows that all other streamlines are orbits. Solving for x with a constant ψ value gives a parameterization of the left (−) and right (+) halves of the orbits generated by (7): x=± [Lln(Uy/ψ)]2−y2 y 6=0 (12) q The analogous approach for making F~ and ψ functions of time is to take the time derivative of the functional inverse of the integral of the reciprocal of the velocity as a function of y segmented at x=0. Because the orbits and their velocities are continuous, Uy/ψ along the orbits must be either greater or less than 1. Figure 11 shows that 0<Uy/ψ≤1. x2=[L ln(yU/y )]2-y2 4 L=1.5 U/y =.5 3.5 L=1.5 U/y =1. L=1.5 U/y =2. 3 L=0.7 U/y =.5 L=0.7 U/y =1. 2.5 L=0.7 U/y =2. x 2 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 y Figure 11 Let w(y) be the magnitude of velocity: 2 2 2 w(y) =F (x (y),y) +F (x (y),y) x ψ y ψ [L2+y2]/L2−2y2 L L2ln(yU/C)2+y2 2 =U exp 2 L2l.n(hyUp/C)2+y2/L i (cid:16) p (cid:17) It seems unlikely that dy/w(y) will be integrable and evenless likely that the displacementfield will be expressibleinclosedform. WhilethevorticityofF~ is0wheny =0,itisnon-zeroawayfromthex-axis. This makes sense because ink-marbling develops turbulent eddies when the rod is moved too rapidly. Because this rotational system is dissipative, reversibility of strokes is not guaranteed. 7 Figure 12 Figure 13 Figure 14 Figure 15 Figure 5 shows that with iteration, the velocity field generates a reasonable approximation to the displacement field. It was found that ⌈tU/L⌉ = ⌈λ/L⌉ iterations produce deformations that are visually acceptable. Figure 12 shows a marbling deformation having stroke-length λ = 1 (half the height of the image) and L=0.15. With homeomorphisms we can render either by filling contours computed from the boundaries of ink- drops, or by finding the color mapped from each point of the display raster (using the inverse homeo- morphism). For the exact homeomorphisms (drops, lines) the geometries rendered are the same. But the imperfect reversibility of the stroke approximation leads to geometric artifacts. Executing a stroke from B~ to E~ followed by a stroke from E~ to B~ does not completely undo the first stroke’s deformation as seen in Figure 13. Figure 14 shows a pattern with five radial strokes. The pattern is asymmetrical because the strokes were applied in the sequence 0,2,4,1,3 and each stroke affects the whole space. Figure 15 shows the raster- renderingofthe same five-strokepattern. Unlike Figure 14,the peaks ofthe fourthwhite bandreachnearly all the way in. 8 7. Reversibility Averaging the velocity field and the negative of the velocity field with the stroke beginning [0,0] and end [tU,0] points reversed results in a velocity field which is more nearly reversible: r = (x−x )2+y2 s= (x−x )2+y2 tU =x −x B E E B p p U y2 −r y2 −s F = 1− exp + 1− exp x 2 rL L sL L (cid:20)(cid:18) (cid:19) (cid:18) (cid:19) (cid:21) (13) U (x−x )y −r (x−x )y −s B E F = exp + exp y 2 rL L sL L (cid:20) (cid:21) Figure16showsadisplacementgraphresultingfromfourevenlyspacedapplicationsofthenewvelocity field to a square grid. Compared wth Figure 5, the point of the vee is rounded and the tail has slightly more displacement. A more roundedpoint is not inconsistent with the effect of a stylus having finite width. Figures 17 and 18 show displacement graphs resulting from two and one applications respectively of the velocity field with the same total displacement; both are noticeably distorted. Figure 16 Figure 5 Figure 17 Figure 18 9 8. Application This section generalizes equations (6) and (13) to work at any angle and any stroke length larger than 0. Vector function Q~ maps point P~ to its new position as a result of a stroke from B~ to E~. For (6) the transform is: 2 rL−y −r xy −r Q~ P~,B~,E~,L =P~ +λ exp , exp rL L rL L (cid:16) (cid:17) (cid:20) (cid:21) λ= E~ −B~ (cid:13) (cid:13) N~ =(cid:13)E~ −B~(cid:13) λ (cid:13) (cid:13) (cid:16) (cid:17). C~ = B~ +E~ 2 (cid:16) (cid:17). r = P~ −C~ (cid:13) (cid:13) x=(cid:13)N~ · P~ −(cid:13) C~ (cid:13) (cid:13) (cid:16) (cid:17) y =N~ × P~ −C~ (cid:16) (cid:17) For (13) the transform is: λ rL−y2 −r sL−y2 −s x y −r x y −s Q~ P~,B~,E~,L =P~ + exp + exp , B exp + E exp 2 rL L sL L rL L sL L (cid:16) (cid:17) (cid:20) (cid:21) λ= E~ −B~ (cid:13) (cid:13) N~ =(cid:13)E~ −B~(cid:13) λ (cid:13) (cid:13) (cid:16) (cid:17). r = P~ −B~ (cid:13) (cid:13) s=(cid:13)P~ −E~(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) x =N(cid:13)~ · P~ −(cid:13) B~ B (cid:13) (cid:13) (cid:16) (cid:17) x =N~ · P~ −E~ E (cid:16) (cid:17) y =N~ × P~ −B~ (cid:16) (cid:17) The number of segments for the computation is n = E~ −B~ L and the increment vector is I~ = E~ −B~ n. l(cid:13)(cid:13) (cid:13)(cid:13). m (cid:13) (cid:13) (cid:16) (cid:17). P~ ←Q~ P~,B~,B~ +I~,L (cid:16) (cid:17) P~ ←Q~ P~,B~ +I~,B~ +2I~,L (14) . (cid:16) (cid:17) . . P~ ←Q~ P~,B~ +(n−1)I~,B~ +nI~,L (cid:16) (cid:17) Figure 19 shows the reversible short stroke; the point is not as sharp as Figure 12. Figure 20 shows the reversible short stroke followed by the reverse stroke; cancellation, while better than Figure 13, is not complete. Figure 21 shows the contour-filling and Figure 22 shows the raster-rendering versions of the five radial stroke pattern. The difference between them is less than between Figures 14 and 15. 10