|
13 | 13 | "In the case of high values of the Lamé parameter $ \\lambda $, the pressure field $ p $ can be defined as:\n", |
14 | 14 | "\n", |
15 | 15 | "$$\n", |
16 | | - "\n", |
17 | 16 | "p := -\\lambda \\nabla \\cdot u \\hspace{1cm} \\text{ in } \\Omega\n", |
18 | | - "\n", |
19 | 17 | "$$\n", |
20 | 18 | "\n", |
21 | 19 | "As $u \\in H^1(\\Omega)$, the pressure field $p \\in L^2(\\Omega)$.\n", |
22 | 20 | "\n", |
23 | 21 | "## The strong formulation of the problem\n", |
24 | 22 | "The strong formulation of the problem can be expressed as follows:\n", |
25 | | - "Find $ (u, p) \\in V \\times L^2(\\Omega) $ such that\n", |
| 23 | + "Find $(u, p) \\in V \\times L^2(\\Omega)$ such that\n", |
26 | 24 | "\n", |
27 | 25 | "$$\n", |
28 | | - "\n", |
29 | 26 | "\\begin{align}\n", |
30 | 27 | " -\\nabla \\cdot \\sigma(u,p) &= f & \\text{in } & \\Omega \\\\\n", |
31 | 28 | " u &= 0 & \\text{on } & \\partial \\Omega_D \\\\\n", |
32 | 29 | " \\sigma(u,p) \\cdot n &= g_T & \\text{on } & \\partial \\Omega_T \\\\\n", |
33 | 30 | " \\nabla \\cdot u + \\frac{1}{\\lambda} p &= 0 & \\text{in } & \\Omega\n", |
34 | 31 | "\\end{align}\n", |
35 | | - "\n", |
36 | 32 | "$$\n", |
37 | 33 | "\n", |
38 | 34 | "## The Variational Formulation (Pressure-Displacement Formulation)\n", |
39 | | - "To derive the weak formulation of the problem, we multiply the first equation by a test function $ v \\in V $ and integrate by parts over the domain $ \\Omega $ and we add the fourth equation multiplied by a test function $ q \\in L^2(\\Omega) $ integrated over the domain $ \\Omega $:\n", |
| 35 | + "To derive the weak formulation of the problem, we multiply the first equation by a test function $v \\in V$ and integrate by parts over the domain $\\Omega$ and we add the fourth equation multiplied by a test function $q \\in L^2(\\Omega)$ integrated over the domain $\\Omega$:\n", |
40 | 36 | "\n", |
41 | 37 | "$$\n", |
42 | | - "\n", |
43 | 38 | "\\boxed{\n", |
44 | 39 | "\\begin{aligned}\n", |
45 | 40 | "&\\text{Find } (u, p) \\in V \\times L^2(\\Omega) \\text{ such that:} \\\\\n", |
46 | 41 | "&\\qquad \\tilde{a}((u, p), (v, q)) = l(v, q) \\quad \\forall (v, q) \\in V \\times L^2(\\Omega)\n", |
47 | 42 | "\\end{aligned}\n", |
48 | 43 | "}\n", |
49 | | - "\n", |
50 | 44 | "$$\n", |
51 | 45 | "\n", |
52 | 46 | "with\n", |
53 | 47 | "\n", |
54 | 48 | "$$\n", |
55 | | - "\n", |
56 | 49 | "\\begin{aligned}\n", |
57 | 50 | "&\\tilde{a} : \n", |
58 | 51 | "\\begin{cases}\n", |
|
65 | 58 | "v \\longmapsto \\int_\\Omega f \\cdot v \\, dx + \\int_{\\partial\\Omega_N} t_N \\cdot v \\, ds\n", |
66 | 59 | "\\end{cases}\n", |
67 | 60 | "\\end{aligned}\n", |
68 | | - "\n", |
69 | 61 | "$$" |
70 | 62 | ] |
71 | 63 | }, |
|
469 | 461 | "metadata": {}, |
470 | 462 | "outputs": [], |
471 | 463 | "source": [ |
472 | | - "p_simulated = -lambda_ * compute_divergence(u1, u2, u3, x_vals, y_vals, z_vals)\n", |
473 | | - "exact_p = -lambda_ * np.pi * np.sin(np.pi * X) * np.sin(np.pi * Y) * np.cos(np.pi * Z)\n", |
474 | | - "error_p = p_simulated - exact_p\n", |
475 | | - "\n", |
476 | 464 | "fig = plt.figure(figsize=(24, 6 * len(z_plane_list)))\n", |
477 | 465 | "for i, z_plane in enumerate(z_plane_list):\n", |
478 | 466 | " list_index = int(z_plane * (h - 1))\n", |
479 | 467 | "\n", |
480 | 468 | " # Simulated pressure\n", |
481 | 469 | " ax1 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 1, projection='3d')\n", |
482 | | - " ax1.plot_surface(X[:, :, list_index], Y[:, :, list_index], p_simulated[:, :, list_index], cmap='viridis')\n", |
| 470 | + " ax1.plot_surface(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index], cmap='viridis')\n", |
483 | 471 | " ax1.set_title('Simulated Pressure on z={:.2f} plane'.format(z_plane))\n", |
484 | 472 | " ax1.set_xlabel('X-axis')\n", |
485 | 473 | " ax1.set_ylabel('Y-axis')\n", |
486 | | - " ax1.set_zlabel('p_simulated')\n", |
| 474 | + " ax1.set_zlabel('p_sim')\n", |
487 | 475 | " ax1.set_xlim(0, x_lim)\n", |
488 | 476 | " ax1.set_ylim(0, y_lim)\n", |
489 | | - " ax1.set_zlim(p_simulated[:, :, list_index].min(), p_simulated[:, :, list_index].max())\n", |
| 477 | + " ax1.set_zlim(p_sim[:, :, list_index].min(), p_sim[:, :, list_index].max())\n", |
490 | 478 | "\n", |
491 | 479 | " # Exact pressure\n", |
492 | 480 | " ax2 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 2, projection='3d')\n", |
|
501 | 489 | "\n", |
502 | 490 | " # Error in pressure\n", |
503 | 491 | " ax3 = fig.add_subplot(len(z_plane_list), 3, 3 * i + 3, projection='3d')\n", |
504 | | - " mappable = ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], error_p[:, :, list_index], cmap='viridis')\n", |
| 492 | + " mappable = ax3.plot_surface(X[:, :, list_index], Y[:, :, list_index], p_sim[:, :, list_index] - exact_p[:, :, list_index], cmap='viridis')\n", |
505 | 493 | " ax3.set_title('Error in Pressure on z={:.2f} plane'.format(z_plane))\n", |
506 | 494 | " ax3.set_xlabel('X-axis')\n", |
507 | 495 | " ax3.set_ylabel('Y-axis')\n", |
508 | 496 | " ax3.set_zlabel('Error in pressure')\n", |
509 | 497 | " ax3.set_xlim(0, x_lim)\n", |
510 | 498 | " ax3.set_ylim(0, y_lim)\n", |
511 | | - " ax3.set_zlim(error_p[:, :, list_index].min(), error_p[:, :, list_index].max())\n", |
| 499 | + " ax3.set_zlim((p_sim[:, :, list_index] - exact_p[:, :, list_index]).min(), (p_sim[:, :, list_index] - exact_p[:, :, list_index]).max())\n", |
512 | 500 | " fig.colorbar(mappable, ax=ax3, shrink=0.5, aspect=5)\n", |
513 | 501 | "\n", |
514 | 502 | "plt.tight_layout()\n", |
|
0 commit comments