$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$
There is a deep connection between holomorphic functions and harmonic functions in dimension 2. Namely, if $f(z)$ is a holomorphic function then $u(x,y)=\Re f(x+iy)$ and $v(x,y)=\Im f(x+iy)$ are harmonic functions.
Conversely, if $\Omega$ is a simply connected domain, then each harmonic function is a real part of some holomorphic function $f(z)$ (and we do not distinguish $(x,y)$ and $z=x+iy$) (if $\Omega$ is connected but not simply connected, $f(z)$ could be multivalued).
Such relations are due to the fact, that \begin{align} &\Delta =4\partial_z \bar{\partial}_z, \label{eq-14.4.1}\\ &\partial_z=\frac{1}{2}(\partial_x-i\partial_y), \label{eq-14.4.2}\\ &\bar{\partial}_z=\frac{1}{2}(\partial_x+i\partial_y) &&=\partial_{\bar{z}}. \label{eq-14.4.3} \end{align}
The one-to-one maps $z\mapsto w(z)$ with holomorphic $w(z)$, s.th. $w'(z)$ does not vanish, are called conformal maps. These are maps, which preserve angles and orientation. Then \begin{equation} \Delta_z= |w'(z)|^2 \Delta_w \label{eq-14.4.4} \end{equation} Because of this solving Dirichlet problem \begin{align} &\Delta u =0 &&\text{in } \Omega, \label{eq-14.4.5}\\ &u|_\Gamma = \phi, \label{eq-14.4.6} \end{align} where $\Gamma=\partial\Omega$ is the boundary of $\Omega$ could be reduced to finding a conformal map of $\Omega$ into some standard domain $\mathcal{D}$, for which we can solve such problem (f.e. half-plane or a unit disk).
Further, for such domains Neumann problem for (\ref{eq-14.4.5}): \begin{align} &\frac{\partial u}{\partial \nu}|_\Gamma = \psi, \label{eq-14.4.7} \end{align} where $\nu$ is an inner unit normal to $\Gamma$, could be reduced to Dirichlet problem. Indeed, let us parametrize $\Gamma$ by $s$, the length of the counter-clockwise path along $\Gamma$ from some point $z_0$ to the current point $z$. Then \begin{equation} \frac{\partial u}{\partial \nu}=-\frac{\partial v}{\partial s}, \label{eq-14.4.8} \end{equation} where $v=\Im f(z)$ while $u=\Re f(z)$. Therefore, we first find $\phi(s)= \Im f(z)$ from $\phi'(s)=-\psi(s)$. $\phi(s)$ is defined up to a constant by integration. However, to have it properly defined as a single-valued function, we need \begin{equation} \int_\Gamma \psi(s)\,ds=0, \label{eq-14.4.9} \end{equation} which is a necessary and sufficient condition for $u$, satisfying Neumann problem to exist. Then we find $v(x,y)$ as a solution to $\Delta v=0$ in $\Omega$ $v|_\gamma =\phi$, and then $u(x,y)$. Recall that \begin{equation} \partial_x u=\partial_y v,\qquad \partial_y u=-\partial_x v. \label{eq-14.4.10} \end{equation}
Remark 1. Classical aerodynamics (2D subsonic aerodynamics) is a basically theory of 2D Laplace equation and thus a complex variables theory.
Biharmonic equation \begin{equation} \Delta^2 u=0 \label{eq-14.4.11} \end{equation} is not invariant with respect to conformal maps, but theory of complex variables is useful too: \begin{equation} u =\Re \bigl(\bar{z}f(z)+g(z)\bigr), \label{eq-14.4.12} \end{equation} satisfies (\ref{eq-14.4.12}), provided $f,g$ are holomorphic functions. Conversely, if $\Omega$ is simply connected then each solution to (\ref{eq-14.4.11}) is in the form (\ref{eq-14.4.12}).