# PGD for Solving the Biharmonic Equation

Guang Tao Xu, Francisco Chinesta, Adrien Leygue, Michel Visonneau

► To cite this version:

Guang Tao Xu, Francisco Chinesta, Adrien Leygue, Michel Visonneau. PGD for Solving the Biharmonic Equation. ASME 2012 11th Biennial Conference on Engineering Systems Design and Analysis, Jul 2012, Nantes, France. pp.219-223, 10.1115/ESDA2012-82484. http://onlinelibrary.wiley.com/doi/abs/10.1115/ESDA2012-82484

HAL Id: hal-01202443

https://hal.science/hal-01202443

Submitted on 18 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PGD FOR SOLVING THE BIHARMONIC EQUATION

GuangTao XU

Institut de Recherche en Genie Civil et Mecanique
Laboratoire de Mcanique des Fluides
Ecole Centrale de NANTES
Nantes 44321
FRANCE

Email: guangtao.xu@ec-nantes.fr

Francisco CHINESTA

Institut de Recherche en Genie Civil et Mecanique
Ecole Centrale de NANTES
Nantes 44321
FRANCE

Email: francisco.chinesta@ec-nantes.fr

Adrien LEYGUE

Institut de Recherche en Genie Civil et Mecanique
Ecole Centrale de NANTES
Nantes 44321
FRANCE

Email: adrien.leygue@ec-nantes.fr

Michel VISONNEAU

Laboratoire de Mecanique des Fluides
Ecole Centrale de NANTES
Nantes 44321
FRANCE

Email: michel.visonneau@ec-nantes.fr

ABSTRACT

Biharmonic problem has been raised in many research fields, such as elasticity problem in plate geometries or the Stokes flow problem formulated by using the stream function. The fourth order partial differential equation can be solved by applying many techniques. When using finite elements C¹ continuity must be assured. For this purpose Hermite interpolations constitute an appealing choice, but it imply the consideration of many degrees of freedom at each node with the consequent impact on the resulting discrete linear problem. Spectral approaches allow exponential convergence whilst a single degree of freedom is needed. However, the enforcement of boundary conditions re-mains a tricky task. In this paper we propose a separated representation of the stream function which transform the 2D solution in a sequence of 1D problems, each one be solved by using a spectral approximation.

INTRODUCTION

The Biharmonic problem [1] has been raised in many research fields [2], such as in elasticity problem which dealing with the transverse displacements of elastic plates [3] and in 2D flows when using the stream function [4].

The Finite Element Method - FEM - [5] and the Boundary Element Method - BEM - [6] can be used to solve biharmonic equations. Other methods and variants exist, e.g. [1].

The spectral method has been widely used in the solution of Partial Differential Equations - PDE -, in particular high order PDEs, e.g. [7].

The Chebyshev spectral collocation method [8, 9] has been traditionally used to solve biharmonic problems. Its main advantage lies in the fact that it only needs a degree of freedom per node and it exhibits exponential convergence rates.

In this paper spectral collocation schemes are combined with the PGD technique [10, 11] that allows a separated representation of the fields involved in the model, and then in our case, transform the solution of a 2D model into the solution of few 1D problems.
THE BIHARMONIC EQUATION

We consider the biharmonic equation:

$$ \Delta^2 u = f \quad \text{in } \Omega_x \times \Omega_y = ]-1, 1[ \times ]-1, 1[ \quad (1) $$

where

$$ \Delta^2 = \frac{\partial^4}{\partial x^4} + 2 \frac{\partial^4}{\partial x^2 \partial y^2} + \frac{\partial^4}{\partial y^4} \quad (2) $$

subjected the following boundary conditions:

$$ u = f_1 \quad \text{on } \Gamma \quad (3) $$

and

$$ \frac{\partial u}{\partial n} = f_2 \quad \text{on } \Gamma \quad (4) $$

where $u$ is the solution of the biharmonic equation. In plate theory $u$ represents the transverse displacement and in flow simulations it represents the stream function, from which the velocity can be calculated from:

$$ \left\{ \begin{array}{l} v_x = \frac{\partial u}{\partial y} \\ v_y = -\frac{\partial u}{\partial x} \end{array} \right. \quad (5) $$

PGD FOR BIHARMONIC EQUATION

In this section we illustrate the biharmonic problem by using the separated representation within the PGD framework.

The aim of the method is to compute $N$ couples of functions $(X_i(x), Y_i(y))$, $i = 1, \dots, N$ such that $X_i(x)$, $i = 1, \dots, N$ and $Y(y)$, $i = 1, \dots, N$ are defined in 1D domains. The 2D solution reads:

$$ u(x,y) = \sum_{i=1}^{N} X_i(x) \cdot Y_i(y) \quad (6) $$

The weak form of problem (1) writes: Find $u(x,y)$ verifying the boundary conditions (3) and (4) such that

$$ \int_{\Omega_x} \int_{\Omega_y} u^*(x,y) (\Delta^2 u(x,y) - f(x,y)) dx \cdot dy = 0 \quad (7) $$

for all the functions $u^*(x,y)$ in an appropriate functional space.

We now compute the functions involved in the separated representation. We suppose that the set of functional couples $(X_i(x), Y_i(y))$, $i = 1, \dots, n$ with $1 \le n < N$ are already known (they have been previously computed) and at the present iteration we search the enrichment couple $(R(x), S(y))$ by applying an alternating directions fixed-point algorithm which after convergence will constitute the next functional couple $(X_{n+1}, Y_{n+1})$. Hence at the present iteration, $n+1$, we assume the separated representation

$$ u(x,y) \approx \sum_{i=1}^{n} X_i(x) \cdot Y_i(y) + R(x) \cdot S(y) \quad (8) $$

The weighting function $u^*(x,y)$ is then assumed as

$$ u^*(x,y) = R^*(x) \cdot S(y) + R(x) \cdot S^*(y) \quad (9) $$

Introducing the trial and test function into the weak form it results

$$ \begin{aligned} & \int_{\Omega} [R^*(x) \cdot S(y) + R(x) \cdot S^*(y)] \\ & \qquad \left[ \frac{\partial^4 R(x)}{\partial x^4} \cdot S(y) + 2 \frac{\partial^2 R(x)}{\partial x^2} \cdot \frac{\partial^2 S(y)}{\partial y^2} + R(x) \cdot \frac{\partial^4 S(y)}{\partial y^4} \right] dx \cdot dy \\ & = \int_{\Omega} [R^*(x) \cdot S(y) + R(x) \cdot S^*(y)] \cdot [f(x,y) - \\ & \qquad \sum_{i=1}^{n} \left( \frac{\partial^4 X_i(x)}{\partial x^4} \cdot Y_i(y) + 2 \frac{\partial^2 X_i(x)}{\partial x^2} \cdot \frac{\partial^2 Y_i(y)}{\partial y^2} + X_i(x) \frac{\partial^4 Y_i(y)}{\partial y^4} \right)] dx \cdot dy \end{aligned} \quad (10) $$

First, we suppose that $R(x)$ is known, implying that $R^*(x) = 0$. Thus, equation (10) reads

$$ \int_{\Omega_y} S^*(y) [\alpha_{Rx} S(y) + 2\beta_{Rx} \frac{\partial^2 S(y)}{\partial y^2} + \gamma_{Rx} \frac{\partial^4 S(y)}{\partial y^4}] dy = \int_{\Omega_y} S^*(y) [\eta_{Rx}(y) - \sum_{i=1}^{n} (\alpha_{Rx}^i Y_i(y) + 2\beta_{Rx}^i \frac{\partial^2 Y_i(y)}{\partial y^2} + \gamma_{Rx}^i \frac{\partial^4 Y_i(y)}{\partial y^4})] dy \quad (11) $$

where

$$ \begin{aligned} \alpha_{Rx} &= \int_{\Omega_x} R(x) \frac{\partial R^4(x)}{\partial x^4} dx \\ \beta_{Rx} &= \int_{\Omega_x} R(x) \frac{\partial R^2(x)}{\partial x^2} dx \\ \gamma_{Rx} &= \int_{\Omega_x} R(x) R(x) dx \\ \alpha_{Rx}^i &= \int_{\Omega_x} R(x) \frac{\partial X_i^4(x)}{\partial x^4} dx \\ \beta_{Rx}^i &= \int_{\Omega_x} R(x) \frac{\partial X_i^2(x)}{\partial x^2} dx \\ \gamma_{Rx}^i &= \int_{\Omega_x} R(x) X_i(x) dx \\ \eta_{Rx}(y) &= \int_{\Omega_x} R(x) f(x,y) dx \end{aligned} \quad (12) $$

As the weak formulation is satisfied for all $S^*(y)$, we can come back to its associated strong form:

$$ \begin{aligned} & \alpha_{Rx} S(y) + 2\beta_{Rx} \frac{\partial^2 S(y)}{\partial y^2} + \gamma_{Rx} \frac{\partial^4 S(y)}{\partial y^4} = \\ & \eta_{Rx}(y) - \sum_{i=1}^{n} (\alpha_{Rx}^i Y_i(y) + 2\beta_{Rx}^i \frac{\partial^2 Y_i(y)}{\partial y^2} + \gamma_{Rx}^i \frac{\partial^4 Y_i(y)}{\partial y^4}) \end{aligned} \quad (13) $$
This fourth order equation will be solved by using a pseudo-spectral Chebyshev method.

Now, from the function $S(y)$ just computed, we search $R(x)$, in this case, $S(y)$ being known, $S^*(y)$ vanishes and Eq. (10) reads:

$$ \int_{\Omega_x} R^*(x)[\alpha_{Sy} \frac{\partial^4 R(x)}{\partial x^4} + 2\beta_{Sy} \frac{\partial^2 R(x)}{\partial x^2} + \gamma_{Sy} R(x)]dx = \\ \int_{\Omega_x} R^*(x)\eta_{Sy}(x)dx - \\ \int_{\Omega_x} R^*(x)[\sum_{i=1}^n (\alpha_{Sy}^i \frac{\partial^4 X_i(x)}{\partial x^4} + 2\beta_{Sy}^i \frac{\partial^2 X_i(x)}{\partial x^2} + \gamma_{Sy}^i X_i(x))]dx \quad (14) $$

where

$$ \alpha_{Sy} = \int_{\Omega_x} S(y)S(y)dy \\ \beta_{Sy} = \int_{\Omega_y} S(y)\frac{\partial S^2(y)}{\partial y^2}dy \\ \gamma_{Sy} = \int_{\Omega_y} S(y)\frac{\partial S^4(y)}{\partial y^4}dy \\ \alpha_{Sy}^i = \int_{\Omega_x} S(y)Y_i(y)dy \\ \beta_{Sy}^i = \int_{\Omega_y} S(y)\frac{\partial Y_i^2(y)}{\partial y^2}dy \\ \gamma_{Sy}^i = \int_{\Omega_y} S(y)\frac{\partial Y_i^4(y)}{\partial y^4}dy \\ \eta_{Sy}(x) = \int_{\Omega_y} S(y)f(x,y)dy \quad (15) $$

whose strong form reads

$$ \alpha_{Sy} \frac{\partial^4 R(x)}{\partial x^4} + 2\beta_{Sy} \frac{\partial^2 R(x)}{\partial x^2} + \gamma_{Sy} R(x) = \\ \eta_{Sy}(x) - \sum_{i=1}^n (\alpha_{Sy}^i \frac{\partial^4 X_i(x)}{\partial x^4} + 2\beta_{Sy}^i \frac{\partial^2 X_i(x)}{\partial x^2} + \gamma_{Sy}^i X_i(x)) \quad (16) $$

that will be solved again by using a pseudo-spectral Chebyshev method.

These two steps continue repeat until reaching the fixed point. If we denote the functions $R(x)$ at the present and previous iteration as $R^p(x)$ and $R^{p-1}(x)$, respectively, and the same for the function $S(y)$, $S^p(y)$ and $S^{p-1}(y)$, the error at present iteration can be defined from:

$$ e = \int_{\Omega_x \times \Omega_y} (R^p(x) \cdot S^p(y) - R^{p-1}(x) \cdot S^{p-1}(y))^2 dx \cdot dy \le \varepsilon \quad (17) $$

where $\varepsilon$ is a small enough parameter.

After the convergence we can define the next functional couple: $X_{n+1} = R$ and $Y_{n+1} = S$.

The enrichment procedure must continue until reaching the convergence, that can be evaluated from the error $E$:

$$ E = \frac{\|\Delta^2 u - f(x, y)\|}{\|f(x, y)\|} \le \tilde{\varepsilon} \quad (18) $$

with $\tilde{\varepsilon}$ another small enough parameter.

## PSEUDO-SPECTRAL COLLOCATION DISCRETIZATION

We assume the general form of a 1D fourth order differential equation:

$$ a \frac{d^4 u}{dx^4} + b \frac{d^2 u}{dx^2} + cu = g(x) \quad (19) $$

The unknown function $u(x)$ is approximated in $\Omega_x = ]-1, 1[$ from:

$$ u(x) = \sum_{i=1}^{i=M} \alpha_i \cdot T_i(x) \quad (20) $$

where $M$ denotes the number of nodes considered on $\Omega_x$, whose coordinates are given by

$$ x_i = \cos \left( \frac{(i-1) \cdot \pi}{M-1} \right), \quad i = 1, \dots, M \quad (21) $$

The interpolants $T_i(x)$ verify the Kroenecker delta property, i.e. $T_i(x_k) = \delta_{ik}$.

At each node $k$, $3 \le k \le M-2$ (the remaining 4 nodes will be used for enforcing the boundary conditions) the discrete equations writes:

$$ a \cdot \sum_{i=1}^{i=M} \alpha_i \cdot \frac{dT_i^4}{dx^4}|_{x_k} + b \cdot \sum_{i=1}^{i=M} \alpha_i \cdot \frac{dT_i^2}{dx^2}|_{x_k} + c \cdot \alpha_k = f(x_k) \quad (22) $$

When we assume that the first modes of the separated representation verified the boundary conditions (3) and (4), functions $R(x)$ and $S(y)$ are subjected to homogeneous Dirichlet and Neumann conditions. Thus, we should enforce $u(x_1) = u(x_M) = 0$ and $\frac{du}{dx}|_{x_1} = \frac{du}{dx}|_{x_M} = 0$. This conditions results in:

$$ \left\{ \begin{array}{l} \alpha_1 = 0 \\ \sum_{i=1}^{i=M} \alpha_i \cdot \frac{dT_i(x)}{dx}|_{x_1} = 0 \\ \alpha_M = 0 \\ \sum_{i=1}^{i=M} \alpha_i \cdot \frac{dT_i(x)}{dx}|_{x_M} = 0 \end{array} \right. \quad (23) $$

## NUMERICAL EXAMPLE

Let us consider the plate problem

$$ \Delta^2 u(x, y) = f(x, y) \quad \text{in } \Omega_x \times \Omega_y = ]-1, 1[ \times ]-1, 1[ \quad (24) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. EXACT SOLUTION</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. PGD BASED SPECTRAL SOLUTION</figcaption></figure>

with

$$f(x,y) = 4\cos(\pi x)\cos(\pi y) + \cos(\pi x) + \cos(\pi y) \quad (25)$$

The boundary conditions write

$$u = 0 \quad \text{on } \Gamma \quad (26)$$

and

$$\frac{\partial u}{\partial n} = 0 \quad \text{on } \Gamma \quad (27)$$

The exact solution is given by

$$u = \frac{1}{\pi^4}(1 + \cos(\pi x))(1 + \cos(\pi y)) \quad (28)$$

which is shown in Figure 1 and that serves as reference.

The solution computed by using the separated representation within the PGD framework with $M = 100$ nodes in each direction is shown in Fig. 2. Figure 3 depicts the main modes involved in the separated representation. The error with respect to the reference solution (exact solution) is depicted in Fig. 4, where the error was computed at each node. As the exact solution can be expressed from 3 functional couples, the error when considering more modes is in the order of $10^{-12}$ as noticed in Figure 5.

## CONCLUSION

In this work we analyzed the possibility of using separated representations for solving high order partial differential equa-

<figure><img src="image_4.png" /><figcaption>FIGURE 3. PGD MODES</figcaption></figure>

tions, as is the case of the biharmonic equation. The first results seem indicate that PGD and spectral techniques can be efficiently combined. The tricky point concerns the enforcement of the boundary conditions, that is, how to cote the first modes of the separated representation in order to account for the two boundary conditions known in the whole domain boundary.

## REFERENCES

[1] Marin, L., and Lesnic, D., 2005. “The method of fundamental solutions for inverse boundary value problems associated with the two-dimensional biharmonic equation”. *Mathematical and Computer Modelling*, **42**(3-4), pp. 261–278.

[2] Erturk, E., and Dursun, B., 2007. “Numerical solutions of 2-d steady incompressible flow in a driven skewed cavity”. *Journal of Applied Mathematics and Mechanics*, **87**, pp. 377–392.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. ERROR OF THE PGD SOLUTION</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. ERROR VS THE NUMBER OF MODES</figcaption></figure>

[3] Li, Z.-C., Lee, M.-G., Chiang, J. Y., and Liu, Y., 2011. “The trefftz method using fundamental solutions for biharmonic equations”. *Journal of Computational and Applied Mathematics*, **235**(15), pp. 4350–4367.

[4] Montlaur, A., Fernandez-Mendez, S., and Huerta, A., 2008. “Discontinuous galerkin methods for the stokes equations using divergence-free approximations”. *International Journal for Numerical Methods in fluids*, **57**, April, pp. 1071–1092.

[5] Gudi, T., Nataraj, N., and K.pani, A., 2007. “Mixed discontinuous galerkin finite element method for the biharmonic equation”. *Journal of Scientific Computing*, **37**(2), April, pp. 139–161.

[6] Mai-duy, N., Trancong, T., and Tanner, R., 2006. “A domain-type boundary-integral-equation method for two-dimensional biharmonic dirichlet problem”. *Engineering Analysis with Boundary Elements*, **30**(10), October, pp. 809–817.

[7] Mai-duy, N., and Tanner, R., 2007. “A spectral collo-

cation method based on integrated chebyshev polynomials for two-dimensional biharmonic boundary-value problems”. *Journal of Computational and Applied Mathematics*, **201**, pp. 30–47.

[8] Martinez, J. J., and Esperanca, P. D. T. T., 2007. “A chebyshev collocation spectral method for numerical simulation of incompressible flow problems”. *Journal of the Brazilian Society of Mechanical Sciences and Engineering*, **29**(3), pp. 317–328.

[9] Li, J., Sun, Y., and Yu, Y., 2008. “Iterative and direct chebyshev collocation spectral methods for one-dimensional radiative heat transfer”. *International Journal of Heat and Mass Transfer*, **51**(25-26), December, pp. 5887–5894.

[10] Chinesta, F., Ammar, A., Leygue, A., and Keunings, R., 2011. “An overview of the proper generalized decomposition with applications in computational rheology”. *Journal of Non-Newtonian Fluid Mechanics*, **166**, January, pp. 578–592.

[11] Chinesta, F., Ammar, A., and Cueto, E., 2011. “Proper generalized decomposition of multiscale models”. *International Journal of Numerical Methods in Engineering*, **83**(8-9), pp. 1114–1132.