We have measured the temperature and magnetic-field dependences of the sound velocity for one longitudinal and two transverse waves in the low field phase (LFP) and the high field phase (HFP) of nuclear spin ordered bcc solid 3He crystals with a single magnetic domain along the melting curve. From sound velocity measurements for various crystal orientations as a function of the sound propagation direction, we determined the elastic stiffness constants, cij(T,B). In the LFP with tetragonal symmetry for the nuclear spin structure, we extracted six nuclear spin elastic stiffness constants Δcijl(T,0.06 T) from the temperature dependence of the sound velocity at 0.06 T and Δcijl(0.5 mK,B) from the magnetic-field dependence of sound velocity at 0.5 mK. In the HFP with cubic symmetry for the nuclear spin structure, we extracted three Δcijh(T,0.50 T) at 0.50 T and Δcijh(0.5 mK,B) at 0.5 mK. At the first-order magnetic phase transition from the LFP to the HFP at the lower critical field Bc1, large jumps in sound velocities were observed for various crystal directions and we extracted three Δcijtotal|Bc1. Using the thermodynamic relation between Δcij and the change in the internal energy for the exchange interaction in this system, ΔUex(T,B), Δcij are related to the generalized second-order Grüneisen constants ΓijX≡∂2ln X/∂εi∂εj as Δcij(T,B)=ΓijXΔUex(T,B), where X represents some physical quantity which depends on the molar volume and εj is the j-th component of a strain tensor. In the LFP, the Δcijl(T,0.06 T) were proportional to T4, and Δcijl(0.5 mK,B) were proportional to B2. We extracted Γijsl for the spin wave velocity in the LFP, sl, from Δcijl(T,0.06 T) and Γij1/Xl for the inverse susceptibility, 1/χl from Δcijl(0.5 mK,B). In the HFP, Δcijh(T,0.50 T) were proportional to T4 and Δcijh(0.5 mK,ΔB) were proportional to ΔB(≡ B − Bc1). We obtained Γijsh for the spin wave velocity in the HFP, sh, from Δcijh(T,0.50 T) and ΓijBc1 for Bc1 from Δcijh(0.5 mK,ΔB). The values obtained for Γijsl and Γij1/Xl were compared with the Multiple Spin Exchange model (MSE) with three parameters by using analytic expressions for sl and χl. The three-parameter MSE does not agree with the observed Δcij in the LFP.