Finding the stationary points of a functional

variational formulation
1

f = \frac{K_1}{2}[\frac{\partial\phi}{\partial x}]^2 + \frac{K_2}{2}[\frac{\partial\phi}{\partial z}+q_0]^2

Now I have a functional which is the total energy given by:

F = \displaystyle \int_{0 }^{L} \int_{0 }^{\Lambda}f dx dz

Now I want to find the stationary point
\delta F = 0

The boundary conditions are given as

\phi(x,0) = \frac{\pi x}{\Lambda}

\phi(x,L) = \frac{(2u+1)\pi}{2}

given that u is an integer.

q_0 and L and \Lambda are constants.

I am using finite element analysis to solve for \phi computationally

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

pi =3.141

Step 1: Define the Problem (Geometry and Material Parameters)

L_k = 14e-6  # Length of the domain in the y direction
d = 14e-6 # Length of the domain in the x direction
lam = 60e-6
K1 = Constant(7.49e-12)  # Splay elastic constant
K2 = Constant(3.81e-12)  # Twist elastic constant
q0 = Constant(0)

# Create mesh and define function space
dim_x = 64
dim_y = 64
# Step 2: Define Function Space
mesh = RectangleMesh(Point(0, 0), Point(d, L_k),dim_x-1,dim_y-1)  # Define the rectangular mesh
V = FunctionSpace(mesh, "CG", 1)  # Vector function space for the director field

# Define the trial and test functions
phi = TrialFunction(V)
v = TestFunction(V)




# Boundary condition at y=0 (linear with respect to x)
bx0 = Expression("pi*x[0]/lam", degree=1,lam=lam,pi=pi)
bc_0 = DirichletBC(V, bx0, "near(x[1], 0)")
#Step 7 : Apply Boundary Conditions

# Boundary condition at y=L (takes only odd multiples of pi/2)
u_k = Constant(0)

theta_L = Expression("pi/2*(2*u_k+1)", pi=pi,u_k=u_k, degree=1)
bc_L = DirichletBC(V, theta_L, "near(x[1], 14e-6)")

bcs= [bc_0,bc_L]

L = (K1 / 2) * (grad(phi)[0] ** 2) + (K2 / 2) * (grad(phi)[1] + q0) ** 2
F = L * dx

F_var = derivative(F, phi,v)

problem = NonLinearVariationalProblem(F_var, phi, bcs)

solver = NonLinearVariationalSolver(problem)
solver.solve()

But unfortunately I get an error message like this

Invalid coefficient type for
Traceback (most recent call last):

File ~/opt/anaconda3/envs/fenicsproject/lib/python3.11/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
  exec(code, globals, locals)

File ~/Desktop/LIQUID CRYSTALS/codebase/weak_euler_lagrangian.py:54
  F_var = derivative(F, phi,v)

File ~/opt/anaconda3/envs/fenicsproject/lib/python3.11/site-packages/dolfin/fem/formmanipulations.py:82 in derivative
  return ufl.derivative(form, u, du, coefficient_derivatives)

File ~/opt/anaconda3/envs/fenicsproject/lib/python3.11/site-packages/ufl/formoperators.py:280 in derivative
  coefficients, arguments = _handle_derivative_arguments(form, coefficient,

File ~/opt/anaconda3/envs/fenicsproject/lib/python3.11/site-packages/ufl/formoperators.py:230 in _handle_derivative_arguments
  error("Invalid coefficient type for %s" % ufl_err_str(c))

File ~/opt/anaconda3/envs/fenicsproject/lib/python3.11/site-packages/ufl/log.py:172 in error
  raise self._exception_type(self._format_raw(*message))

UFLException: Invalid coefficient type for

Please can someone help understand the issue with the code?
Thank you and if there is something wrong with the mathematics please correct me on that aspect as well.