It is possible to construct Green's functions to nondivergence form equations like this, but with caveats:
- the Green's function will have very low regularity (e.g. not locally bounded away from $x = y$)
- even the existence of a Green's function is nontrivial, and some effort is required to understand in what precise sense it "solves" the PDE written
- Green's functions are not the most typical approach to treating this problem outside of a priori estimates; the lack of symmetry and poor behavior in the $y$ variable in particular make them less relevant
The Green's function
Let us think through what would be required to construct a Green's function for this PDE. We start with some $\{L_\epsilon\}$ whose coefficients are smooth and converge (say almost everywhere) to the coefficients of $L$ as $\epsilon \rightarrow 0$. Solving $L_\epsilon u_\epsilon = f$ for a smooth $f$ with zero boundary conditions is not a problem; we can move the coefficients inside and write this as a divergence-form equation, for example. Let us call $R_\epsilon(f)$ the linear mapping of $f$ to the corresponding solution $u_\epsilon$.
To construct a Green's function, we want to use the Riesz representation theorem: view $R_\epsilon(\cdot)(x)$ as a linear functional on (say) $L^p$, with $x$ fixed. If we know that it is bounded
$$|R_\epsilon(f)(x)|\leq C \|f\|_{L^p}, \qquad \qquad (*)$$
we can write
$$
R_\epsilon(f)(x) = \int G_\epsilon(x, y) f(y) dy
$$
for some function $G_\epsilon$ with $G_\epsilon(x, \cdot) \in L^{\frac{p}{p - 1}}$, the dual space.
We will return to the estimate $(*)$ shortly. The next step would be to take the limit as $\epsilon \rightarrow 0$. This is immediately problematic, as we need to argue that $G_\epsilon(x, y)$ has some compactness, that $R_\epsilon(f)$ converges to a function $R(f)$ which we can make sense of as a solution, and possibly that $R(f)$ is the unique solution with the given data. We are aided somewhat by getting a better estimate, which reads:
$$|R_\epsilon(f)(x) - R_\epsilon(f)(y)|\leq C \|f\|_{L^p}|x - y|^\alpha \qquad \qquad (**)$$
for some $\alpha > 0$. This lets you pass to the limit in $R_\epsilon(f)(\cdot)$ (uniform convergence in $x$) for each fixed $f$, recovering a function $R(f)(\cdot)$ satisfying the same estimates. The Riesz representation theorem can again be used to write $R$ as $R(f)(x) = \int G(x, y) f(y) dy$. The remaining questions, then, are: is $R(f)$ a "meaningful" solution to the PDE, and is it unique?
Unfortunately, it's not clear that $R(f)$ is $C^2$, or has two bounded derivatives: indeed, if the coefficients are only assumed to be bounded, it does not need to be (i.e. there are counterexamples where $R(f)$ fails to have bounded, or $L^p$, or $L^1$, second derivatives; they do exist distributionally in $L^\delta$ for a small $\delta$ by a theorem of F. Lin). We also cannot use distributional notions of solution, as there's no way to integrate by parts on to a test function, you will always have derivatives hitting the coefficients when you try. So possible options are: (1) just stick with $R(f)$ as a solution in the approximating sense that $R_\epsilon(f) \rightarrow R(f)$, (2) talk about solutions in $W^{2, \epsilon}$ which satisfy the PDE a.e., (3) talk about viscosity solutions, which is a concept of weak solution based on the maximum principle. As for uniqueness: it's not known in this generality for any of these classes of solution.
The estimates
The estimates $(*), (**)$ are best possible in $L^\infty$-like spaces. $(*)$ is the Aleksandrov-Backel'man-Pucci estimate, $(**)$ is due to Krylov and Safonov. Both are valid with $p = n$, or with more work for $p > n - \delta$ for a small $\delta$ depending on the ellipticity constants ($n$ is the dimension).
As mentioned, the best "Sobolev space"-type estimate places $P(f) \in W^{2, \delta}$ for some small $\delta$, due to F. Lin.
There are some additional estimates in the $y$ variable on $G$; these are usually of integral type and based on studying the adjoint equation ($L^* u = \partial_{ij} (a_{ij} u)$). The ones I am familiar with are due to P. Bauman, but am not an expert on the topic.
If you assume that $a_{ij}$ are continuous (or have small oscillation), you get a lot more estimates, including ones which can ensure uniqueness.
The books of Gilbarg Trudinger and Caffarelli Cabre are both reasonable references on this topic.
Best Answer
There might be subtleties depending on the exact assumptions that you make, but the general reasoning is as follows.
When $u \in H^1(\Omega)$, its trace on $\partial \Omega$ belongs to $H^{1/2}(\partial\Omega)$. Conversely, if you take a $v \in H^{1/2}(\partial\Omega)$, you can find a $\bar{v} \in H^1(\Omega)$ such that $\bar v = v$ on $\partial \Omega$ and $$ \| \bar{v} \|_{H^1(\Omega)} \leq C \| v \|_{H^{1/2}(\partial\Omega)}. $$ Then you look for $u$ under the form $u = u' + \bar{v}$ so that $u'$ is a solution to $Lu' = f$ with $f = - L\bar{v}$ and $u' = 0$ on $\partial \Omega$. For this problem, you might know that $$ \| u' \|_{H^1(\Omega)} \leq C \| f \|_{H^{-1}(\Omega)} $$ and here you have $$ \|f\|_{H^{-1}(\Omega)} \leq C \| \bar{v} \|_{H^1(\Omega)}. $$ Putting all of this together leads you to the estimate $$ \| u \|_{H^1(\Omega)} \leq C \| v \|_{H^{1/2}(\Omega)}. $$ Now the $H^{1/2}=W^{1/2,2}$ fractional Sobolev norm is bounded above by the $W^{1,\infty}$ norm which you use in your assumption. So, essentially, what you are looking for is indeed valid.
As mentioned above, there are subtleties depending on the exact regularity of $\partial\Omega$ and on the coefficients of $L$.