Summary: The Bohr phenomenon for analytic functions of the form \(f(z)=\sum_{n=0}^\infty a_nz^n\), first introduced by H. Bohr [Proc. Lond. Math. Soc. (2) 13, 1–5 (1914; JFM 44.0289.01)], deals with finding the largest radius \(r_f\), \(0<r_f<1\), such that the inequality \(\sum_{n=0}^\infty |a_nz^n|\leq 1\) holds whenever the inequality \(|f(z)|\leq 1\) holds in the unit disk \(\mathbb{D}=\{z\in\mathbb{C}:|z|<1\}\). The exact value of this largest radius known as Bohr radius, which has been established to be \(r_f=1/3\). The Bohr phenomenon [Y. A. Muhanna, Complex Var. Elliptic Equ. 55, No. 11, 1071–1078 (2010; Zbl 1215.46018)] for harmonic functions \(f\) of the form \(f(z)=h(z)+\overline{g(z)}\), where \(h(z)=\sum_{n=0}^\infty a_nz^n\) and \(g(z)=\sum_{n=1}^\infty b_nz^n\) is to find the largest radius \(r_f\), \(0<r_f<1\) such that \[ \sum\limits_{n=1}^\infty(|a_n|+|b_n|)|z|^n\leq d(f(0),\partial f(\mathbb{D})) \] holds for \(|z|\leq r_f\), here \(d(f(0),\partial f(\mathbb{D}))\) denotes the Euclidean distance between \(f(0)\) and the boundary of \(f(\mathbb{D})\). In this paper, we investigate the Bohr radius for several classes of harmonic functions in the unit disk \(\mathbb{D}\).