In the case of X-states there is an analytical solution http://arxiv.org/abs/1501.06900

For matrices in general form you can try maximal-correlation-direction measurement http://arxiv.org/abs/1009.1476

Also, if you could restrict yourself to the case of von Neumann projective measurements you can consider measurement in the basis {U(theta,phi)|0>,  U(theta,phi)|1>},  where U(theta,phi) is unitary rotation matrix of the Bloch sphere. Then the problem reduces to maximization of mutual information, which depends only on two angels.

Thank you Evgeny Kiktenko. My states are not X-states.  I'll try these and get back if I need more discussion.

quantum discord is defined in terms of the quantum mutual information. For a classical case, 

I (A; B) = H (A) + H (B) - H (A,B)

J (A; B) = H (A) - H (A|B)

D=I (A; B)-J (A; B)

For a quantum case, the quantum physics analogy for the three terms are used – S(ρA) the von Neumann entropy, S(ρ) the joint quantum entropy and S(ρA|ρB) the conditional quantum entropy, respectively, for probability density function ρ

I (ρ) = S (ρA) + S (ρB) - S (ρ)

JA (ρ) = S (ρB) - S (ρB|ρA)

The notation J represents the part of the correlations that can be attributed to classical correlations and varies in dependence on the chosen eigenbasis; therefore, in order for the quantum discord to reflect the purely nonclassical correlations independently of basis, it is necessary that J first be maximized over the set of all possible projective measurements onto the eigenbasis

DA (ρ) = I (ρ) -max{MjA} J{MjA} (ρ) = S (ρA) - S(ρ) + \min{MjA} S (ρB |{MjA} )

where {MjA} are the measurement operators. For two qubit states, you can choose the Bell states, Werner state or isotropic state and calculate the quantum discord for this case.