boundary conditions with non-constant coefficients in cell centred finite volume method

• Suppose am solving the heat conduction equation in 1d with Dirichlet boundary conditions. The thermal conductivity $k$ is a non constant function. So $-(k(x)u'(x))' = f(x)$

• The value of $k$ enters the discretization only through its values at the faces of the control volume $k_{i+1/2}, \; k_{i-1/2}$, and not at the cell centres $k_i$.

• But most books and notes available online suggest that we store the value of $k$ at cell centres a priori and use it in the computation.

1. If $k(x)$ has jumps across the faces, I understand one uses the harmonic average. What about control volumes near the boundary ? What is the standard way to approximate $k(x)$ there; harmonic averaging with neighbour is not possible. Can I simply use the cell centre value for these cells ?

I am unable to see how approximating $k(x)$ will eventually affect the order of accuracy of the method in these cases.

My question is about accuracy, not exactly the same as How should non-constant coefficients be treated with finite-volume first order upwind scheme?,

EDIT: When $k(x)$ is smooth, by expanding the terms $k_{i+1/2}, \; k_{i-1/2}$ in a Taylor's series around $k_i$ in the standard 3 point scheme, I see that using the cell centre value directly is atleast $O(h^2)$ or second order accurate. But that would violate the conservation property of the scheme since $k_{i+1/2}^+ \neq k_{i+1/2}^-$. But we can use linear interpolation with $k_{i}, \; k_{i+1}$ to calculate $k_{i+1/2}$ and maintain $O(h^2)$ accuracy.

But in case of discontinuous $k$, I am still puzzled by order of accuracy of the harmonic mean approximation for $k_{i+1/2}, \; k_{i-1/2}$, and how to calculate $k_{i+1/2}, \; k_{i-1/2}$ at the boundary cells.