To describe the anisotropy of sheets a sixth-order polynomial type 3D yield function is proposed. The yield function is constructed as a sum of several components of the Cazacu-Barlat function (2001) which was derived as an extension of the J2-J3 Drucker yield criterion (1949) to orthotropy using the linear transformation of the stress deviator. In this framework of modeling, the convexity of the yield locus is perfectly guaranteed. The model was validated by comparing the numerical predictions of planar anisotropy of r-values and flow stress directionality, as well as the shape of yield loci, with the corresponding experimental data on several types of steel sheets (high r-valued IF steel and SPCE, and high strength steel sheets of 440-980MPa TS grades). For most of steel sheets, the model of the sum of two J2 components, which involve eight anisotropic coefficients, is sufficient for the description of their anisotropies. For the description of the Bauchinger effect and cyclic workhardening characteristic, Yoshida-Uemori kinematic hardening model (2002a, 2002b, 2003) was employed, which includes a new proposal to describe non-saturation type workhardening.