Complex dynamics of phase changes occur when the alloy solution's temperature suddenly drops below a critical value. The well-known Cahn-Hilliard model shows that a system of fourth-order parabolic partial equations controls this intricate process. However, the Cahn-Hilliard equation with more than four component phases in three dimensions has not been solved to our best knowledge. In this work, the negative chemical potential, namely the first variation of free energy, was convoluted with a sixth-order accurate Laplacian kernel and used as the gradient flow in a projected gradient method. Also, we calculated the Lagrangian multiplies of Gibbs n-simplex phase constraint in a nested loop. Numerical examples illustrate that the proposed method can reveal the nucleation, separation, and growth of grains for alloys with up to 16 component phases in three dimensions. When the number of component phases is large than four, we found small grains are usually dissolved and redeposited onto larger ones. However, the separated phases twist into a highly interconnected structure in the binary and ternary alloys.