Warning: what you are setting out to do is achievable, but the calculations are quite bothersome. I will offer a static approximation for the reaction forces, and then the outline of the dynamic "proper" solution. Once you have the reaction forces with respect to a given angle(and some other parameters if you choose the hard way), you can iterate for small timesteps to obtain the solution (like you suggested).
General intro
Let's use the points O, A, C in the image below. Also please notice that I corrected the the normal reaction to point A and the direction of friction. The friction is opposing motion, and the motion is along the ramp.
Finally, if $\phi$ is the angle that the bottom side of the box forms with the horizontal plane, then you can find the coordinates of the point A, C with respect to angle $\phi$, by using the rotation matrix
$$R(\phi) =\begin{bmatrix}\cos(\phi)& -\sin(\phi) & 0 \\ \sin(\phi)&\cos(\phi)& 0
\\0 & 0 &1\end{bmatrix}$$
For example for point A, the coordinate when the box is completely flat with respect to O would be $A_0=[l_p, 0,0]^T$ (from now on, for simplicity, i'll be lazy and skip the formal notation $[]^T$ for simplicity and just use $A_0=[l_p, 0,0]$ for a vector).
This will yield that the coordinates of A wrt to $\phi$ are:
$$A(a) = R(a)A_0 = [l_p\cos(\phi), l_p\sin(\phi),0]$$
Similarly you can calculate the coordinates for C with respect to point O.
$$C(a) = R(\phi)C_0 = \left[\frac{l_p}{2}\cos(\phi)- \frac{h_p}{2}\sin(\phi), \frac{l_p}{2}\sin(\phi) + \frac{h_p}{2}\cos(\phi)\right]$$
Now you can easily write the equations of dynamic equilibrium to calculate reaction forces and accelerations.
$$\sum F_x = m a_{C,x} \rightarrow \mu_c F_{n,c} + \mu_r F_{n,r}\cos(\theta) - F_{n,r}\sin(\theta) =m a_{C,x}$$
$$\sum F_y = m a_{C,y} \rightarrow F_{n,c} + \mu_r F_{n,r}\sin(\theta) + F_{n,r}\cos(\theta) - F_g =m a_{C,y}$$
$$\sum M_O = I \ddot{a} \rightarrow \vec{r}_{c} \times \vec{F}_{g} + \vec{r}_{OA} \times\left(\overrightarrow{\mu_r F_{n,r}} + \vec{F}_{n,r}\right) =I \ddot{\phi} + r_{OC}\times (m \vec{a}_C)$$
Again, I am being very lazy and I:
- denote $\overrightarrow{\mu_r F_{n,r}}$ for the vector of the friction force on the ramp
- use the vector form for a force moment $\vec{r}\times\vec{F}$, which is a shorthand for:
$$\vec{r}\times\vec{F}= \left|\begin{matrix}\vec{i} &\vec{j}&\vec{k}
\\ r_x & r_y & r_z\\ F_x & F_y & F_z \end{matrix}\right|$$
STATIC approximation
At this point, there is no appreciable acceleration, you can get away with setting $a_{C,x},a_{C,y}, \ddot{\phi}$ equal to zero and you just solve the system and estimate the reaction forces.
That will be a lot easier, and you'd be within 10% error (if the assumptions are valid).
If that's not good enough, then proceed to the next section( at your own peril :-) ).
Truly Dynamic
AS you might have noticed you now have 3 equations with (Seemingly) 5 unknowns ($F_{n,c}, F_{n,r}, a_{C,x},a_{C,y}, \ddot{\phi} $). The word Seemingly is stressed because you also need to know $\dot{\phi}$ which doesn't come into play yet but it will in a moment.
So in order to solve the system you'd need at least another 3 equations. What you can use is the relative velocity and acceleration between points O,and A. Basically, the velocity and acceleration at O should only have a horizontal component, while similarly the components velocity and acceleration of point A should comply with the following constraint $\tan\theta=\frac{a_{A,y}}{a_{A,x}}$.
So for example, the relative velocity of A with respect to O is:
$$\begin{cases}
\vec{v}_A = \vec{v}_O + \dot{\phi}\times \vec{r}_{AO}
\end{cases} $$
which - considering $\tan\theta=\frac{v_{A,y}}{v_{A,x}}$ - can be expanded to
$$\begin{cases}
v_{A,x} = v_{O,x} + \left(\dot{\phi}\times \vec{r}_{OA}\right)_x \\
v_{A,y} = v_{O,y} + \left(\dot{\phi}\times \vec{r}_{OA}\right)_y
\end{cases} \rightarrow
\begin{cases}
v_{A,x} = v_{O,x} + \left(\dot{\phi}\times \vec{r}_{OA}\right)_x \\
v_{A,x}\tan(\theta) = \left(\dot{\phi}\times \vec{r}_{OA}\right)_y
\end{cases}
$$
So in the above equation you can estimate the $\dot{\phi}$ if you know the velocity of $v_{O,x}$. However $v_{O,x}$ should be known or can be obtained if you know the kinetic energy of the system at any point (you'd need to calculate first the magnitude of $v_C$, and then use a similar equation to relate $v_A$ to $v_C$. In the simplest case you can start from velocity equal to 0.
The same procedure can be done for the acceleration:
$$\vec{a}_A = \vec{A}_O + \vec{a}_{A|O}= \vec{A}_O + \vec{\dot{\phi}}\times(\vec{\dot{\phi}}\times \vec{r}_{OA}) + \vec{\ddot{\phi}}\times r_{OA}
$$
The point here is to calculate $\ddot{\phi}$, which can be done in a similar way to the $\dot{\phi}$.
If you got $\dot{\phi}$ and $\ddot{\phi}$ then you can relate the acceleration of the center of mass to the acceleration at either points A, O.
Now you should have all available equations to solve the problem. The problem is writing them down in a clear and concise manner and then solve the system. .
