Commit 37d99ac4bcbcf7b78a0d203b2f7a4d782d5041b3

wlsqm_bchandling.lyx: small update vol 3
docs/wlsqm_bchandling.lyx
(833 / 248)
  
11891189\end_layout
11901190
11911191\begin_layout Itemize
1192for Euler flows,
1193\emph on
1194a priori
1195\emph default
1196 there is not much difference between stationary obstacles and axially moving
1197 obstacles, because the tangential obstacle BC is of the perfect-slip type.
1198\end_layout
1199
1200\begin_layout Itemize
1201but also
1202\emph on
1203a priori
1204\emph default
1205, from the viewpoint of the structure model, the addition of the nonlinear
1206 term into the flow model (compared with potential flow) may change the
1207 pressure field drastically.
1208\end_layout
1209
1210\begin_deeper
1211\begin_layout Itemize
1212indeed, the velocity fields of Euler flows typically look completely different
1213 from those of potential flows.
1214 Although the Euler model neglects viscosity, it already captures the asymmetry
1215 between upstream and downstream (which comes from the advection term) that
1216 is an important feature of actual physical flows.
1217\end_layout
1218
1219\end_deeper
1220\begin_layout Itemize
1221the solver needs a way to save/load flow fields to provide for custom ICs
1222\end_layout
1223
1224\begin_layout Itemize
11921225practical notes on using the WLSQM solver for this:
11931226\end_layout
11941227
12631263
12641264\begin_deeper
12651265\begin_layout Itemize
1266if the obstacle does not rotate, the point clouds for the boundary problems
1267 can be preprocessed (to generate the matrices)
1266if the obstacle does not rotate (or undergoes only small vibrations), the
1267 point clouds for the boundary problems can be preprocessed (to generate
1268 the matrices)
12681269\end_layout
12691270
1271\begin_layout Itemize
1272in the case of small vibrations, use the reference configuration
1273\end_layout
1274
12701275\end_deeper
12711276\begin_layout Itemize
12721277from the solution of the boundary problems, we really only need
13001300\end_layout
13011301
13021302\begin_layout Itemize
1303the derivatives are expressed in the
1303the derivatives in the
13041304\begin_inset Quotes eld
13051305\end_inset
13061306
1307fi
1308\begin_inset Quotes erd
1309\end_inset
1310
1311 array are expressed in the
1312\begin_inset Quotes eld
1313\end_inset
1314
13071315wrong
13081316\begin_inset Quotes erd
13091317\end_inset
13451345\end_deeper
13461346\end_deeper
13471347\begin_layout Itemize
1348set up the inside-domain nodes the usual way; the sets
1348set up the inside-domain nodes the usual way
1349\end_layout
1350
1351\begin_deeper
1352\begin_layout Itemize
1353the sets
13491354\begin_inset Formula $x_{k}$
13501355\end_inset
13511356
13581358\end_layout
13591359
13601360\begin_layout Itemize
1361\begin_inset Formula $f_{k}$
1362\end_inset
1363
1364 needs only function value data, this is available directly also for boundary
1365 nodes
1366\end_layout
1367
1368\end_deeper
1369\begin_layout Itemize
13611370set up the boundary nodes so that each problem instance has its
13621371\begin_inset Formula $x_{i}$
13631372\end_inset
13791379\begin_inset Formula $x_{k}$
13801380\end_inset
13811381
1382 inside the domain (the sets
1382 inside the domain
1383\end_layout
1384
1385\begin_deeper
1386\begin_layout Itemize
1387the sets
13831388\begin_inset Formula $x_{k}$
13841389\end_inset
13851390
1386 contain no boundary nodes)
1391 contain no boundary nodes
13871392\end_layout
13881393
13891394\begin_layout Itemize
1390at each timestep:
1395use the local
1396\begin_inset Formula $(\tau,n)$
1397\end_inset
1398
1399 coordinate system to express
1400\begin_inset Formula $x_{k}$
1401\end_inset
1402
1403; we can choose
1404\begin_inset Formula $x_{i}=(0,0)$
1405\end_inset
1406
1407 for each problem instance
13911408\end_layout
13921409
1410\end_deeper
1411\begin_layout Itemize
1412at each timestep:
1413\series bold
1414###NO, don't use this algorithm — correct algorithm provided further below.
1415 But the general notes here are fine.
1416\end_layout
1417
13931418\begin_deeper
13941419\begin_layout Itemize
13951420first update advection for the inside-domain DOFs
21982198\begin_inset Formula $\boldsymbol{u}$
21992199\end_inset
22002200
2201, equation
2202\begin_inset CommandInset ref
2203LatexCommand eqref
2204reference "eq:u-pressure"
2205
2206\end_inset
2207
22012208, written at the obstacle boundary:
22022209\begin_inset Formula
22032210\[
22602260
22612261\begin_deeper
22622262\begin_layout Itemize
2263All these relations hold at the beginning of the timestep,
2264\begin_inset Formula $t=t_{0}$
2265\end_inset
2266
2267.
2268\end_layout
2269
2270\begin_layout Itemize
22632271In practice, to mitigate discretization error (especially at curved boundaries),
2264 on input it is preferable to take
2272 on input it is preferable to use a trick: take
22652273\begin_inset Formula $u_{\tau}=\mathrm{sgn}(\widehat{\boldsymbol{\tau}}\cdot\boldsymbol{u})\left|\boldsymbol{u}\right|$
22662274\end_inset
22672275
22822282\end_inset
22832283
22842284 representation).
2285 On output, use the updated
2285\end_layout
2286
2287\begin_layout Itemize
2288If the representation of
2289\begin_inset Formula $\boldsymbol{u}$
2290\end_inset
2291
2292 at the boundaries is stored in
2293\begin_inset Formula $(\tau,n)$
2294\end_inset
2295
2296 coordinates, this input conversion is equivalent to a do-nothing.
2297\end_layout
2298
2299\begin_layout Itemize
2300On output, use the updated
22862301\begin_inset Formula $u_{\tau}^{*}$
22872302\end_inset
22882303
23872387
23882388\end_deeper
23892389\begin_layout Itemize
2390Consider now the case of moving obstacles.
2391 The local time derivative of the velocity field in (
2390Just to be sure, consider also the source/sink subproblem for
2391\begin_inset Formula $\rho$
2392\end_inset
2393
2394 at the obstacle boundary:
2395\begin_inset Formula
2396\[
2397\frac{\partial\rho}{\partial t}\vert_{\Gamma_{N}}=-\rho\nabla\cdot\boldsymbol{u}=-\rho(\frac{\partial u_{\tau}}{\partial\tau}+\frac{\partial u_{n}}{\partial n})
2398\]
2399
2400\end_inset
2401
2402The (physically appropriate) BCs tell us nothing about the components of
2403
2404\begin_inset Formula $\nabla\cdot\boldsymbol{u}$
2405\end_inset
2406
2407; thus this equation produces no new constraints.
2408\end_layout
2409
2410\begin_deeper
2411\begin_layout Itemize
2412This also holds in the case of moving obstacles.
2413\end_layout
2414
2415\end_deeper
2416\begin_layout Itemize
2417Summary: for the coupled subproblems, with
2418\begin_inset Formula $\boldsymbol{g}\equiv0$
2419\end_inset
2420
2421, the obstacle boundary conditions
2422\series bold
2423for stationary obstacles
2424\series default
2425 are
2426\begin_inset Formula
2427\[
2428\frac{\partial u_{\tau}}{\partial t}\vert_{\Gamma_{N}}=-\frac{\alpha}{\rho}\widehat{\boldsymbol{\tau}}\cdot\nabla\rho\;,
2429\]
2430
2431\end_inset
2432
2433
2434\begin_inset Formula
2435\[
2436u_{n}=0\;,
2437\]
2438
2439\end_inset
2440
2441
2442\begin_inset Formula
2443\[
2444(\widehat{\boldsymbol{n}}\cdot\nabla\rho)\vert_{\Gamma_{N}}=0\;.
2445\]
2446
2447\end_inset
2448
2449
2450\end_layout
2451
2452\begin_deeper
2453\begin_layout Itemize
2454In fact, if the second and third equation are applied first, the first equation
2455 is then just the PDE itself (the coupled subproblem for
2456\begin_inset Formula $\boldsymbol{u}$
2457\end_inset
2458
2459, equation
2460\begin_inset CommandInset ref
2461LatexCommand eqref
2462reference "eq:u-pressure"
2463
2464\end_inset
2465
2466), written at the obstacle boundary.
2467\end_layout
2468
2469\begin_layout Itemize
2470The second equation is the original BC for
2471\begin_inset Formula $\boldsymbol{u}$
2472\end_inset
2473
2474.
2475\end_layout
2476
2477\begin_layout Itemize
2478The third equation is a consequence of
2479\begin_inset Formula $u_{n}=\text{\mathrm{const.}}$
2480\end_inset
2481
2482 (whence
2483\begin_inset Formula $\partial u_{n}/\partial t=0$
2484\end_inset
2485
2486) and the PDE
2487\begin_inset CommandInset ref
2488LatexCommand eqref
2489reference "eq:u-pressure"
2490
2491\end_inset
2492
2493.
2494\end_layout
2495
2496\begin_layout Itemize
2497In the solver, at the boundary,
2498\begin_inset Formula $\widehat{\boldsymbol{\tau}}\cdot\nabla\rho$
2499\end_inset
2500
2501 is conveniently available in the local
2502\begin_inset Formula $(\tau,n)$
2503\end_inset
2504
2505 representation.
2506\end_layout
2507
2508\end_deeper
2509\begin_layout Itemize
2510Conversion of vectors from
2511\begin_inset Formula $(\tau,n)$
2512\end_inset
2513
2514 components to
2515\begin_inset Formula $(x,y)$
2516\end_inset
2517
2518 components.
2519\end_layout
2520
2521\begin_deeper
2522\begin_layout Itemize
2523Once we have the new tangential velocity
2524\begin_inset Formula $u_{\tau}^{*}$
2525\end_inset
2526
2527, in order to obtain the new
2528\begin_inset Formula $u_{x}^{*}$
2529\end_inset
2530
2531 and
2532\begin_inset Formula $u_{y}^{*}$
2533\end_inset
2534
2535, we have the linear equation system
2536\begin_inset Formula
2537\begin{align*}
2538cu_{x}^{*}+su_{y}^{*} & =u_{\tau}^{*}\;,\\
2539-su_{x}^{*}+cu_{y}^{*} & =u_{n}^{*}\;,
2540\end{align*}
2541
2542\end_inset
2543
2544where
2545\begin_inset Formula $u_{n}^{*}=0$
2546\end_inset
2547
2548.
2549 In matrix form,
2550\begin_inset Formula
2551\[
2552\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\;,
2553\]
2554
2555\end_inset
2556
2557where
2558\begin_inset Formula
2559\[
2560\boldsymbol{A}=\left[\begin{array}{rr}
2561c & s\\
2562-s & c
2563\end{array}\right]\;,\quad\boldsymbol{x}=\left[\begin{array}{c}
2564u_{x}^{*}\\
2565u_{y}^{*}
2566\end{array}\right]\;,\quad\boldsymbol{b}=\left[\begin{array}{c}
2567u_{\tau}^{*}\\
2568u_{n}^{*}
2569\end{array}\right]
2570\]
2571
2572\end_inset
2573
2574Taking note that since
2575\begin_inset Formula $\boldsymbol{A}$
2576\end_inset
2577
2578 represents a rotation, it is orthogonal,
2579\begin_inset Formula $\boldsymbol{A}^{-1}=\boldsymbol{A}^{\mathrm{T}}$
2580\end_inset
2581
2582:
2583\begin_inset Formula
2584\[
2585\boldsymbol{A}^{\mathrm{T}}\boldsymbol{A}=\left[\begin{array}{rr}
2586c & -s\\
2587s & c
2588\end{array}\right]\left[\begin{array}{rr}
2589c & s\\
2590-s & c
2591\end{array}\right]=\left[\begin{array}{rr}
2592c^{2}+s^{2} & cs-sc\\
2593sc-cs & s^{2}+c^{2}
2594\end{array}\right]=\left[\begin{array}{rr}
25951 & 0\\
25960 & 1
2597\end{array}\right]\;,
2598\]
2599
2600\end_inset
2601
2602so in general,
2603\begin_inset Formula
2604\[
2605\left[\begin{array}{c}
2606u_{x}^{*}\\
2607u_{y}^{*}
2608\end{array}\right]=\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{b}=\boldsymbol{A}^{\mathrm{T}}\boldsymbol{b}=\left[\begin{array}{rr}
2609c & -s\\
2610s & c
2611\end{array}\right]\left[\begin{array}{c}
2612u_{\tau}^{*}\\
2613u_{n}^{*}
2614\end{array}\right]
2615\]
2616
2617\end_inset
2618
2619Thus, explicitly for the velocity,
2620\begin_inset Formula
2621\[
2622u_{x}^{*}=cu_{\tau}^{*}\;,\quad u_{y}^{*}=su_{\tau}^{*}\;,
2623\]
2624
2625\end_inset
2626
2627since
2628\begin_inset Formula $u_{n}^{*}=0$
2629\end_inset
2630
2631 by the boundary condition
2632\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\boldsymbol{u}\vert_{\Gamma_{N}}=0$
2633\end_inset
2634
2635.
2636\end_layout
2637
2638\end_deeper
2639\begin_layout Itemize
2640Consider now
2641\series bold
2642moving obstacles
2643\series default
2644.
2645 (Rotation, translation and deformation are allowed.) The local time derivative
2646 of the flow velocity field in arbitrary orthonormal (
23922647\begin_inset Formula $\tau,n$
23932648\end_inset
23942649
26892689\begin_inset Formula $\boldsymbol{g}\equiv0$
26902690\end_inset
26912691
2692, we obtain the equations
2692, we obtain the tangential and normal components of equation
2693\begin_inset CommandInset ref
2694LatexCommand eqref
2695reference "eq:u-pressure"
2696
2697\end_inset
2698
2699 as
26932700\begin_inset Formula
26942701\begin{align*}
26952702\widehat{\boldsymbol{\tau}}\cdot\frac{\partial\boldsymbol{u}}{\partial t}=\frac{\partial u_{\tau}}{\partial t}+u_{n}\frac{\partial\theta}{\partial t} & =-\frac{\alpha}{\rho}\widehat{\boldsymbol{\tau}}\cdot\nabla\rho\;,\\
27052705
27062706\end_inset
27072707
2708Accounting for the fact that
2708Accounting for the fact that at the boundary,
27092709\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot(\boldsymbol{u}-\boldsymbol{u}_{\mathrm{obst}})\equiv0$
27102710\end_inset
27112711
27232723
27242724\end_inset
27252725
2726On the RHS in both equations, in an explicit time-integration scheme we
2727 use only
2728\series bold
2726
2727\end_layout
2728
2729\begin_deeper
2730\begin_layout Itemize
2731This accounts for both normal-direction acceleration and rotation of the
2732 obstacle.
2733\end_layout
2734
2735\begin_deeper
2736\begin_layout Itemize
2737This reduces correctly.
2738 If the obstacle does not move (
2739\begin_inset Formula $\theta$
2740\end_inset
2741
2742 constant in time (but still local to each point on the surface!),
2743\begin_inset Formula $u_{n,\mathrm{obst}}\equiv0$
2744\end_inset
2745
2746), we again have
2747\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\nabla\rho=0$
2748\end_inset
2749
2750.
2751 Note also
2752\begin_inset Formula $u_{n}=u_{n,\mathrm{obst}}=0$
2753\end_inset
2754
2755.
2756 The first equation reduces to the previous PDE for the tangential velocity
2757 component
2758\begin_inset Formula $u_{\tau}$
2759\end_inset
2760
2761.
2762\end_layout
2763
2764\end_deeper
2765\begin_layout Itemize
2766The LHS will be obtained for the same point of time
2767\begin_inset Formula $t$
2768\end_inset
2769
2770 from which data is fed to the RHS.
2771 If we use
2772\begin_inset Quotes eld
2773\end_inset
2774
27292775old
2730\series default
2731 data.
2732 All quantities on the RHS are known.
2776\begin_inset Quotes erd
2777\end_inset
2778
2779 data at
2780\begin_inset Formula $t=t_{0}$
2781\end_inset
2782
2783:
27332784\end_layout
27342785
27352786\begin_deeper
28132813\end_layout
28142814
28152815\begin_layout Itemize
2816This accounts for both normal-direction acceleration and rotation of the
2817 obstacle, but also the result is still expressed in the
2816The result is also expressed in the
28182817\begin_inset Quotes eld
28192818\end_inset
28202819
28282828 system.
28292829\end_layout
28302830
2831\end_deeper
28312832\begin_layout Itemize
2833For FSI problems,
2834\begin_inset Formula $\partial\theta/\partial t$
2835\end_inset
2836
2837 and
2838\begin_inset Formula $u_{n,\mathrm{obst}}$
2839\end_inset
2840
2841 may be challenging to obtain for
2842\begin_inset Formula $t>t_{0}$
2843\end_inset
2844
2845 before the timestep is complete (these are needed for RK2 integration of
2846
2847\begin_inset Formula $u_{\tau}$
2848\end_inset
2849
2850).
2851 Must model somehow.
2852\end_layout
2853
2854\begin_deeper
2855\begin_layout Itemize
2856Fixed point iteration possible (IG = use the values at
2857\begin_inset Formula $t=t_{0}$
2858\end_inset
2859
2860), but very slow, since it must iterate between fluid and structure solvers.
2861\end_layout
2862
2863\begin_layout Itemize
2864Trivial model (values from the structural model remain constant across the
2865 timestep) also possible.
2866 Fast, but not very accurate.
2867 Maybe good enough?
2868\end_layout
2869
2870\begin_layout Itemize
2871Structure solver must in any case provide not only
2872\begin_inset Formula $\boldsymbol{u}_{\mathrm{obst}}$
2873\end_inset
2874
2875 (velocity, not displacement!) and
2876\begin_inset Formula $\theta$
2877\end_inset
2878
2879, but also
2880\begin_inset Formula $\partial\theta/\partial t$
2881\end_inset
2882
2883 and
2884\begin_inset Formula $\partial\boldsymbol{u}_{\mathrm{obst}}/\partial t$
2885\end_inset
2886
2887.
2888\end_layout
2889
2890\begin_layout Itemize
2891If the structural discretization uses a first-order form with displacement
2892 and velocity DOFs,
2893\begin_inset Formula $\theta$
2894\end_inset
2895
2896 and
2897\begin_inset Formula $\partial\theta/\partial t$
2898\end_inset
2899
2900 are easy to obtain.
2901 The acceleration can be computed for each structural DOF from the first-order
2902 ODE system (to get it as accurately as possible).
2903 The result can be converted to obstacle acceleration in a format that the
2904 fluid solver understands.
2905\end_layout
2906
2907\end_deeper
2908\begin_layout Itemize
28322909If we want to have
28332910\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\nabla\rho$
28342911\end_inset
29852985\begin_inset Formula $\boldsymbol{u}$
29862986\end_inset
29872987
2988 at the boundary is expressed in
2988 at the boundary is stored in
29892989\begin_inset Formula $(\tau,n)$
29902990\end_inset
29912991
30413041
30423042\begin_deeper
30433043\begin_layout Itemize
3044This is the case for the papermaking application.
3045\end_layout
3046
3047\begin_layout Itemize
30483044Each surface segment oscillates around its local reference position and
30493045 orientation.
30503046 In this case we can probably get away with not actually updating the point
30613061\begin_inset Quotes erd
30623062\end_inset
30633063
3064 at
3065\begin_inset Formula $\theta=\theta_{0}$
3066\end_inset
3067
3068.
3064 (computed for the reference configuration).
30693065\end_layout
30703066
30713067\begin_layout Itemize
30723068This gives major savings in computation time, because the boundary point
30733069 clouds can be preprocessed just like for classical rigid stationary obstacles.
3070
30743071\end_layout
30753072
30763073\begin_layout Itemize
3077The obstacle then remains geometrically stationary, but its vibration is
3078 seen by the boundary conditions, thereby affecting the flow.
3074This makes the fluid solver just as fast as for the rigid obstacle case.
3075 Useful for realtime applications in FSI of structures undergoing small
3076 vibrations?
3077\end_layout
3078
3079\begin_layout Itemize
3080In the fluid model, the obstacle then remains geometrically stationary,
3081 but its vibration is seen by the no-penetration boundary conditions, thereby
3082 affecting the flow.
30793083 This is an adaptation (for Euler flows) of the approach that was used for
30803084 potential flows in our earlier studies.
30813085\end_layout
30823086
3087\begin_layout Itemize
3088This can be used in the papermaking application having small transverse
3089 deflections.
3090 The axial motion of the paper web is not seen by the Euler flow surrounding
3091 it.
3092\end_layout
3093
3094\begin_layout Itemize
3095Even for a Navier–Stokes flow, in this application the structural geometry
3096 remains stationary in the Eulerian frame, and only the
3097\begin_inset Quotes eld
3098\end_inset
3099
3100zero relative velocity
3101\begin_inset Quotes erd
3102\end_inset
3103
3104 BC is affected (since the obstacle surface moves axially).
3105\end_layout
3106
3107\begin_layout Itemize
3108Could test this with a simpler problem, e.g.
3109 a kinematically oscillating rigid cylinder.
3110 This gives a one-way coupling, with a trivial structural model driving
3111 the fluid model.
3112 (With applications to 2D musical instruments and/or spatial audio processing?
3113 Record the numerical pressure response at some point(s) in the domain.
3114 Probably too expensive computationally, though; Helmholtz is much cheaper
3115 than an Euler flow.)
3116\end_layout
3117
30833118\end_deeper
30843119\begin_layout Itemize
30853120For compressible Navier–Stokes, doing this is even easier: then at the obstacle
31883188\end_inset
31893189
31903190 are constant in space.
3191 The result can then be rotated into
3191 The result can then be rotated into the local
31923192\begin_inset Formula $(\tau,n)$
31933193\end_inset
31943194
32103210\end_deeper
32113211\end_deeper
32123212\begin_layout Itemize
3213Just to be sure, consider also the source/sink subproblem for
3214\begin_inset Formula $\rho$
3215\end_inset
3216
3217 at the obstacle boundary:
3218\begin_inset Formula
3219\[
3220\frac{\partial\rho}{\partial t}\vert_{\Gamma_{N}}=-\rho\nabla\cdot\boldsymbol{u}=-\rho(\frac{\partial u_{\tau}}{\partial\tau}+\frac{\partial u_{n}}{\partial n})
3221\]
3222
3223\end_inset
3224
3225The (physically appropriate) BCs tell us nothing about the components of
3226
3227\begin_inset Formula $\nabla\cdot\boldsymbol{u}$
3228\end_inset
3229
3230; thus this equation produces no new constraints.
3231\end_layout
3232
3233\begin_deeper
3234\begin_layout Itemize
3235This also holds in the case of moving obstacles.
3236\end_layout
3237
3238\end_deeper
3239\begin_layout Itemize
3240Summary: for the coupled subproblems, with
3241\begin_inset Formula $\boldsymbol{g}\equiv0$
3242\end_inset
3243
3244, the obstacle boundary conditions are
3245\begin_inset Formula
3246\[
3247\frac{\partial u_{\tau}}{\partial t}\vert_{\Gamma_{N}}=-\frac{\alpha}{\rho}\widehat{\boldsymbol{\tau}}\cdot\nabla\rho\;,
3248\]
3249
3250\end_inset
3251
3252
3253\begin_inset Formula
3254\[
3255u_{n}=0\;,
3256\]
3257
3258\end_inset
3259
3260
3261\begin_inset Formula
3262\[
3263(\widehat{\boldsymbol{n}}\cdot\nabla\rho)\vert_{\Gamma_{N}}=0\;.
3264\]
3265
3266\end_inset
3267
3268
3269\end_layout
3270
3271\begin_deeper
3272\begin_layout Itemize
3273In fact, if the second and third equation are applied first, the first equation
3274 is then just the PDE itself (the coupled subproblem for
3275\begin_inset Formula $\boldsymbol{u}$
3276\end_inset
3277
3278), written at the obstacle boundary.
3279\end_layout
3280
3281\begin_layout Itemize
3282The second equation is the original BC.
3283\end_layout
3284
3285\begin_layout Itemize
3286The third equation is a consequence of
3287\begin_inset Formula $u_{n}=\text{\mathrm{const.}}$
3288\end_inset
3289
3290, whence
3291\begin_inset Formula $\partial u_{n}/\partial t=0$
3292\end_inset
3293
3294, and the PDE.
3295\end_layout
3296
3297\begin_layout Itemize
3298In the solver, at the boundary,
3299\begin_inset Formula $\widehat{\boldsymbol{\tau}}\cdot\nabla\rho$
3300\end_inset
3301
3302 is conveniently available in the local
3303\begin_inset Formula $(\tau,n)$
3304\end_inset
3305
3306 representation.
3307\end_layout
3308
3309\end_deeper
3310\begin_layout Itemize
33113213Inflow boundaries.
33123214 Here we have
33133215\begin_inset Formula $\boldsymbol{u}=\boldsymbol{u}_{\infty}$
33543354We have obtained a separate update formula for the boundary nodes.
33553355\end_layout
33563356
3357\begin_deeper
33573358\begin_layout Itemize
3359
3360\series bold
3361But we don't need it — we already have the BC update above.
3362 (See algorithm below.)
3363\end_layout
3364
3365\end_deeper
3366\begin_layout Itemize
33583367On the other hand, the first line says this is 1D advection of the field
33593368
33603369\begin_inset Formula $\rho$
35063506\end_deeper
35073507\end_deeper
35083508\begin_layout Itemize
3509advection with MacCormack is
3510\begin_inset Formula $O((\Delta t)^{2})$
3511\end_inset
3512
3513.
3514\end_layout
3515
3516\begin_layout Itemize
3517stiff part could be integrated using RK2, also
3518\begin_inset Formula $O((\Delta t)^{2})$
3519\end_inset
3520
3521.
3522\end_layout
3523
3524\begin_layout Itemize
3525boundary nodes can easily get the old value, but how about the
3526\begin_inset Formula $u_{\tau}$
3527\end_inset
3528
3529 update that gets only
3530\begin_inset Formula $\partial u_{\tau}/\partial t$
3531\end_inset
3532
3533 and needs to integrate in time? It seems we must interleave the RK2 integration
3534 of the inside-domain part with the boundary update?
3535\end_layout
3536
3537\begin_deeper
3538\begin_layout Itemize
3539Let's do this...
3540 RK2:
3541\begin_inset Formula
3542\begin{align*}
3543\boldsymbol{k}_{1} & =\boldsymbol{f}(\boldsymbol{u}_{n},\,t_{n})\;,\\
3544\boldsymbol{k}_{2} & =\boldsymbol{f}(\boldsymbol{u}_{n}+\beta\,\Delta t\,\boldsymbol{k}_{1},\,t_{n}+\beta\,\Delta t)\;,\\
3545\boldsymbol{u}_{n+1} & =\boldsymbol{u}_{n}+\Delta t\left[(1-\frac{1}{2\beta})\boldsymbol{k}_{1}+\frac{1}{2\beta}\boldsymbol{k}_{2}\right]\;,
3546\end{align*}
3547
3548\end_inset
3549
3550where
3551\begin_inset Formula $\beta\in(0,1]$
3552\end_inset
3553
3554.
3555 Classical choices
3556\begin_inset Formula $\beta=1/2$
3557\end_inset
3558
3559 (explicit midpoint method),
3560\begin_inset Formula $\beta=2/3$
3561\end_inset
3562
3563 (Ralston's method),
3564\begin_inset Formula $\beta=1$
3565\end_inset
3566
3567 (Heun's method / explicit trapezoid rule).
3568\end_layout
3569
3570\begin_layout Itemize
3571Exiting advection step, entering stiff part update.
3572\end_layout
3573
3574\begin_layout Itemize
3575At the beginning of a timestep, enforce boundary conditions for
3576\begin_inset Formula $u_{n}$
3577\end_inset
3578
3579 and
3580\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\nabla\rho$
3581\end_inset
3582
3583 at
3584\begin_inset Formula $t=t_{0}$
3585\end_inset
3586
3587
3588\end_layout
3589
3590\begin_layout Itemize
3591Compute
3592\begin_inset Formula $\boldsymbol{k}_{1}$
3593\end_inset
3594
3595 simultaneously in the interior and on the boundaries
3596\end_layout
3597
3598\begin_deeper
3599\begin_layout Itemize
3600This needs data only from
3601\begin_inset Formula $t=t_{0}$
3602\end_inset
3603
3604
3605\end_layout
3606
3607\begin_layout Itemize
3608We have also
3609\begin_inset Formula $u_{\tau}$
3610\end_inset
3611
3612 at
3613\begin_inset Formula $t=t_{0}$
3614\end_inset
3615
3616
3617\end_layout
3618
3619\begin_layout Itemize
3620\begin_inset Formula $\boldsymbol{k}_{1}$
3621\end_inset
3622
3623 contains
3624\begin_inset Formula $\partial u_{\tau}/\partial t$
3625\end_inset
3626
3627 at
3628\begin_inset Formula $t=t_{0}$
3629\end_inset
3630
3631
3632\end_layout
3633
3634\end_deeper
3635\begin_layout Itemize
3636Repeat for
3637\begin_inset Formula $\boldsymbol{k}_{2}$
3638\end_inset
3639
3640
3641\end_layout
3642
3643\begin_deeper
3644\begin_layout Itemize
3645As input, we have
3646\begin_inset Quotes eld
3647\end_inset
3648
3649candidate data
3650\begin_inset Quotes erd
3651\end_inset
3652
3653
3654\begin_inset Formula $\boldsymbol{u}_{n}+\beta\,\Delta t\,\boldsymbol{k}_{1}$
3655\end_inset
3656
3657 (with
3658\begin_inset Formula $\beta=1$
3659\end_inset
3660
3661, this is the forward Euler prediction)
3662\end_layout
3663
3664\begin_layout Itemize
3665Enforce boundary conditions for
3666\begin_inset Formula $u_{n}$
3667\end_inset
3668
3669 and
3670\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\nabla\rho$
3671\end_inset
3672
3673 at
3674\begin_inset Formula $t=t_{0}+\beta\,\Delta t$
3675\end_inset
3676
3677
3678\end_layout
3679
3680\begin_layout Itemize
3681\begin_inset Formula $\boldsymbol{k}_{2}$
3682\end_inset
3683
3684 contains an approximation of
3685\begin_inset Formula $\partial u_{\tau}/\partial t$
3686\end_inset
3687
3688 at
3689\begin_inset Formula $t_{0}+\beta\,\Delta t$
3690\end_inset
3691
3692
3693\end_layout
3694
3695\begin_layout Itemize
3696How about FSI? Must model the behavior of the structure across the timestep
3697 before we actually have it...
3698 (
3699\begin_inset Formula $\boldsymbol{u}_{\mathrm{obst}}$
3700\end_inset
3701
3702 and
3703\begin_inset Formula $\theta$
3704\end_inset
3705
3706 not yet available at
3707\begin_inset Formula $t_{0}+\beta\,\Delta t$
3708\end_inset
3709
3710)
3711\end_layout
3712
3713\end_deeper
3714\begin_layout Itemize
3715RK2 also outputs updated
3716\begin_inset Formula $u_{\tau}$
3717\end_inset
3718
3719
3720\end_layout
3721
3722\begin_layout Itemize
3723\begin_inset Formula $u_{\tau}$
3724\end_inset
3725
3726 behaves otherwise like an inside-domain DOF, but with a different update
3727 formula.
3728\end_layout
3729
3730\begin_layout Itemize
3731After computing the last
3732\begin_inset Formula $\boldsymbol{u}_{n+1}$
3733\end_inset
3734
3735 for this stiff part update, enforce boundary conditions for
3736\begin_inset Formula $u_{n}$
3737\end_inset
3738
3739 and
3740\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\nabla\rho$
3741\end_inset
3742
3743 at
3744\begin_inset Formula $t_{0}+\Delta t$
3745\end_inset
3746
3747
3748\end_layout
3749
3750\begin_deeper
3751\begin_layout Itemize
3752May be important, because the advection step comes next
3753\end_layout
3754
3755\end_deeper
3756\end_deeper
3757\begin_layout Itemize
3758
3759\series bold
3760The stuff below here is not needed.
3761\end_layout
3762
3763\begin_layout Itemize
35093764From the advection subproblem for
35103765\begin_inset Formula $\boldsymbol{u}$
35113766\end_inset
40544054
40554055\end_deeper
40564056\end_deeper
4057\begin_layout Itemize
4058Once we have the new tangential velocity
4059\begin_inset Formula $u_{\tau}^{*}$
4060\end_inset
4061
4062, in order to obtain the new
4063\begin_inset Formula $u_{x}^{*}$
4064\end_inset
4065
4066 and
4067\begin_inset Formula $u_{y}^{*}$
4068\end_inset
4069
4070, we have the linear equation system
4071\begin_inset Formula
4072\begin{align*}
4073cu_{x}^{*}+su_{y}^{*} & =u_{\tau}^{*}\;,\\
4074-su_{x}^{*}+cu_{y}^{*} & =u_{n}^{*}\;,
4075\end{align*}
4076
4077\end_inset
4078
4079where
4080\begin_inset Formula $u_{n}^{*}=0$
4081\end_inset
4082
4083.
4084 In matrix form,
4085\begin_inset Formula
4086\[
4087\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\;,
4088\]
4089
4090\end_inset
4091
4092where
4093\begin_inset Formula
4094\[
4095\boldsymbol{A}=\left[\begin{array}{rr}
4096c & s\\
4097-s & c
4098\end{array}\right]\;,\quad\boldsymbol{x}=\left[\begin{array}{c}
4099u_{x}^{*}\\
4100u_{y}^{*}
4101\end{array}\right]\;,\quad\boldsymbol{b}=\left[\begin{array}{c}
4102u_{\tau}^{*}\\
4103u_{n}^{*}
4104\end{array}\right]
4105\]
4106
4107\end_inset
4108
4109Taking note that since
4110\begin_inset Formula $\boldsymbol{A}$
4111\end_inset
4112
4113 represents a rotation, it is orthogonal,
4114\begin_inset Formula $\boldsymbol{A}^{-1}=\boldsymbol{A}^{\mathrm{T}}$
4115\end_inset
4116
4117:
4118\begin_inset Formula
4119\[
4120\boldsymbol{A}^{\mathrm{T}}\boldsymbol{A}=\left[\begin{array}{rr}
4121c & -s\\
4122s & c
4123\end{array}\right]\left[\begin{array}{rr}
4124c & s\\
4125-s & c
4126\end{array}\right]=\left[\begin{array}{rr}
4127c^{2}+s^{2} & cs-sc\\
4128sc-cs & s^{2}+c^{2}
4129\end{array}\right]=\left[\begin{array}{rr}
41301 & 0\\
41310 & 1
4132\end{array}\right]\;,
4133\]
4134
4135\end_inset
4136
4137so
4138\begin_inset Formula
4139\[
4140\left[\begin{array}{c}
4141u_{x}^{*}\\
4142u_{y}^{*}
4143\end{array}\right]=\boldsymbol{x}=\boldsymbol{A}^{-1}\boldsymbol{b}=\boldsymbol{A}^{\mathrm{T}}\boldsymbol{b}=\left[\begin{array}{rr}
4144c & -s\\
4145s & c
4146\end{array}\right]\left[\begin{array}{c}
4147u_{\tau}^{*}\\
4148u_{n}^{*}
4149\end{array}\right]
4150\]
4151
4152\end_inset
4153
4154Thus, explicitly,
4155\begin_inset Formula
4156\[
4157u_{x}^{*}=cu_{\tau}^{*}\;,\quad u_{y}^{*}=su_{\tau}^{*}\;,
4158\]
4159
4160\end_inset
4161
4162since
4163\begin_inset Formula $u_{n}^{*}=0$
4164\end_inset
4165
4166 by the boundary condition
4167\begin_inset Formula $\widehat{\boldsymbol{n}}\cdot\boldsymbol{u}\vert_{\Gamma_{N}}=0$
4168\end_inset
4169
4170.
4171\end_layout
4172
41734057\end_deeper
41744058\begin_layout Standard
41754059The advection problems can be handled by the semi-Lagrangean approach with